diff --git a/tools/WW3_mesh_read.m b/tools/WW3_mesh_read.m
new file mode 100755
index 000000000..8bc80df28
--- /dev/null
+++ b/tools/WW3_mesh_read.m
@@ -0,0 +1,91 @@
+function [tri,x,y,z] = WW3_mesh_write(file,plott)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This function reads WW3 grid including nodes              %
+% (longitude,latitude,depth) and element connections        %
+%  (triangles)                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%      Ali Abdolali August 2018 ali.abdolali@noaa.gov       %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%if plot==1, then plot the mesh
+
+% no of nodes is mentioned in 5th row and first column
+
+N_n      = dlmread(file,'',[5-1 1-1 5-1 1-1]);
+N_e      = dlmread(file,'',[7+N_n 0 7+N_n 0]);
+
+node_id     = dlmread(file,'',[5 0 4+N_n 0]);
+element_id     = dlmread(file,'',[8+N_n 0 8+N_e+N_n-1 0]);
+
+nodes       = dlmread(file,'',[5 1 4+N_n 3]);
+TRI    = dlmread(file,'',[8+N_n 0 8+N_e+N_n-1 8]);
+
+%------- 2D Geometry
+
+two_d_nodes = nodes(:,1:2);
+elem_type   = TRI(:,3);
+
+%--- find the starting indices of 2D elements
+two_ind = 1;
+for i = 1:N_e
+    if(elem_type(i) ~= 2)
+        two_ind = two_ind+1;
+    end
+end
+%----------------------------------------------
+
+%two_d_elements(1:N_e-two_ind,1:3) = 0;
+   two_d_elements(:,1:3) = TRI(N_e-two_ind+2:N_e,7:9);
+   
+   tri=two_d_elements;
+   x=nodes(:,1);
+   y=nodes(:,2);
+   z=nodes(:,3);
+   
+   
+   %%
+   
+   if plott==1;
+width=880;  % Width of figure for movie [pixels]
+height=700;  % Height of figure of movie [pixels]
+left=700;     % Left margin between figure and screen edge [pixels]
+bottom=200;  % Bottom margin between figure and screen edge [pixels]
+
+figure
+set(gcf,'Position', [left bottom width height])
+
+
+cmap = colormap;
+
+trisurf(tri,x,y,z);
+shading interp
+view(2)
+hcb=colorbar
+caxis([min(z) max(z)])
+title(hcb,'[m]','fontsize',14);
+colormap(jet)
+hold on
+
+hold on
+ 
+sz=40;
+
+ 
+xlabel('Longitude ^o ');
+ylabel('Latitude ^o ');
+
+title(['depth (m)'],'fontsize',15)
+
+
+
+ylim([min(y)-0.1*((max(y)-min(y))) max(y)+0.1*((max(y)-min(y)))]);
+xlim([min(x)-0.1*((max(x)-min(x))) max(x)+0.1*((max(x)-min(x)))]);
+
+grid off
+box on
+
+
+ print('-dpng',[file,'.png']);
+   end
+ 
\ No newline at end of file
diff --git a/tools/WW3_mesh_write.m b/tools/WW3_mesh_write.m
new file mode 100755
index 000000000..4927b2126
--- /dev/null
+++ b/tools/WW3_mesh_write.m
@@ -0,0 +1,89 @@
+function [filename]=WW3_mesh_write(tri,x,y,h,OB_ID,filename,plott)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This function writes WW3 grid including nodes             %
+% (longitude,latitude,depth), open bounday nodes            %
+% and element connections (triangles)                       %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%      Ali Abdolali August 2018 ali.abdolali@noaa.gov       %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% tri: triangle (nelem,3);
+% x: longitude (nnode,1);
+% y: laritude (nnode,1);
+% h: depth (nnode,1);
+% OB_ID: Open boundary nodes ID
+% filename: mesh name
+% Plot it, if plott=1
+
+node(:,1)=(1:length(x));
+node(:,2)=x;
+node(:,3)=y;
+node(:,4)=h;
+
+fileID = fopen(filename,'w');
+fprintf(fileID,'%s\n', '$MeshFormat');
+fprintf(fileID,'%s\n', '2 0 8');
+fprintf(fileID,'%s\n', '$EndMeshFormat');
+fprintf(fileID,'%s\n', '$Nodes');
+fprintf(fileID,'%d\n', length(node(:,1)));
+
+for i=1:length(node(:,1))
+    fprintf(fileID,['%d %s %5.5f %s %5.5f %s %5.5f\n'], node(i,1),'', node(i,2),'',node(i,3),'',node(i,4));
+end
+fprintf(fileID,'%s\n', '$EndNodes');
+fprintf(fileID,'%s\n', '$Elements');
+fprintf(fileID,'%d\n', length(tri(:,1))+length(OB_ID));
+m=0;
+for i=1:length(OB_ID)
+    m=m+1;
+    fprintf(fileID,['%d %s %d %s %d %s %d %s %d %s %d\n'], m,'',15,'',2,'',0,'',0,'',OB_ID(i));
+end
+
+for i=1:length(tri(:,1))
+    m=m+1;
+    fprintf(fileID,['%d %s %d %s %d %s %d %s %d %s %d %s %d %s %d %s %d\n'], m,'',2,'',3,'',0,'',i,'',0,'',tri(i,1),'',tri(i,2),'',tri(i,3));
+end
+fprintf(fileID,'%s', '$EndElements');
+fclose(fileID);
+
+if plott==1
+width=880;  % Width of figure for movie [pixels]
+height=700;  % Height of figure of movie [pixels]
+left=700;     % Left margin between figure and screen edge [pixels]
+bottom=200;  % Bottom margin between figure and screen edge [pixels]
+ 
+
+
+    figure
+set(gcf,'Position', [left bottom width height])
+
+
+
+trisurf(tri,node(:,2),node(:,3),node(:,4));
+shading interp
+
+view(2);
+axis equal
+hold on
+p1=scatter(node(OB_ID,2),node(OB_ID,3),'xk')
+hCbar=colorbar
+colormap(jet)
+legend([p1],'Open Boundary Nodes')
+xlim([min(node(:,2))-(max(node(:,2))-min(node(:,2)))/10 max(node(:,2))+(max(node(:,2))-min(node(:,2)))/10])
+ylim([min(node(:,3))-(max(node(:,3))-min(node(:,3)))/10 max(node(:,3))+(max(node(:,3))-min(node(:,3)))/10])
+
+if min(node(:,4))~=max(node(:,4))
+caxis([min(node(:,4)) max(node(:,4))])
+end
+
+hold on
+xlabel('Longitude ^{\circ} ');
+ylabel('Latitude ^{\circ} ');
+axis equal
+box on
+grid off
+axis on
+title(['WW3 Grid - depth [m]'],'fontsize',15)
+
+print(gcf,'-dpng',[filename,'.png'],'-r900');
+end
diff --git a/tools/readfort14.m b/tools/readfort14.m
new file mode 100755
index 000000000..72826c56f
--- /dev/null
+++ b/tools/readfort14.m
@@ -0,0 +1,171 @@
+function [EToV,VX,B,opedat,boudat,title] = readfort14( finame, read_bou )
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This function reads fort.14 file including nodes          %
+% (longitude,latitude,depth), open bounday nodes            %
+% and element connections (triangles)                       %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%      Ali Abdolali August 2018 ali.abdolali@noaa.gov       %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+disp('read adcirc fort.14') ;
+
+tic
+if ( nargin == 0 )
+    finputname = 'fort.14'
+else
+    finputname = finame ;
+end
+
+fid = fopen(finputname) ;
+
+agrid = fgetl(fid) ;
+disp(agrid) ;
+title = agrid ;
+
+N = fscanf(fid,'%g %g',2) ;
+
+% Val = zeros(N(2),4) ;
+
+%
+% for i = 1: N(2)
+%    Val(i,1:4) = fscanf(fid,'%d %g %g %g \n', 4 ) ;
+% end
+%
+% Nov 15, 2012, improve reading efficiency
+Val = fscanf(fid,'%d %g %g %g \n', [4 N(2)])' ;
+
+%Val = Val(iv,:) ;
+
+% idx = zeros(N(1),5) ;
+%
+% for i = 1: N(1)
+%    idx(i,1:5) = fscanf(fid,'%d %d %d %d %d \n', 5) ;
+% end
+%
+% Nov 15, 2012, improve reading efficient
+idx = fscanf(fid,'%d %d %d %d %d \n', [5 N(1)])' ;
+
+%idx = idx(iv,:) ;
+
+% VX = zeros(N(2),2) ;
+% B = zeros(N(2),2) ;
+% EToV = zeros(N(1),3) ;
+
+% Arrange it to a Nodal DG input
+VX = Val(:,2:3) ;
+B  = Val(:,4) ;
+EToV = idx(:,3:5) ;
+
+if(read_bou)
+    % Read in boundary
+    % Open boundary
+    msgline = fgetl(fid) ;
+    nope = sscanf(msgline,'%d %*s') ;
+    
+    msgline = fgetl(fid) ;
+    neta = sscanf(msgline,'%d %*s') ;
+    
+    nvdll = zeros(nope,1) ;
+    ibtypee = zeros(nope,1) ;
+    nbdv = zeros(neta,nope) ;
+    % tic
+    for i = 1: nope
+        msgline = fgetl(fid) ;
+        
+        [varg] = sscanf(msgline,'%d %*s \n') ;
+        nvdll(i) = varg ;
+        ibtypee(i) = 0 ;
+        
+        %
+        % for k = 1: nvdll(i)
+        %
+        %    % % nbdv(k,i) = fscanf(fid,'%d \n')
+        %    % % fscanf(fid,'%d \n')
+        %    msgline = fgetl(fid) ;
+        %
+        %    nbdv(k,i) = str2num(msgline) ;
+        % end
+        %
+        % Nov 25, 2012, improve reading efficiency
+        nbdv(1:nvdll(i),i) = fscanf(fid,'%d \n', nvdll(i) ) ;
+    end
+    % toc
+    
+    % ocean boundary
+    opedat.nope = nope ;
+    opedat.neta = neta ;
+    opedat.nvdll = nvdll ;
+    opedat.ibtypee = ibtypee ;
+    opedat.nbdv = nbdv ;
+    
+    % land boundary
+    msgline = fgetl(fid) ;
+    nbou = sscanf(msgline,'%d %*s') ;
+    
+    msgline = fgetl(fid) ;
+    nvel = sscanf(msgline,'%d %*s') ;
+    
+    nvell = zeros(nbou,1) ;
+    ibtype = zeros(nbou,1) ;
+    nbvv = sparse(nvel,nbou) ;
+    ibconn = sparse(nvel,nbou) ;
+    barinht = sparse(nvel,nbou) ;
+    barincfsb = sparse(nvel,nbou) ;
+    barincfsp = sparse(nvel,nbou) ;
+    
+    % tic
+    for i = 1: nbou
+        msgline = fgetl(fid) ;
+        
+        [varg,~] = sscanf(msgline,'%d %d %*s \n') ;
+        nvell(i) = varg(1) ;
+        ibtype(i) = varg(2) ;
+        
+        switch ( ibtype(i) )
+            case {0,1,2,10,11,12,20,21,22,30,60,61,101,52}
+                % Nov 15, 2012, improve reading efficiency
+                nbvv(1:nvell(i),i) = fscanf(fid,'%d \n', nvell(i) ) ;
+            case  {3, 13, 23}
+                disp('3 13 23')
+                val = fscanf(fid,'%g %g %g \n', [3 nvell(i)] )  ;
+                nbvv(1:nvell(i),i) = val(1,:) ;
+            case  {4, 24}
+                %disp('4 24')
+                val = fscanf(fid,'%g %g %g %g %g \n', [5 nvell(i)] )  ;
+                nbvv(1:nvell(i),i) = val(1,:) ;
+                ibconn(1:nvell(i),i) = val(2,:) ;
+                barinht(1:nvell(i),i) = val(3,:) ;
+                barincfsb(1:nvell(i),i) = val(4,:) ;
+                barincfsp(1:nvell(i),i) = val(5,:) ;
+            case  {5, 25}
+                %disp('5 25')
+                val = fscanf(fid,'%g % g %g %g %g %g %g %g \n', [8 nvell(i)] ) ;
+                nbvv(1:nvell(i),i) = val(1,:) ;
+                %otherwise
+                %    msgline = fgetl(fid) ;
+        end
+    end
+    % toc
+    
+    % land boundary
+    boudat.nbou = nbou ;
+    boudat.nvel = nvel ;
+    boudat.nvell = nvell ;
+    boudat.ibtype = ibtype ;
+    boudat.nbvv = nbvv ;
+    
+    if ( sum(ibtype == 24) > 0 ||  sum(ibtype == 4) > 0 )
+        boudat.ibconn = ibconn ;
+        boudat.barinht = barinht ;
+        boudat.barincfsb = barincfsb ;
+        boudat.barincfsp = barincfsp ;
+    end
+else
+    opedat = []; 
+    boudat  = []; 
+end
+
+fclose(fid) ;
+
+return
+