WAVEWATCH III  beta 0.0.1
Functions/Subroutines
w3sbt8md Module Reference

Contains routines for computing dissipation by viscous fluid mud using Dalrymple and Liu (1978) "Thin Model". More...

Functions/Subroutines

subroutine w3sbt8 (AC, H_WDEPTH, S, D, IX, IY)
 Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978). More...
 
subroutine csinh (C, CS)
 Complex hyperbolic sin (sinh). More...
 
subroutine ccosh (C, CC)
 Complex hyperbolic cos (cosh). More...
 

Detailed Description

Contains routines for computing dissipation by viscous fluid mud using Dalrymple and Liu (1978) "Thin Model".

Author
M. Orzech
W. E. Rogers
Date
21-Nov-2013
Copyright
Copyright 2009-2022 National Weather Service (NWS), National Oceanic and Atmospheric Administration. All rights reserved. WAVEWATCH III is a trademark of the NWS. No unauthorized use without permission.

Function/Subroutine Documentation

◆ ccosh()

subroutine w3sbt8md::ccosh ( complex, intent(in)  C,
complex, intent(out)  CC 
)

Complex hyperbolic cos (cosh).

Parameters
[in]C
[out]CC
Author
NA
Date
NA

Definition at line 524 of file w3sbt8md.F90.

524  COMPLEX, INTENT(IN) :: C
525  COMPLEX, INTENT(OUT) :: CC
526  x = real(c)
527  y = aimag(c)
528  cc = cmplx(cosh(x) * cos(y), sin(y) * sinh(x))
529  RETURN

Referenced by w3sbt8().

◆ csinh()

subroutine w3sbt8md::csinh ( complex, intent(in)  C,
complex, intent(out)  CS 
)

Complex hyperbolic sin (sinh).

Parameters
[in]C
[out]CS
Author
NA
Date
NA

Definition at line 505 of file w3sbt8md.F90.

505  COMPLEX, INTENT(IN) :: C
506  COMPLEX, INTENT(OUT) :: CS
507  x = real(c)
508  y = aimag(c)
509  cs = cmplx(sinh(x) * cos(y), sin(y) * cosh(x))
510  RETURN

Referenced by w3sbt8().

◆ w3sbt8()

subroutine w3sbt8md::w3sbt8 ( real, dimension(nspec), intent(in)  AC,
real, intent(in)  H_WDEPTH,
real, dimension(nspec), intent(out)  S,
real, dimension(nspec), intent(out)  D,
integer, intent(in)  IX,
integer, intent(in)  IY 
)

Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978).

"Thin Model" (adapted from Erick Rogers code by Mark Orzech, NRL).

Parameters
[in]ACAction density spectrum (1-D).
[in]H_WDEPTHMean water depth.
[out]SSource term (1-D version).
[out]DDiagonal term of derivative (1-D version).
[in]IX
[in]IY
Author
M. Orzech
W. E. Rogers
Date
21-Nov-2013

Definition at line 112 of file w3sbt8md.F90.

