-
Notifications
You must be signed in to change notification settings - Fork 0
/
SeaWaves_multipleheight.m
83 lines (56 loc) · 2.07 KB
/
SeaWaves_multipleheight.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
function [Umi,Tpseai,HS,F,PWTOT,QsWslope_seaTOT]=SeaWaves(A,angle,hwSea_lim,range,wind,VEG,ndir,Nhseawave,dx,z,msl,tempdeltaMSL,ws1,fTide,extraHseaimposed,addextrafetch,extrafetch);
rhos=2650;ss=1.5728;
z=z-tempdeltaMSL;%temporary shift to change MSL at ever tide, trick!!!
[N,M]=size(A);
Lbasin=0;%1000/dx;
Fetchlim=0;%max(50,dx*2);%dx*2;%600;%dx*2*10;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
PWTOT=0;
QsWslope_seaTOT=0;
for i=1:Nhseawave
dHW=max(0,msl+range/2-z);%water depth at MHW
h=max(0,dHW-range*(i-1)/(Nhseawave-1));
MASK=0*A+1;MASK(h<=hwSea_lim | A==0 | VEG==1)=0;
%The standard way
if addextrafetch==0
F=calculatefetch(MASK,ndir,dx,angle);
elseif addextrafetch==1
F=calculatefetchWITHEXTRAS(MASK,ndir,dx,angle,extrafetch,Lbasin,MASK);
end
Fo=F;F(Fo<=Fetchlim)=0;
%diffuse the fetch field
alphadiffusefetch=1; %messo 10 for the VCR wave validation 10;%0;%%%QUESTO ERA 1 FINO AD APRILE 23 2018!!!!!
F=diffusefetchPERPEND(MASK,F,alphadiffusefetch,dx,angle);
F(Fo<=Fetchlim | MASK==0)=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=find(Fo>Fetchlim & h>hwSea_lim & F>0 & MASK==1);%h>dlo & %a=find(Fo>dx*2);%h>dlo & %a=find(h>dlo);
D=h(a);
Ff=F(a);
TP=0*h;HS=0*h;
[Hs,Tp]=YeV(Ff,wind,D);%[Hs,Tp]=YeV(Ff,wind,min(3,D)); %TRUCCO PER EVITARE LARGE WAVES IN CHANELS
HS(a)=Hs;TP(a)=Tp;TP(TP==0)=1;
if extraHseaimposed==1
HS=getextraHsea(HS,h,MASK,dx,N,M);
end
HS(h<hwSea_lim)=0;
TP(h<hwSea_lim)=2;
kwave=0*h;kk=wavek(1./TP(a),h(a));kwave(a)=kk;
kwave(kwave==0 | h<hwSea_lim)=1;
Um=pi*HS./(TP.*sinh(kwave.*h));
cg=(2*pi./kwave./TP)*0.5.*(1+2*kwave.*h./(sinh(2*kwave.*h)));
PW=cg*1030*9.8.*HS.^2/16;
Um(h<hwSea_lim)=0;
PW(h<hwSea_lim)=0;
Um(MASK==0)=0;
PW(MASK==0)=0;
Umi(:,:,i)=Um;
Tpseai(:,:,i)=TP;
PWTOT=PWTOT+PW;
QsWslope_sea=WaveSedimentTransport(HS,h,kwave,rhos,N,M,TP,dx,ss,ws1,hwSea_lim,NaN);
QsWslope_sea(HS==0)=0;
QsWslope_seaTOT=QsWslope_seaTOT+QsWslope_sea;
%figure;imagesc(HS);caxis([0 0.5]);pause
%figure;imagesc(F);caxis([0 10000]);pause
end
PWTOT=PWTOT/Nhseawave;
QsWslope_seaTOT=QsWslope_seaTOT/Nhseawave;