-
Notifications
You must be signed in to change notification settings - Fork 0
/
SwanPropvelS.ftn90
335 lines (335 loc) · 13 KB
/
SwanPropvelS.ftn90
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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
subroutine SwanPropvelS ( cad , cas , ux2 , uy2 , &
dep1 , dep2 , cax , cay , &
kwave , cgo , spcsig, iddlow, &
iddtop, ecos , esin , coscos, &
sincos, sinsin, rdx , rdy , &
dhdx , dhdy , dkdx , dkdy , &
duxdx , duxdy , duydx , duydy )
!
! --|-----------------------------------------------------------|--
! | Delft University of Technology |
! | Faculty of Civil Engineering and Geosciences |
! | Environmental Fluid Mechanics Section |
! | P.O. Box 5048, 2600 GA Delft, The Netherlands |
! | |
! | Programmer: Marcel Zijlema |
! --|-----------------------------------------------------------|--
!
!
! SWAN (Simulating WAves Nearshore); a third generation wave model
! Copyright (C) 1993-2024 Delft University of Technology
!
! This program is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with this program. If not, see <http://www.gnu.org/licenses/>.
!
!
! Authors
!
! 40.59: Erick Rogers
! 40.80: Marcel Zijlema
! 41.02: Marcel Zijlema
! 41.06: Gerbrant van Vledder
! 41.07: Marcel Zijlema
! 41.35: Casey Dietrich
! 41.60: Marcel Zijlema
! 41.67: Marcel Zijlema
!
! Updates
!
! 40.80, July 2007: New subroutine
! 40.59, August 2007: muddy bottom included
! 41.02, February 2009: adaption of velocities in case of diffraction
! 41.06, March 2009: add option of limitation of velocity in theta-direction
! 41.07, August 2009: add option of alternative formula for computation of
! wave transport velocity in theta-direction based on
! (x,y)-derivatives of the wave number
! (see Holthuijsen (2007), page 210, footnote 4)
! 41.35, March 2012: add option of limitation on csigma and ctheta
! 41.60, July 2015: more accurate computation of gradients of depth or wave number for turning rate
! 41.67, August 2017: more accurate computation of gradients of ambient currents for transport velocities
!
! Purpose
!
! computes wave transport velocities of energy in spectral space
!
! Modules used
!
use ocpcomm4
use swcomm2
use swcomm3
use swcomm4
use m_diffr
use SwanGriddata
use SwanGridobjects
use SwanCompdata
!
implicit none
!
! Argument variables
!
integer, intent(in) :: iddlow ! minimum direction bin that is propagated within a sweep
integer, intent(in) :: iddtop ! maximum direction bin that is propagated within a sweep
!
real, dimension(MDC,MSC), intent(out) :: cad ! wave transport velocity in theta-direction
real, dimension(MDC,MSC), intent(out) :: cas ! wave transport velocity in sigma-direction
real, dimension(MDC,MSC,ICMAX), intent(in) :: cax ! wave transport velocity in x-direction
real, dimension(MDC,MSC,ICMAX), intent(in) :: cay ! wave transport velocity in y-direction
real, dimension(MSC,ICMAX), intent(in) :: cgo ! group velocity
real, dimension(MDC), intent(in) :: coscos ! help array containing cosine to power 2 of spectral directions
real, dimension(nverts), intent(in) :: dep1 ! water depth at previous time level
real, dimension(nverts), intent(in) :: dep2 ! water depth at current time level
real, intent(in) :: dhdx ! derivative of depth in x-direction
real, intent(in) :: dhdy ! derivative of depth in y-direction
real, dimension(MSC), intent(in) :: dkdx ! derivative of wave number in x-direction
real, dimension(MSC), intent(in) :: dkdy ! derivative of wave number in y-direction
real, intent(in) :: duxdx ! derivative of ux2 to x
real, intent(in) :: duxdy ! derivative of ux2 to y
real, intent(in) :: duydx ! derivative of uy2 to x
real, intent(in) :: duydy ! derivative of uy2 to y
real, dimension(MDC), intent(in) :: ecos ! help array containing cosine of spectral directions
real, dimension(MDC), intent(in) :: esin ! help array containing sine of spectral directions
real, dimension(MSC,ICMAX), intent(in) :: kwave ! wave number
real, dimension(MDC), intent(in) :: sincos ! help array containing sine * cosine of spectral directions
real, dimension(MDC), intent(in) :: sinsin ! help array containing sine to power 2 of spectral directions
real, dimension(nverts), intent(in) :: ux2 ! ambient velocity in x-direction at current time level
real, dimension(nverts), intent(in) :: uy2 ! ambient velocity in y-direction at current time level
real, dimension(2), intent(in) :: rdx ! first component of contravariant base vector rdx(b) = a^(b)_1
real, dimension(2), intent(in) :: rdy ! second component of contravariant base vector rdy(b) = a^(b)_2
real, dimension(MSC), intent(in) :: spcsig ! relative frequency bins
!
! Local variables
!
integer :: icell ! cell index
integer :: id ! loop counter over direction bins
integer :: iddum ! counter in directional space
integer, save :: ient = 0 ! number of entries in this subroutine
integer :: is ! loop counter over frequency bins
integer :: iv1 ! first index in computational stencil
integer :: iv2 ! second index in computational stencil
integer :: iv3 ! third index in computational stencil
!
real :: alpha ! upper limit of CFL restriction
real, dimension(3) :: cd ! coefficients for computing cad
real, dimension(10) :: cs ! coefficients for computing cas
real, dimension(3) :: dloc ! local depth at vertices
real :: fac ! a factor
real :: fac2 ! another factor
real :: frlim ! frequency range in which limit on velocity in theta-direction is applied
real :: kd ! help variable, wave number times water depth
real :: pp ! power of the frequency dependent limiter on refraction
!
type(celltype), dimension(:), pointer :: cell ! datastructure for cells with their attributes
type(verttype), dimension(:), pointer :: vert ! datastructure for vertices with their attributes
!
! Structure
!
! Description of the pseudo code
!
! Source text
!
if (ltrace) call strace (ient,'SwanPropvelS')
!
! point to vertex and cell objects
!
vert => gridobject%vert_grid
cell => gridobject%cell_grid
!
! set velocities to zero
!
cad = 0.
cas = 0.
!
! if no frequency shift and no refraction, return
!
if ( (ITFRE == 0 .or. (.not.DYNDEP .and. ICUR == 0)) .and. IREFR == 0 ) return
!
! initialize coefficients
!
cd = 0.
cs = 0.
!
iv1 = vs(1)
iv2 = vs(2)
iv3 = vs(3)
!
dloc(1) = dep2(iv1)
dloc(2) = dep2(iv2)
dloc(3) = dep2(iv3)
!
! if at least one vertex is dry, return
!
if ( dloc(1) <= DEPMIN .or. dloc(2) <= DEPMIN .or. dloc(3) <= DEPMIN ) return
!
! compute wave transport velocity in sigma-direction
!
if ( ITFRE /= 0 .and. (DYNDEP .or. ICUR /= 0) ) then
!
! compute time derivative of water depth
!
if ( DYNDEP ) cs(2) = ( dep2(iv1) - dep1(iv1) ) * RDTIM
!
! compute coefficients depending on ambient currents
!
if ( ICUR /= 0 ) then
!
cs(3) = ux2(iv1) * dhdx
cs(4) = uy2(iv1) * dhdy
!
endif
!
do is = 1, MSC
!
! compute frequency-dependent coefficients
!
kd = min(30.,kwave(is,1) * dep2(iv1))
!
cs(1) = kwave(is,1) * spcsig(is) / sinh (2.* kd)
cs(5) = -cgo(is,1) * kwave(is,1)
if ( IDIFFR /= 0 ) cs(5) = cs(5)*DIFPARAM(iv1)
!
cs( 6) = cs(1) * cs(2)
cs( 7) = cs(1) * (cs(3)+cs(4))
cs( 8) = cs(5) * duxdx
cs( 9) = cs(5) * (duxdy+duydx)
cs(10) = cs(5) * duydy
!
do iddum = iddlow-1, iddtop+1
id = mod ( iddum - 1 + MDC , MDC ) + 1
!
cas(id,is) = cs(6)
if ( ICUR /= 0 ) cas(id,is) = cas(id,is) + cs(7) + coscos(id)*cs(8) + sincos(id)*cs(9) + sinsin(id)*cs(10)
!
enddo
enddo
!
! limit velocity using Courant number
!
if ( int(PNUMS(33)) == 1 ) then
!
alpha = PNUMS(34)
!
do is = 1, MSC
!
fac2 = alpha * FRINTF * spcsig(is)
!
do iddum = iddlow-1, iddtop+1
id = mod ( iddum - 1 + MDC , MDC ) + 1
!
fac = fac2 * ( abs((rdx(1)+rdx(2))*cax(id,is,1)) + abs((rdy(1)+rdy(2))*cay(id,is,1)) )
!
if ( abs(cas(id,is)) > fac ) cas(id,is) = cas(id,is) * fac / abs(cas(id,is))
!
enddo
!
enddo
!
endif
!
endif
!
! compute wave transport velocity in theta-direction
!
if ( IREFR /= 0 ) then
!
do is = 1, MSC
!
! compute frequency-dependent coefficients
!
if ( int(PNUMS(32)) == 0 ) then
!
kd = min(30.,kwave(is,1) * dep2(iv1))
!
cd(1) = spcsig(is) / sinh (2.* kd)
!
cd(2) = cd(1) * dhdx
cd(3) = cd(1) * dhdy
!
else
!
cd(1) = -cgo(is,1) / kwave(is,1)
!
cd(2) = cd(1) * dkdx(is)
cd(3) = cd(1) * dkdy(is)
!
endif
!
do iddum = iddlow-1, iddtop+1
id = mod ( iddum - 1 + MDC , MDC ) + 1
!
cad(id,is) = esin(id)*cd(2) - ecos(id)*cd(3)
if ( IDIFFR /= 0 ) cad(id,is) = cad(id,is)*DIFPARAM(iv1) - DIFPARDX(iv1)*cgo(is,1)*esin(id) + DIFPARDY(iv1)*cgo(is,1)*ecos(id)
if ( ICUR /= 0 ) cad(id,is) = cad(id,is) + sincos(id)*(duxdx-duydy) + sinsin(id)*duydx - coscos(id)*duxdy
!
enddo
!
enddo
!
! adapt velocity in case of spherical coordinates
!
if ( KSPHER > 0 ) then
!
cd(1) = tan(DEGRAD*(vert(iv1)%attr(VERTY) + YOFFS)) / REARTH
!
do id = 1, MDC
cd(2) = cd(1) * ecos(id)
do is = 1, MSC
cad(id,is) = cad(id,is) - cd(2)*(cax(id,is,1)*ecos(id) + cay(id,is,1)*esin(id))
enddo
enddo
!
endif
!
! limit velocity in some frequency range if requested
!
if ( int(PNUMS(29)) == 1 ) then
!
frlim = PI2*PNUMS(26)
pp = PNUMS(27)
!
do is = 1, MSC
!
fac = min(1.,(spcsig(is)/frlim)**pp)
!
do id = 1, MDC
cad(id,is) = fac*cad(id,is)
enddo
enddo
!
endif
!
! limit velocity using Courant number
!
if ( int(PNUMS(35)) == 1 ) then
!
alpha = PNUMS(36)
!
fac2 = alpha * DDIR
!
do is = 1, MSC
!
do iddum = iddlow-1, iddtop+1
id = mod ( iddum - 1 + MDC , MDC ) + 1
!
fac = fac2 * ( abs((rdx(1)+rdx(2))*cax(id,is,1)) + abs((rdy(1)+rdy(2))*cay(id,is,1)) )
!
if ( abs(cad(id,is)) > fac ) cad(id,is) = cad(id,is) * fac / abs(cad(id,is))
!
enddo
!
enddo
!
endif
!
endif
!
end subroutine SwanPropvelS