112  !/
113  !/ +-----------------------------------+
114  !/ | WAVEWATCH III NOAA |
115  !/ | M. Orzech NRL |
116  !/ | W. E. Rogers NRL |
117  !/ | FORTRAN 90 |
118  !/ | Last update : 21-Nov-2013 |
119  !/ +-----------------------------------+
120  !/
121  !/ 20-Dec-2004 : Origination. ( version 3.06
122  !/ 23-Jun-2006 : Formatted for submitting code for ( version 3.09 )
123  !/ inclusion in WAVEWATCH III.
124  !/
125  ! 1. Purpose :
126  !
127  ! Compute dissipation by viscous fluid mud using Dalrymple and Liu (1978)
128  ! "Thin Model" (adapted from Erick Rogers code by Mark Orzech, NRL).
129  !
130  ! 2. Method :
131  !
132  ! 3. Parameters :
133  !
134  ! Parameter list
135  ! ----------------------------------------------------------------
136  ! AC R.A. I Action density spectrum (1-D)
137  ! H_WDEPTH Real I Mean water depth.
138  ! S R.A. O Source term (1-D version).
139  ! D R.A. O Diagonal term of derivative (1-D version).
140  ! ----------------------------------------------------------------
141  !
142  ! 4. Subroutines used :
143  !
144  ! Name Type Module Description
145  ! ----------------------------------------------------------------
146  ! STRACE Subr. W3SERVMD Subroutine tracing.
147  ! ----------------------------------------------------------------
148  !
149  ! 5. Called by :
150  !
151  ! Name Type Module Description
152  ! ----------------------------------------------------------------
153  ! W3SRCE Subr. W3SRCEMD Source term integration.
154  ! W3EXPO Subr. N/A Point output post-processor.
155  ! GXEXPO Subr. N/A GrADS point output post-processor.
156  ! ----------------------------------------------------------------
157  !
158  ! 6. Error messages :
159  !
160  ! None.
161  !
162  ! 7. Remarks :
163  !
164  ! Cg_mud calculation could be improved by using dsigma/dk instead
165  ! of n*C. The latter is a "naive" method and its accuracy has
166  ! not been confirmed.
167  !
168  ! 8. Structure :
169  !
170  ! See source code.
171  !
172  ! 9. Switches :
173  !
174  ! !/S Enable subroutine tracing.
175  !
176  ! 10. Source code :
177  !
178  !/ ------------------------------------------------------------------- /
179  USE w3gdatmd, ONLY: nk,sig,nspec,mapwn
180  USE w3idatmd, ONLY: mudt, mudv, mudd, inflags1
181  USE constants, ONLY: pi,grav,dwat,nu_water
182  USE w3odatmd, ONLY: ndse
183  USE w3servmd, ONLY: extcde
184  USE w3dispmd, ONLY: wavnu1
185 #ifdef W3_S
186  USE w3servmd, ONLY: strace
187 #endif
188  !/
189  IMPLICIT NONE
190  !/
191  !/ ------------------------------------------------------------------- /
192  !/ Parameter list
193  !/
194  REAL, INTENT(IN) :: H_WDEPTH ! water depth
195  REAL, INTENT(IN) :: AC(NSPEC) ! action density
196  INTEGER, INTENT(IN) :: IX, IY
197  REAL, INTENT(OUT) :: S(NSPEC), D(NSPEC)
198  !/
199  !/ ------------------------------------------------------------------- /
200  !/ Local parameters
201  !/
202 #ifdef W3_S
203  INTEGER, SAVE :: IENT = 0
204 #endif
205 
206  COMPLEX :: K
207  COMPLEX :: SHH
208  COMPLEX :: CHH
209  COMPLEX :: SHD
210  COMPLEX :: CHD
211  COMPLEX :: LAM1
212  COMPLEX :: LAM2
213  COMPLEX :: CHLAM2
214  COMPLEX :: SHLAM2
215  COMPLEX :: A1
216  COMPLEX :: A2
217  COMPLEX :: A3
218  COMPLEX :: B1
219  COMPLEX :: B2
220  COMPLEX :: B3
221  COMPLEX :: B4
222  COMPLEX :: C0
223  COMPLEX :: TESTALF1
224  COMPLEX :: TESTALF2
225  COMPLEX :: ALF1
226  COMPLEX :: ALF2
227  COMPLEX :: PSI1
228  COMPLEX :: PSI2
229  COMPLEX :: M1
230  COMPLEX :: M0
231  COMPLEX :: C(4,3)
232  COMPLEX :: C41A
233  COMPLEX :: C42A
234  COMPLEX :: C43A
235  COMPLEX :: B(4)
236  COMPLEX :: CD
237  COMPLEX :: HH
238  COMPLEX :: DD
239  COMPLEX :: GG
240  COMPLEX :: FM1
241  COMPLEX :: KM1
242  COMPLEX :: FP
243  COMPLEX :: I
244  COMPLEX :: F
245  COMPLEX :: KMUD
246 
247  REAL :: BET0
248  REAL :: KINVISW
249  REAL :: RHOW
250  REAL :: EXPH
251  REAL :: A
252  REAL :: K_UNMUD
253  REAL :: SIGMA ! radian frequency (rad)
254  REAL :: SMUDWD(NK) ! dissipation due to mud
255  REAL :: KMIMAG(NK) ! imag part of kmud
256  REAL :: KD
257  REAL :: KDCUTOFF
258  REAL :: CWAVE
259  REAL :: ZTMP
260  REAL :: NWAVE_MUD
261  REAL :: CG_MUD
262  REAL :: KCHECK
263  REAL :: KTHRESHOLD
264  REAL :: RHOM
265  REAL :: KINVISM
266  REAL :: THICKM
267  REAL :: CG_UNMUD
268 
269  INTEGER :: ICOUNT
270  INTEGER :: IK
271  INTEGER :: IS
272 
273  parameter(i=(0.,1.))
274 
275  !/
276  !/ ------------------------------------------------------------------- /
277  !/
278 #ifdef W3_S
279  CALL strace (ient, 'W3SBT8')
280 #endif
281  !
282  ! 0. Initializations ------------------------------------------------ *
283  !
284  ! Dalrymple and Liu, Waves over soft muds: 1978.
285  ! Thin layer solution.
286  ! Matlab code provided by Tony Dalrymple
287  ! Converted to Fortran by Erick Rogers
288 
289  ! Initialize properties from mud fields if available
290  IF (inflags1(-2))THEN
291  rhom = mudd(ix,iy)
292  ELSE
293  WRITE(ndse,*)'RHOM NOT SET'
294  CALL extcde ( 1 )
295  ENDIF
296  IF (inflags1(-1)) THEN
297  thickm = mudt(ix,iy)
298  ELSE
299  WRITE(ndse,*)'THICKM NOT SET'
300  CALL extcde ( 2 )
301  ENDIF
302  IF (inflags1(0)) THEN
303  kinvism = mudv(ix,iy)
304  ELSE
305  WRITE(ndse,*)'KINVISM NOT SET'
306  CALL extcde ( 3 )
307  ENDIF
308 
309  rhow=dwat ! Density of seawater
310  kinvisw=nu_water
311  kdcutoff = 10.0
312  kthreshold=1.0e-9
313 
314  a=1.0
315 
316  ! initialize matrix diagonal contributions
317  d = 0.0
318  s = 0.0
319 
320  IF ( thickm>0.0 .AND. rhom>0.0 .AND. kinvism>0.0 ) THEN
321 
322  smudwd = 0.0
323 
324  ! *** loop over frequencies
325  DO ik = 1,nk
326 
327  sigma = sig(ik)
328 
329  ! un-muddy wave number, to start things off
330  CALL wavnu1(sigma,h_wdepth,k_unmud,cg_unmud)
331  k=k_unmud
332 
333  ! start iterative loop
334 
335  DO icount=1,20 ! *** May need more ***
336 
337  CALL csinh(k*h_wdepth,shh)
338  CALL ccosh(k*h_wdepth,chh)
339  CALL csinh(k*thickm,shd)
340  CALL ccosh(k*thickm,chd)
341 
342  ! define lambdas
343  lam1=sqrt(k*k-i*sigma/kinvisw)
344  lam2=sqrt(k*k-i*sigma/kinvism)
345 
346  ! define hyperbolics on lamda2, lamda1
347  CALL ccosh(lam2*thickm,chlam2)
348  CALL csinh(lam2*thickm,shlam2)
349 
350  ! define exp decay
351  exph=exp(-lam1*h_wdepth)
352 
353  ! define a1, a2, a3
354  a1=-lam2*shd/k+shlam2
355  a2=-chd+chlam2
356  a3=shd-lam2*shlam2/k
357 
358  ! define b1, b2, b3, b4
359  b1=lam1*shh/k-chh
360  b2=lam2*a2*shh/k+a1*chh
361  b3=-a3*shh+a2*chh
362  b4=lam1*shh/k+chh
363 
364  ! define c0
365  c0=b4*exph-(lam1*lam1+k*k)/(2*k*k)
366 
367  ! define beta0
368  bet0=-exph/c0
369 
370  ! define alfa1, alfa2
371  testalf1=-rhom*kinvism*(lam2*lam2+k*k)*(-lam2/k)*chd/k &
372  -2*rhom*kinvism*lam2*chlam2-(rhom-rhow)*grav*(i*(a1)/sigma)
373  testalf2=-rhom*kinvism*(lam2*lam2+k*k)*(-1)*shd/k &
374  -2*rhom*kinvism*lam2*shlam2-(rhom-rhow)*grav*(i*(a2)/sigma)
375  alf1=-testalf1
376  alf2=-testalf2
377 
378  ! define psi1, psi2
379  psi1=2*k*(-lam2/k)*shd+(lam2*lam2+k*k)*shlam2/k
380  psi2=2*k*(-1)*chd+(lam2*lam2+k*k)*chlam2/k
381 
382  ! define M1, MO
383  m1=i*rhow*sigma/k-2*rhow*kinvisw*k
384  m0=b1+(lam1*lam1+k*k)*exph/(k*k)
385 
386  ! matrix coefficients (eq. 22)
387  c(1,1)=lam1*(bet0*m0+1)*shh/k+(bet0*m0-1)*chh+m0/c0+exph
388  c(1,2)=(lam1*bet0*b2+lam2*a2)*shh/k+(bet0*b2+a1)*chh+b2/c0
389  c(1,3)=(lam1*bet0*b3/k-a3)*shh+(bet0*b3+a2)*chh+b3/c0
390 
391  ! matrix coefficients (eq. 23)
392  c(2,1)=lam1*(bet0*m0+1)*m1*chh/k+(bet0*m0-1)*m1*shh &
393  -2*rhow*kinvisw*lam1*m0/c0+2*rhow*kinvisw*lam1*exph
394  c(2,2)=(lam1*bet0*b2+lam2*a2)*m1*chh/k+(bet0*b2+a1)*m1*shh &
395  -2*rhow*kinvisw*lam1*b2/c0
396  c(2,3)=(lam1*bet0*b3/k-a3)*m1*chh+(bet0*b3+a2)*m1*shh &
397  -2*rhow*kinvisw*lam1*b3/c0
398 
399  ! matrix coefficients (eq. 21)
400  c(3,1)=2*k*rhow*kinvisw*(bet0*m0-1)+rhow*kinvisw &
401  *(lam1*lam1+k*k)*(1-bet0*m0)/k
402  c(3,2)=2*k*rhow*kinvisw*(bet0*b2+a1)-rhow*kinvisw &
403  *(lam1*lam1+k*k)*bet0*b2/k-i*rhom*kinvism*psi1
404  c(3,3)=2*k*rhow*kinvisw*(bet0*b3+a2)-rhow*kinvisw &
405  *(lam1*lam1+k*k)*bet0*b3/k-i*rhom*kinvism*psi2
406 
407  ! matrix coefficients (eq.19)
408  c41a=lam1*m1*(bet0*m0+1)/k+2*rhow*kinvisw*lam1 &
409  +2*rhow*kinvisw*lam1*bet0*m0
410  c42a=m1*(lam1*bet0*b2+lam2*a2)/k &
411  +2*rhow*kinvisw*lam1*bet0*b2+alf1
412  c43a=m1*(lam1*bet0*b3/k-a3)+2*rhow*kinvisw*lam1*bet0*b3+alf2
413 
414  ! method 1
415  c(4,1)=c41a*c(3,1)/c41a-c(3,1)
416  c(4,2)=c42a*c(3,1)/c41a-c(3,2)
417  c(4,3)=c43a*c(3,1)/c41a-c(3,3)
418 
419  ! force terms......righthand side
420  b(1)=-i*sigma*a
421  b(2)=rhow*grav*a
422  b(3)=0
423  b(4)=0
424  ! coefficients
425  cd=-(c(3,3)-c(3,2)*c(4,3)/c(4,2))/c(3,1)
426  hh=b(2)/(c(2,1)*cd -c(2,2)*(c(4,3)/c(4,2))+c(2,3))
427  dd=cd*hh
428  gg=-c(4,3)*hh/c(4,2)
429 
430  ! find k
431  f=c(1,1)*dd+c(1,2)*gg+c(1,3)*hh-b(1)
432 
433  IF(icount.EQ.1)THEN
434  fm1=f
435  km1=k
436  k=k*(.995)+.001*i
437  kcheck=100.0
438  ELSE
439  kcheck=abs(imag(k)-imag(km1))
440  IF((f.EQ.fm1).OR.(k.EQ.km1).OR.(kcheck<kthreshold))THEN
441  ! notes: I have noticed that if iterations are not stopped early enough, NaNs result
442  kmud=k
443  EXIT
444  END IF
445  fp=(f-fm1)/(k-km1)
446  km1=k
447  fm1=f
448  k=k-0.8*f/fp
449  ENDIF
450  END DO
451 
452  kmud=k
453 
454  ! KD calc: not that important: just used to determine whether we make the mud calculation
455  kd = real(kmud) * h_wdepth
456 
457  IF ( kd .LT. kdcutoff ) THEN
458 
459  ! Notes: It would be better to have CG_MUD stored for each freq.
460  cwave=sigma/real(kmud)
461  ztmp=2.0*real(kmud)*h_wdepth
462  IF(ztmp.LT.70)THEN
463  ztmp=sinh(ztmp)
464  ELSE
465  ztmp=1.0e+30
466  ENDIF
467  nwave_mud=0.5*(1.0+2.0*real(kmud)*h_wdepth/ztmp)
468  cg_mud=nwave_mud*cwave
469  ! --- compute fluid mud-induced wave dissipation
470  smudwd(ik)=2.0*imag(kmud)*cg_mud
471 
472  ! --- store imaginary part of KMUD
473  kmimag(ik)=imag(kmud)
474 
475  END IF ! IF ( KD .LT. KDCUTOFF ) THEN
476 
477  END DO ! IK
478 
479  ! *** store the results in the DIAGONAL array D ***
480  DO is = 1,nspec
481  d(is) = -smudwd(mapwn(is))
482  END DO
483 
484  END IF ! IF ( THICKM>0.0 & RHOM>0.0 & KINVISM>0.0 ) THEN
485 
486  s = d * ac
487 
488  RETURN
489  !/
490  !/ End of W3SBT8 ----------------------------------------------------- /
491  !/

References ccosh(), csinh(), constants::dwat, w3servmd::extcde(), constants::grav, w3idatmd::inflags1, w3gdatmd::mapwn, w3odatmd::ndse, w3gdatmd::nk, w3gdatmd::nspec, constants::nu_water, constants::pi, w3gdatmd::sig, w3servmd::strace(), and w3dispmd::wavnu1().

Referenced by gxexpo(), and w3srcemd::w3srce().

constants::pi
real, parameter pi
PI Value of Pi.
Definition: constants.F90:71
w3idatmd::inflags1
logical, dimension(:), pointer inflags1
Definition: w3idatmd.F90:260
w3gdatmd::sig
real, dimension(:), pointer sig
Definition: w3gdatmd.F90:1234
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3servmd
Definition: w3servmd.F90:3
constants::nu_water
real, parameter nu_water
NU_WATER Kinematic viscosity of water (m2/s).
Definition: constants.F90:67
w3odatmd
Definition: w3odatmd.F90:3
w3sbt8md::csinh
subroutine csinh(C, CS)
Complex hyperbolic sin (sinh).
Definition: w3sbt8md.F90:505
constants::dwat
real, parameter dwat
DWAT Density of water (kg/m3).
Definition: constants.F90:62
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3idatmd
Define data structures to set up wave model input data for several models simultaneously.
Definition: w3idatmd.F90:16
w3gdatmd::mapwn
integer, dimension(:), pointer mapwn
Definition: w3gdatmd.F90:1231
constants
Define some much-used constants for global use (all defined as PARAMETER).
Definition: constants.F90:20
w3gdatmd
Definition: w3gdatmd.F90:16
w3dispmd::wavnu1
subroutine wavnu1(SI, H, K, CG)
Definition: w3dispmd.F90:85
w3dispmd
Definition: w3dispmd.F90:3
w3sbt8md::ccosh
subroutine ccosh(C, CC)
Complex hyperbolic cos (cosh).
Definition: w3sbt8md.F90:524
constants::grav
real, parameter grav
GRAV Acc.
Definition: constants.F90:61