clear all;close all;
% jcwarner 19August2014
% a list of procedures and notes
% used to create the Sandy app files
% for roms and swan. The wrf files
% are assumed to already be created.
%
% cd Projects/Sandy
% This is a place to create all the necessary files.
%
% Step 0.  We are drawing heavily on the work of the scientists at
% NCAR who already created the WRF input files. These were provided
% by Dave Gill and Cindy Bruyere. They provided a link to a WRF
% tutorial here:
%
%If you are interested in the WRF-only processing (allowing you to define the atmospheric domain), please take a look at the following web page:
%http://www2.mmm.ucar.edu/wrf/OnLineTutorial/compilation_tutorial.php
%
%
% Step 1:  ROMS grid
%
% As a first cut, you can use the wrf grid. 
% Load the wrf grid to get coastline data

% netcdf_load('D:\MATLAB\matlab_ROM\coawst\sandy\testforaqua\wrfinput_d01')
% figure
% pcolorjw(XLONG,XLAT,double(1-LANDMASK))
% hold on

% pick 4 corners for roms parent grid
Istr=96; Iend=156;
Jstr=5 ; Jend=49;

%step 3: pick number of points in the grid
% numx=503;
% numy=333;

numx=604;
numy=500;

% Make a matrix of the lons and lats
% You should not need to change this.
% start in lower left and go clockwise
% corner_lon=double([XLONG(Istr,Jstr) XLONG(Istr,Jend) XLONG(Iend,Jend) XLONG(Iend,Jstr) ]);
% corner_lat=double([XLAT(Istr,Jstr) XLAT(Istr,Jend) XLAT(Iend,Jend) XLAT(Iend,Jstr) ]);
corner_lon=double([Istr Istr Iend Iend ]);
corner_lat=double([Jstr Jend Jend Jstr ]);

x=[1 numx; 1 numx]; y=[1 1; numy numy];
z=[corner_lon(1) corner_lon(2) corner_lon(4) corner_lon(3)];
F = TriScatteredInterp(x(:),y(:),z(:));
[X,Y]=meshgrid(1:numx,1:numy);
lon=F(X,Y).';
%
x=[1 numx; 1 numx]; y=[1 1; numy numy];
z=[corner_lat(1) corner_lat(2) corner_lat(4) corner_lat(3)];
F = TriScatteredInterp(x(:),y(:),z(:));
lat=F(X,Y).';
% plot(lon,lat,'k-')
% plot(lon',lat','k-')
%
% Call generic grid creation
%
roms_grid='roms_grid_tccip10wrf.nc';
rho.lat=lat;
rho.lon=lon;
rho.depth=zeros(size(rho.lon))+100;  % for now just make zeros
rho.mask=zeros(size(rho.lon));       % for now just make zeros
spherical='T';
%projection='stereographic';
projection='mercator';
save temp_jcw33.mat rho spherical projection
eval(['mat2roms_mw(''temp_jcw33.mat'',''',roms_grid,''');'])
!del temp_jcw33.mat
%User needs to edit roms vbaariables
disp(['     '])
disp(['Created roms grid -->   ',roms_grid])
disp(['     '])
disp(['You need to edit depth in ',roms_grid])
disp(['     '])
%
% Step 2: Masking
%
% base this on the WRF mask.
% F = TriScatteredInterp(double(XLONG(:)),double(XLAT(:)), ...
% %                        double(1-LANDMASK(:)),'nearest');
%  roms_mask=F(lon,lat);
% % figure
% % pcolorjw(lon,lat,roms_mask)
% 
% water = double(roms_mask);
% u_mask = water(1:end-1,:) & water(2:end,:);
% v_mask= water(:,1:end-1) & water(:,2:end);
% psi_mask= water(1:end-1,1:end-1) & water(1:end-1,2:end) & water(2:end,1:end-1) & water(2:end,2:end);
% ncwrite('roms_grid_d01.nc','mask_rho',roms_mask);
% ncwrite('roms_grid_d01.nc','mask_u',double(u_mask));
% ncwrite('roms_grid_d01.nc','mask_v',double(v_mask));
% ncwrite('roms_grid_d01.nc','mask_psi',double(psi_mask));

%
% Step 3: Bathymetry
spherical='T';
%projection='lambert conformal conic';
projection='mercator';
ROMSTOOLS_dir = '';
hmin = 50;hmax_coast = 500;rtarget = 0.4;
n_filter_deep_topo=4;
n_filter_final=2;
topofile = [ROMSTOOLS_dir,'D:\MATLAB\matlab_ROM\Roms_tools/Topo/etopo2.nc'];
grdname='D:\MATLAB\matlab_ROM\coawst\sandy\testforaqua\roms_grid_tccip10wrf.nc';
h=add_topo(grdname,topofile);
maskr=h>0;
maskr=process_mask(maskr);
[masku,maskv,maskp]=uvp_masks(maskr');
nc=netcdf(grdname,'write');
nc{'h'}(:)=h;
nc{'mask_u'}(:)=masku';
nc{'mask_v'}(:)=maskv';
nc{'mask_psi'}(:)=maskp';
nc{'mask_rho'}(:)=maskr;
 %nc{'angle'}(:)=angle;
close(nc);
%
%  Smooth the topography
%
nc=netcdf(grdname,'write');
h=nc{'h'}(:);
 maskr=nc{'mask_rho'}(:);
%
h=smoothgrid(h,maskr,hmin,hmax_coast,...
             rtarget,n_filter_deep_topo,n_filter_final);
ncwrite(grdname,'mask_rho',maskr');
% ncwrite(grdname,'mask_u',double(u_mask));
% ncwrite(grdname,'mask_v',double(v_mask));
% ncwrite(grdname,'mask_psi',double(psi_mask));
% %
%  Write it down
%
disp(' ')
disp(' Write it down...')
nc{'h'}(:)=h;
close(nc);
netcdf_load(roms_grid)
% figure
% pcolorjw(lon_rho,lat_rho,h)
% hnew(isnan(hnew))=50;
% ncwrite(roms_grid,'h',hnew);
%
goon=0;
if goon==1


%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  coupling to SWAN
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Step 15:  create swan grids: run roms2swan twice
roms2swan('roms_grid_d01.nc')


COAWST\Tools\mfiles\swan_forc\ww3_swan_input.m

end







