ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
bbl_mod Module Reference

Functions/Subroutines

subroutine, public bblm (ng, tile)
 
subroutine mb_bbl_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, h, z_r, z_w, angler, zobot, hwave, uwave_rms, dwave, pwave_bot, rho, u, v, bottom, ubot, vbot, ur, vr, bustrc, bvstrc, bustrw, bvstrw, bustrcwmax, bvstrcwmax, bustr, bvstr)
 
subroutine sg_bbl_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, h, z_r, z_w, angler, zobot, hwave, uwave_rms, dwave, pwave_bot, rho, u, v, bottom, iconv, ubot, vbot, ur, vr, bustrc, bvstrc, bustrw, bvstrw, bustrcwmax, bvstrcwmax, bustr, bvstr)
 
subroutine sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, sg_ubouc, sg_mu, sg_epsilon, sg_ro, sg_fofx)
 
subroutine sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
 
subroutine sg_kelvin8m (x, ber, bei, ker, kei, berp, beip, kerp, keip)
 
subroutine sg_kelvin8p (x, ker, kei, ber, bei, kerp, keip, berp, beip)
 
subroutine ssw_bbl_tile (ng, tile, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, nrhs, h, z_r, z_w, angler, zobot, hwave, uwave_rms, dwave, pwave_bot, bedldu, bedldv, bottom, rho, u, v, u_stokes, v_stokes, zeta, iconv, ubot, vbot, ur, vr, bustrc, bvstrc, bustrw, bvstrw, bustrcwmax, bvstrcwmax, bustr, bvstr)
 
subroutine madsen94 (ubr, wr, ucr, zr, phiwc, kn, dstp, ustrc, ustrwm, ustrr, fwc, zoa, dwc_va)
 
subroutine log_interp (kmax, dstp, sg_loc, u_1d, v_1d, z_r_1d, z_w_1d, d50, zapp_loc, zr_sg, ur_sg, vr_sg)
 

Function/Subroutine Documentation

◆ bblm()

subroutine public bbl_mod::bblm ( integer, intent(in) ng,
integer, intent(in) tile )

Definition at line 54 of file mb_bbl.h.

55!***********************************************************************
56!
57 USE mod_param
58 USE mod_bbl
59 USE mod_forces
60 USE mod_grid
61 USE mod_ocean
62 USE mod_sedbed
63 USE mod_stepping
64!
65! Imported variable declarations.
66!
67 integer, intent(in) :: ng, tile
68!
69! Local variable declarations.
70!
71 character (len=*), parameter :: MyFile = &
72 & __FILE__
73!
74# include "tile.h"
75!
76# ifdef PROFILE
77 CALL wclock_on (ng, inlm, 37, __line__, myfile)
78# endif
79 CALL mb_bbl_tile (ng, tile, &
80 & lbi, ubi, lbj, ubj, &
81 & imins, imaxs, jmins, jmaxs, &
82 & nrhs(ng), &
83 & grid(ng) % h, &
84 & grid(ng) % z_r, &
85 & grid(ng) % z_w, &
86 & grid(ng) % angler, &
87 & grid(ng) % ZoBot, &
88# ifdef MB_CALC_UB
89 & forces(ng) % Hwave, &
90# else
91 & forces(ng) % Uwave_rms, &
92# endif
93 & forces(ng) % Dwave, &
94 & forces(ng) % Pwave_bot, &
95 & ocean(ng) % rho, &
96 & ocean(ng) % u, &
97 & ocean(ng) % v, &
98 & sedbed(ng) % bottom, &
99 & bbl(ng) % Ubot, &
100 & bbl(ng) % Vbot, &
101 & bbl(ng) % Ur, &
102 & bbl(ng) % Vr, &
103 & bbl(ng) % bustrc, &
104 & bbl(ng) % bvstrc, &
105 & bbl(ng) % bustrw, &
106 & bbl(ng) % bvstrw, &
107 & bbl(ng) % bustrcwmax, &
108 & bbl(ng) % bvstrcwmax, &
109 & forces(ng) % bustr, &
110 & forces(ng) % bvstr)
111# ifdef PROFILE
112 CALL wclock_off (ng, inlm, 37, __line__, myfile)
113# endif
114!
115 RETURN
mod_bbl
Definition mod_bbl.F:2
mod_bbl::bbl
type(t_bbl), dimension(:), allocatable bbl
Definition mod_bbl.F:62
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_sedbed
Definition sedbed_mod.h:1
mod_sedbed::sedbed
type(t_sedbed), dimension(:), allocatable sedbed
Definition sedbed_mod.h:157
mod_stepping
Definition mod_stepping.F:2
mod_stepping::nrhs
integer, dimension(:), allocatable nrhs
Definition mod_stepping.F:70
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3

References mod_bbl::bbl, mod_forces::forces, mod_grid::grid, mod_param::inlm, mb_bbl_tile(), mod_stepping::nrhs, mod_ocean::ocean, mod_sedbed::sedbed, wclock_off(), and wclock_on().

Referenced by main3d().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ log_interp()

subroutine bbl_mod::log_interp ( integer, intent(in) kmax,
real(r8), intent(in) dstp,
real(r8), intent(in) sg_loc,
real(r8), dimension(1:kmax), intent(in) u_1d,
real(r8), dimension(1:kmax), intent(in) v_1d,
real(r8), dimension(1:kmax), intent(in) z_r_1d,
real(r8), dimension(0:kmax), intent(in) z_w_1d,
real(r8), intent(in) d50,
real(r8), intent(in) zapp_loc,
real(r8), intent(out) zr_sg,
real(r8), intent(out) ur_sg,
real(r8), intent(out) vr_sg )
private

Definition at line 1740 of file ssw_bbl.h.

1745!
1746!=======================================================================
1747! !
1748! Find the near-bottom current velocity in x- and y- directions !
1749! that corresponds to user input elevation to get near-bottom !
1750! current velocity (m) !
1751! !
1752!=======================================================================
1753!
1754 USE mod_param
1755 USE mod_sediment
1756 USE mod_scalars
1757!
1758 implicit none
1759!
1760! Imported variable declarations.
1761!
1762 integer, intent(in) :: kmax
1763 real(r8), intent(in) :: Dstp, sg_loc
1764 real(r8), intent(in) :: u_1d(1:kmax), v_1d(1:kmax)
1765 real(r8), intent(in) :: z_r_1d(1:kmax), z_w_1d(0:kmax)
1766 real(r8), intent(in) :: d50, zapp_loc
1767 real(r8), intent(out) :: Zr_sg
1768 real(r8), intent(out) :: Ur_sg, Vr_sg
1769!
1770! Local variables.
1771!
1772 integer :: k
1773 real(r8) :: z1, z2, Zr
1774 real(r8) :: fac, fac1, fac2
1775!
1776 zr=z_r_1d(1)-z_w_1d(0)
1777!
1778 IF ( sg_loc.ge.zr ) THEN
1779!
1780! If chosen height to get near bottom-current velocity lies
1781! within any vertical level, perform logarithmic interpolation.
1782!
1783 DO k=2,kmax
1784 z1=z_r_1d(k-1)-z_w_1d(0)
1785 z2=z_r_1d(k )-z_w_1d(0)
1786 IF ( ( z1.le.sg_loc ).and.( sg_loc.lt.z2 )) THEN
1787 fac=1.0_r8/log(z2/z1)
1788 fac1=fac*log(z2/sg_loc)
1789 fac2=fac*log(sg_loc/z1)
1790!
1791 ur_sg=fac1*u_1d(k-1)+fac2*u_1d(k)
1792 vr_sg=fac1*v_1d(k-1)+fac2*v_1d(k)
1793 zr_sg=sg_loc
1794 ENDIF
1795 IF ((k.eq.kmax).and.(sg_loc.ge.z2)) THEN
1796 ur_sg=u_1d(k)
1797 vr_sg=v_1d(k)
1798 zr_sg=z2
1799 END IF
1800 END DO
1801 ELSE !IF ( sg_loc.lt.Zr ) THEN
1802 z1=max( 2.5_r8*d50/30.0_r8, zapp_loc, 1.0e-10_r8 )
1803 z2=zr
1804!
1805 IF ( sg_loc.lt.z1 ) THEN
1806!
1807! If chosen height is less than the bottom roughness
1808! perform linear interpolation.
1809!
1810 z1=sg_loc
1811 fac=z1/z2
1812!
1813 ur_sg=fac*u_1d(1)
1814 vr_sg=fac*v_1d(1)
1815 zr_sg=sg_loc
1816!
1817 ELSE !IF ( sg_loc.gt.z1 ) THEN
1818!
1819! If chosen height is less than the bottom cell thickness
1820! perform logarithmic interpolation with bottom roughness.
1821!
1822 fac=1.0_r8/log(z2/z1)
1823 fac2=fac*log(sg_loc/z1)
1824!
1825 ur_sg=fac2*u_1d(1)
1826 vr_sg=fac2*v_1d(1)
1827 zr_sg=sg_loc
1828 END IF
1829 END IF
1830!
1831 RETURN
mod_scalars
Definition mod_scalars.F:2
mod_sediment
Definition sediment_mod.h:1

Referenced by ssw_bbl_tile().

Here is the caller graph for this function:

◆ madsen94()

subroutine bbl_mod::madsen94 ( real(r8), intent(in) ubr,
real(r8), intent(in) wr,
real(r8), intent(in) ucr,
real(r8), intent(in) zr,
real(r8), intent(in) phiwc,
real(r8), intent(inout) kn,
real(r8), intent(in) dstp,
real(r8), intent(out) ustrc,
real(r8), intent(out) ustrwm,
real(r8), intent(out) ustrr,
real(r8), intent(out) fwc,
real(r8), intent(out) zoa,
real(r8), intent(out) dwc_va )
private

Definition at line 1532 of file ssw_bbl.h.

1534!
1535!=======================================================================
1536! !
1537! Grant-Madsen model from Madsen (1994). !
1538! !
1539! On Input: !
1540! !
1541! ubr Rep. wave-orbital velocity amplitude outside WBL (m/s). !
1542! wr Rep. angular wave frequency, 2* pi/T (rad/s). !
1543! ucr Current velocity at height zr (m/s). !
1544! zr Reference height for current velocity (m). !
1545! phiwc Angle between currents and waves at zr (radians). !
1546! kN Bottom roughness height, like Nikuradse k, (m). !
1547! Dstp Total water depth (m). !
1548! !
1549! On Output: !
1550! !
1551! ustrc Current friction velocity, u*c (m/s). !
1552! ustrwm Wave maximum friction velocity, u*wm (m/s). !
1553! ustrr Wave-current combined friction velocity, u*r (m/s). !
1554! fwc Wave friction factor (nondimensional). !
1555! zoa Apparent bottom roughness (m). !
1556! dwc_va Wave-boundary layer thickness (m). !
1557! !
1558!=======================================================================
1559!
1560 USE mod_param
1561 USE mod_scalars
1562!
1563! Imported variable declarations.
1564!
1565 real(r8), intent(in) :: ubr, wr, ucr, zr, phiwc, Dstp
1566 real(r8), intent(inout) :: kN
1567 real(r8), intent(out) :: ustrc, ustrwm, ustrr, fwc, zoa
1568 real(r8), intent(out) :: dwc_va
1569!
1570! Local variable declarations.
1571!
1572 integer, parameter :: MAXIT = 20
1573 integer :: i, nit
1574 integer :: iverbose = 1
1575
1576 real(r8) :: bigsqr, cosphiwc, cukw, diff, lndw, lnln, lnzr
1577 real(r8) :: phicwc, zo
1578 real(r8) :: dval = 99.99_r8
1579
1580 real(r8), dimension(MAXIT) :: Cmu
1581 real(r8), dimension(MAXIT) :: dwc
1582 real(r8), dimension(MAXIT) :: fwci
1583 real(r8), dimension(MAXIT) :: rmu
1584 real(r8), dimension(MAXIT) :: ustrci
1585 real(r8), dimension(MAXIT) :: ustrr2
1586 real(r8), dimension(MAXIT) :: ustrwm2
1587!
1588!-----------------------------------------------------------------------
1589! Compute bottom friction velocities and roughness.
1590!-----------------------------------------------------------------------
1591!
1592! Set special default values.
1593!
1594 ustrc=dval
1595 ustrwm=dval
1596 ustrr=dval
1597# ifdef CRS_FIX
1598 kn=min(kn,0.9_r8*zr)
1599# endif
1600 zoa=kn/30.0_r8
1601 phicwc=phiwc
1602
1603 zo = kn/30.0_r8
1604
1605 IF (ubr.le.0.01_r8) THEN
1606 dwc_va=kn
1607 zoa=dwc_va
1608 fwc=0.0_r8
1609 IF (ucr.le. 0.01_r8) THEN ! no waves or currents
1610 ustrc=0.0_r8
1611 ustrwm=0.0_r8
1612 ustrr=0.0_r8
1613 RETURN
1614 END IF
1615 ustrc=ucr*vonkar/log(zr/zo) ! no waves
1616 ustrwm=0.0_r8
1617 ustrr=ustrc
1618 RETURN
1619 END IF
1620!
1621! Iterate to compute friction velocities, roughness, and wave friction
1622! factor. Notice that the computation of the wave friction factor
1623! has been inlined for efficiency.
1624!
1625 cosphiwc=abs(cos(phiwc))
1626 rmu(1)=0.0_r8
1627 cmu(1)=1.0_r8
1628
1629 cukw=cmu(1)*ubr/(kn*wr)
1630# if defined CRS_FIX
1631!
1632! New fwc CRS calculation
1633!
1634 fwci(1)=cmu(1)*0.3_r8
1635 IF ((cukw.gt.0.352_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1636 fwci(1)=cmu(1)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1637 ELSE IF (cukw.gt.100.0_r8) THEN
1638 fwci(1)=cmu(1)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1639 END IF
1640# else
1641!
1642! Original method fwc calculation
1643!
1644 IF ((cukw.gt.0.2_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1645 fwci(1)=cmu(1)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1646 ELSE IF ((cukw.gt.100.).and.(cukw.le.10000.0_r8)) THEN
1647 fwci(1)=cmu(1)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1648 ELSE IF (cukw.gt.10000.0_r8 ) THEN
1649 fwci(1)=cmu(1)*exp(5.61_r8*10000.0_r8**(-0.109_r8)-7.30_r8)
1650 ELSE
1651 fwci(1)=cmu(1)*0.43_r8
1652 END IF
1653# endif
1654!
1655 ustrwm2(1)=0.5_r8*fwci(1)*ubr*ubr ! Eq 29
1656 ustrr2(1)=cmu(1)*ustrwm2(1) ! Eq 26
1657 ustrr=sqrt(ustrr2(1))
1658 IF (cukw.ge.8.0_r8) THEN
1659 dwc(1)=2.0_r8*vonkar*ustrr/wr
1660# if defined CRS_FIX
1661 dwc(1)=min( 0.9_r8*zr, dwc(1) )
1662# endif
1663 ELSE
1664 dwc(1)=kn
1665 END IF
1666 lnzr=log(zr/dwc(1))
1667 lndw=log(dwc(1)/zo)
1668 lnln=lnzr/lndw
1669 bigsqr=-1.0_r8+sqrt(1.0_r8+((4.0_r8*vonkar*lndw)/ &
1670 & (lnzr*lnzr))*ucr/ustrr)
1671 ustrci(1)=0.5_r8*ustrr*lnln*bigsqr
1672!
1673 i=1
1674 diff=1.0_r8
1675 DO WHILE ((i.lt.maxit).and.(diff.gt.0.000005_r8))
1676 i=i+1
1677 rmu(i)=ustrci(i-1)*ustrci(i-1)/ustrwm2(i-1)
1678 cmu(i)=sqrt(1.0_r8+ &
1679 & 2.0_r8*rmu(i)*cosphiwc+rmu(i)*rmu(i)) ! Eq 27
1680 cukw=cmu(i)*ubr/(kn*wr)
1681# ifdef CRS_FIX
1682!
1683! New fwc CRS calculation
1684!
1685 fwci(i)=cmu(i)*0.3_r8
1686 IF ((cukw.gt.0.352_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1687 fwci(i)=cmu(i)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1688 ELSE IF (cukw.gt.100.0_r8) THEN
1689 fwci(i)=cmu(i)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1690 END IF
1691# else
1692!
1693! Original method fwc calculation
1694!
1695 IF ((cukw.gt.0.2_r8).and.(cukw.le.100.0_r8)) THEN ! Eq 32/33
1696 fwci(i)=cmu(i)*exp(7.02_r8*cukw**(-0.078_r8)-8.82_r8)
1697 ELSE IF ((cukw.gt.100.).and.(cukw.le.10000.0_r8)) THEN
1698 fwci(i)=cmu(i)*exp(5.61_r8*cukw**(-0.109_r8)-7.30_r8)
1699 ELSE IF (cukw.gt.10000.0_r8 ) THEN
1700 fwci(i)=cmu(i)*exp(5.61_r8*10000.0_r8**(-0.109_r8)-7.30_r8)
1701 ELSE
1702 fwci(i)=cmu(i)*0.43_r8
1703 END IF
1704# endif
1705!
1706 ustrwm2(i)=0.5_r8*fwci(i)*ubr*ubr ! Eq 29
1707 ustrr2(i)=cmu(i)*ustrwm2(i) ! Eq 26
1708 ustrr=sqrt(ustrr2(i))
1709!! IF ((Cmu(1)*ubr/(kN*wr)).ge.8.0_r8) THEN ! HGA Why 1?
1710 IF (cukw.ge.8.0_r8) THEN
1711 dwc(i)=2.0_r8*vonkar*ustrr/wr ! Eq 36
1712# if defined CRS_FIX
1713 dwc(i)=min( 0.9_r8*zr, dwc(i) )
1714# endif
1715 ELSE
1716 dwc(i)=kn
1717 END IF
1718 lnzr=log(zr/dwc(i))
1719 lndw=log(dwc(i)/zo)
1720 lnln=lnzr/lndw
1721 bigsqr=-1.0_r8+sqrt(1.0_r8+((4.0_r8*vonkar*lndw)/ &
1722 & (lnzr*lnzr))*ucr/ustrr)
1723 ustrci(i)=0.5_r8*ustrr*lnln*bigsqr ! Eq 38
1724 diff=abs((fwci(i)-fwci(i-1))/fwci(i))
1725 END DO
1726 ustrwm=sqrt(ustrwm2(i))
1727 ustrc=ustrci(i)
1728 ustrr=sqrt(ustrr2(i))
1729 phicwc=phiwc
1730 zoa=exp(log(dwc(i))-(ustrc/ustrr)*log(dwc(i)/zo)) ! Eq 11
1731 dwc_va=dwc(i)
1732 fwc=fwci(i)
1733!
1734 RETURN
mod_scalars::vonkar
real(dp) vonkar
Definition mod_scalars.F:469

References mod_scalars::vonkar.

Referenced by ssw_bbl_tile().

Here is the caller graph for this function:

◆ mb_bbl_tile()

subroutine bbl_mod::mb_bbl_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) h,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) angler,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zobot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) hwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) uwave_rms,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) dwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pwave_bot,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) v,
real(r8), dimension(lbi:ubi,lbj:ubj,mbotp), intent(inout) bottom,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ubot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vbot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ur,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstr )
private

Definition at line 119 of file mb_bbl.h.

136!***********************************************************************
137!
138 USE mod_param
139 USE mod_scalars
140 USE mod_sediment
141!
142 USE bc_2d_mod
143# ifdef DISTRIBUTE
144 USE mp_exchange_mod, ONLY : mp_exchange2d
145# endif
146!
147! Imported variable declarations.
148!
149 integer, intent(in) :: ng, tile
150 integer, intent(in) :: LBi, UBi, LBj, UBj
151 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
152 integer, intent(in) :: nrhs
153!
154# ifdef ASSUMED_SHAPE
155 real(r8), intent(in) :: h(LBi:,LBj:)
156 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
157 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
158 real(r8), intent(in) :: angler(LBi:,LBj:)
159 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
160# ifdef MB_CALC_UB
161 real(r8), intent(in) :: Hwave(LBi:,LBj:)
162# else
163 real(r8), intent(in) :: Uwave_rms(LBi:,LBj:)
164# endif
165 real(r8), intent(in) :: Dwave(LBi:,LBj:)
166 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
167 real(r8), intent(in) :: rho(LBi:,LBj:,:)
168 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
169 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
170
171 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
172
173 real(r8), intent(out) :: Ubot(LBi:,LBj:)
174 real(r8), intent(out) :: Vbot(LBi:,LBj:)
175 real(r8), intent(out) :: Ur(LBi:,LBj:)
176 real(r8), intent(out) :: Vr(LBi:,LBj:)
177 real(r8), intent(out) :: bustrc(LBi:,LBj:)
178 real(r8), intent(out) :: bvstrc(LBi:,LBj:)
179 real(r8), intent(out) :: bustrw(LBi:,LBj:)
180 real(r8), intent(out) :: bvstrw(LBi:,LBj:)
181 real(r8), intent(out) :: bustrcwmax(LBi:,LBj:)
182 real(r8), intent(out) :: bvstrcwmax(LBi:,LBj:)
183 real(r8), intent(out) :: bustr(LBi:,LBj:)
184 real(r8), intent(out) :: bvstr(LBi:,LBj:)
185# else
186 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
187 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
188 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
189 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
190 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
191# ifdef MB_CALC_UB
192 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
193# else
194 real(r8), intent(in) :: Uwave_rms(LBi:UBi,LBj:UBj)
195# endif
196 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
197 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
198 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
199 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
200 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
201
202 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
203
204 real(r8), intent(out) :: Ubot(LBi:UBi,LBj:UBj)
205 real(r8), intent(out) :: Vbot(LBi:UBi,LBj:UBj)
206 real(r8), intent(out) :: Ur(LBi:UBi,LBj:UBj)
207 real(r8), intent(out) :: Vr(LBi:UBi,LBj:UBj)
208 real(r8), intent(out) :: bustrc(LBi:UBi,LBj:UBj)
209 real(r8), intent(out) :: bvstrc(LBi:UBi,LBj:UBj)
210 real(r8), intent(out) :: bustrw(LBi:UBi,LBj:UBj)
211 real(r8), intent(out) :: bvstrw(LBi:UBi,LBj:UBj)
212 real(r8), intent(out) :: bustrcwmax(LBi:UBi,LBj:UBj)
213 real(r8), intent(out) :: bvstrcwmax(LBi:UBi,LBj:UBj)
214 real(r8), intent(out) :: bustr(LBi:UBi,LBj:UBj)
215 real(r8), intent(out) :: bvstr(LBi:UBi,LBj:UBj)
216# endif
217!
218! Local variable declarations.
219!
220 integer :: i, ised, j
221
222 real(r8), parameter :: K1 = 0.6666666666_r8
223 real(r8), parameter :: K2 = 0.3555555555_r8
224 real(r8), parameter :: K3 = 0.1608465608_r8
225 real(r8), parameter :: K4 = 0.0632098765_r8
226 real(r8), parameter :: K5 = 0.0217540484_r8
227 real(r8), parameter :: K6 = 0.0065407983_r8
228
229 real(r8), parameter :: eps = 1.0e-10_r8
230
231 real(r8), parameter :: scf1 = 0.5_r8 * 1.39_r8
232 real(r8), parameter :: scf2 = 0.52_r8
233 real(r8), parameter :: scf3 = 2.0_r8 - scf2
234 real(r8), parameter :: scf4 = 1.2_r8
235 real(r8), parameter :: scf5 = 3.2_r8
236
237 real(r8) :: twopi = 2.0_r8 * pi
238
239 real(r8) :: RHbiomax = 0.006_r8 ! maximum biogenic ripple height
240 real(r8) :: RHmin = 0.001_r8 ! arbitrary minimum ripple height
241 real(r8) :: RLmin = 0.01_r8 ! arbitrary minimum ripple length
242
243 real(r8) :: Ab, Fwave_bot, Kbh, Kbh2, Kdh
244 real(r8) :: angleC, angleW, phiC, phiCW
245 real(r8) :: cff, cff1, cff2, d50, viscosity, wset
246 real(r8) :: RHbio, RHbiofac, RHmax, RLbio
247 real(r8) :: rhoWater, rhoSed
248 real(r8) :: rhgt, rlen, thetw
249 real(r8) :: Znot, ZnotC, Znot_bl, Znot_rip
250 real(r8) :: tau_bf, tau_ex, tau_en, tau_up, tau_w, tau_wb
251 real(r8) :: tau_c, tau_cb, tau_cs, tau_cw, tau_cwb, tau_cws
252
253 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ub
254 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ucur
255 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vcur
256 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr
257 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_mb
258 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_mb
259 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Umag
260 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tauC
261 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tauW
262 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tauCW
263 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: tauCWmax
264
265#include "set_bounds.h"
266!
267!-----------------------------------------------------------------------
268! Initalize stresses due to currents and waves.
269!-----------------------------------------------------------------------
270!
271 rhbiofac=1.0_r8/exp(4.11_r8)
272 DO j=jstrv-1,jend
273 DO i=istru-1,iend
274 tauc(i,j)=0.0_r8
275 taucw(i,j)=0.0_r8
276 taucwmax(i,j)=0.0_r8
277 END DO
278 END DO
279!
280!-----------------------------------------------------------------------
281! Set currents above bed.
282!-----------------------------------------------------------------------
283!
284 DO j=jstrv-1,jend+1
285 DO i=istru-1,iend+1
286 zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
287 ur_mb(i,j)=u(i,j,1,nrhs)
288 vr_mb(i,j)=v(i,j,1,nrhs)
289 END DO
290 END DO
291!
292!=======================================================================
293! Compute bottom stresses.
294!=======================================================================
295!
296 DO j=jstrv-1,jend
297 DO i=istru-1,iend
298 rhbio=0.0_r8
299 rlbio=0.0_r8
300 rlen=bottom(i,j,irlen)
301 rhgt=bottom(i,j,irhgt)
302 rhowater=rho(i,j,1)+1000.0_r8
303 viscosity=0.0013_r8/rhowater ! kinematic viscosity
304!
305!-----------------------------------------------------------------------
306! Compute bed wave orbital velocity (m/s) and excursion amplitude (m)
307! from wind-induced waves. Use Dean and Dalrymple (1991) 6th-degree
308! polynomial to approximate wave number on shoaling water.
309!-----------------------------------------------------------------------
310!
311 fwave_bot=twopi/max(pwave_bot(i,j),0.05_r8)
312# ifdef MB_BBL_CALC_UB
313 kdh=h(i,j)*fwave_bot*fwave_bot/g
314 kbh2=kdh*kdh+ &
315 & kdh/(1.0_r8+kdh*(k1+kdh*(k2+kdh*(k3+kdh*(k4+ &
316 & kdh*(k5+k6*kdh))))))
317 kbh=sqrt(kbh2)
318!
319! Compute bed wave orbital velocity and excursion amplitude.
320!
321 ab=0.5_r8*hwave(i,j)/sinh(kbh)+eps
322 ub(i,j)=fwave_bot*ab
323# else
324 ub(i,j)=uwave_rms(i,j)
325 ab=ub(i,j)/fwave_bot+eps
326# endif
327!
328! Compute bottom current magnitude at RHO-points.
329!
330 ucur(i,j)=0.5_r8*(ur_mb(i,j)+ur_mb(i+1,j))
331 vcur(i,j)=0.5_r8*(vr_mb(i,j)+vr_mb(i,j+1))
332 umag(i,j)=sqrt(ucur(i,j)*ucur(i,j)+vcur(i,j)*vcur(i,j))+eps
333!
334! Compute angle between currents and waves (radians)
335!
336 IF (ucur(i,j).ne.0.0_r8) THEN
337 phic=atan2(vcur(i,j),ucur(i,j))
338 ELSE
339 phic=0.5_r8*pi*sign(1.0_r8,vcur(i,j))
340 END IF
341 phicw=1.5_r8*pi-dwave(i,j)-phic-angler(i,j)
342!
343!-----------------------------------------------------------------------
344! Determine skin roughness from sediment size. Set default logarithmic
345! profile consistent with current-only case.
346! Establish local median grain size for all calculations here and
347! determine local values of critical stresses.
348! Since most parameterizations have been derived ignoring multiple
349! grain sizes, we apply this single d50 also in the case of mixed beds.
350!-----------------------------------------------------------------------
351!
352 d50=bottom(i,j,isd50) ! (m)
353 tau_cb=bottom(i,j,itauc) ! (m2/s2)
354 wset=bottom(i,j,iwsed) ! (m/s)
355 rhosed=bottom(i,j,idens)/rhowater ! (nondimensional)
356!
357! Critical stress (m2/s2) for transition to sheet flow
358! (Li and Amos, 2001, eq. 11).
359!
360 tau_up=0.172_r8*(rhosed-1.0_r8)*g*(d50**0.624_r8)
361!
362! Threshold stress (m2/s2) for break off (Grant and Madsen,1982).
363!
364 tau_bf=0.79_r8*(viscosity**(-0.6_r8))* &
365 & (((rhosed-1.0_r8)*g)**0.3_r8)*(d50**0.9_r8)*tau_cb
366!
367! Set Znot for currents as maximum of user value or grain roughness.
368!
369 znotc=d50/12.0_r8
370 znot=max(zobot(i,j),znotc)
371!
372!-----------------------------------------------------------------------
373! Set tauC (m2/s2) acc. to current-only case (skin friction) [m/s]
374!-----------------------------------------------------------------------
375!
376 cff1=vonkar/log(zr(i,j)/znot)
377 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
378 tauc(i,j)=cff2*umag(i,j)*umag(i,j)
379
380 cff1=vonkar/log(zr(i,j)/znotc)
381 tau_cs=cff1*cff1*umag(i,j)*umag(i,j)
382!
383!-----------------------------------------------------------------------
384! If significant waves (Ub > 0.01 m/s), wave-current interaction case
385! according to Soulsby (1995). Otherwise, tauCWmax = tauC for sediment
386! purposes.
387!-----------------------------------------------------------------------
388!
389 IF (ub(i,j).gt.0.01_r8) THEN
390!
391! Determine skin stresses (m2/s2) for pure waves and combined flow
392! using Soulsby (1995) approximation of the wave friction factor,
393! fw=2*scf1*(Znot/Ab)**scf2 so tauW=0.5fw*Ub**2.
394!
395 tau_w=scf1*((znotc*fwave_bot)**scf2)*(ub(i,j)**scf3)
396!
397! Wave-averaged, combined wave-current stress.(Eqn. 69, Soulsby, 1997).
398!
399 tau_cw=tau_cs* &
400 & (1.0_r8+scf4*((tau_w/(tau_w+tau_cs))**scf5))
401!
402! Maximum of combined wave-current skin stress (m2/s2) component for
403! sediment.(Eqn. 70, Soulsby, 1997).
404!
405 tau_cws=sqrt((tau_cw+tau_w*cos(phicw))**2+ &
406 & (tau_w*sin(phicw))**2)
407!! tauw(i,j)=tau_cws
408 taucwmax(i,j)=tau_cws
409 tauw(i,j)=tau_w
410!
411! Set combined stress for Znot.
412!
413 tau_w=scf1*((znot*fwave_bot)**scf2)*(ub(i,j)**scf3)
414 tau_cw=tauc(i,j)* &
415 & (1.0_r8+scf4*((tau_w/(tau_w+tauc(i,j)))**scf5))
416
417# ifdef MB_Z0BL
418!
419!-----------------------------------------------------------------------
420! Compute bedload roughness for ripple predictor and sediment purposes.
421! At high transport stages, friction depends on thickness of bedload
422! layer. Shear stress due to combined grain and bedload roughness is
423! used to predict ripples and onset of suspension (Li and Amos, 2001).
424!-----------------------------------------------------------------------
425!
426# ifdef MB_CALC_ZNOT
427 tau_ex=max((tau_cws-tau_cb),0.0_r8)
428 cff=(1.0_r8/((rhosed-1.0_r8)*g*d50))
429 znot_bl=17.4_r8*d50*(cff*tau_ex)**0.75_r8
430 znotc=znotc+znot_bl
431# endif
432!
433!-----------------------------------------------------------------------
434! Compute stresses (m2/s2)for sediment purposes, using grain and
435! bedload roughness.
436!-----------------------------------------------------------------------
437!
438 cff1=vonkar/log(zr(i,j)/znotc)
439 tau_c=cff1*cff1*umag(i,j)*umag(i,j)
440 tau_wb=scf1*((znotc*fwave_bot)**scf2)*(ub(i,j)**scf3)
441 tau_cw=tau_c*(1.0_r8+scf4*((tau_wb/(tau_wb+tau_c))**scf5))
442!
443! Maximum of combined wave-current stress (m2/s2) component for
444! sediment purposes.
445!
446 tau_cwb=sqrt((tau_cw+tau_wb*cos(phicw))**2+ &
447 & (tau_wb*sin(phicw))**2)
448 taucwmax(i,j)=tau_cwb
449 tauw(i,j)=tau_wb
450# endif
451
452# ifdef MB_Z0RIP
453!
454!-----------------------------------------------------------------------
455! Determine bedform roughness ripple height (m) and ripple length (m)
456! for sandy beds. Use structure according to Li and Amos (2001).
457!-----------------------------------------------------------------------
458!
459! Check median grain diameter
460!
461 IF (d50.ge.0.000063_r8) THEN
462!
463! Enhanced skin stress if pre-exisiting ripples (Nielsen, 1986).
464!
465 rhmax=0.8_r8*rlen/pi
466 rhgt=max(min(rhmax,rhgt),rhmin)
467 tau_en=max(tau_cws,tau_cws*(rlen/(rlen-pi*rhgt))**2)
468!
469 IF ((tau_cws.lt.tau_cb).and.(tau_en.ge.tau_cb)) THEN
470 rhgt=(19.6_r8*(sqrt(tau_cws/tau_cb))+20.9_r8)*d50
471 rlen=rhgt/0.12_r8 ! local transport
472 ELSE
473 IF ((tau_cws.ge.tau_cb).and.(tau_cwb.lt.tau_bf)) THEN
474 rhgt=(22.15_r8*(sqrt(tau_cwb/tau_cb))+6.38_r8)*d50
475 rlen=rhgt/0.12_r8 ! bed load in eq. range
476 ELSE
477 IF ((tau_cwb.ge.tau_bf).and.(tau_cwb.lt.tau_up)) THEN
478 rlen=535.0_r8*d50 ! break off regime
479 rhgt=0.15_r8*rlen*(sqrt(tau_up)-sqrt(tau_cwb))/ &
480 & (sqrt(tau_up)-sqrt(tau_bf ))
481 ELSE IF (tau_cwb.ge.tau_up) THEN
482 rlen=0.0_r8 ! sheet flow, plane bed
483 rhgt=0.0_r8
484 ELSE
485 rlen=bottom(i,j,irlen)
486 rhgt=bottom(i,j,irhgt)
487 END IF
488 END IF
489 END IF
490 END IF
491# endif
492
493# ifdef MB_Z0BIO
494!
495!-----------------------------------------------------------------------
496! Determine (biogenic) bedform roughness ripple height (m) and ripple
497! length (m) for silty beds, using Harris and Wiberg (2001).
498!-----------------------------------------------------------------------
499!
500! Use 10 cm default biogenic ripple length, RLbio (Wheatcroft 1994).
501!
502! NOTE: For mixed beds take average of silty and sandy bedform
503! roughnesses weighted according to silt and sand fractions.
504!
505 IF (d50.lt.0.000063_r8) THEN
506 rlbio=0.1_r8
507 thetw=tau_cws*(1.0_r8/((rhosed-1.0_r8)*g*d50))
508 rhbio=(thetw**(-1.67_r8))*rlbio*rhbiofac
509 rhgt=min(rhbio,rhbiomax)
510 rlen=rlbio
511 END IF
512# endif
513
514# if defined MB_Z0RIP || defined MB_Z0BIO
515!
516! Ripple roughness using Grant and Madsen (1982) roughness length.
517!
518# ifdef MB_CALC_ZNOT
519 znot_rip=0.92_r8*rhgt*rhgt/(max(rlen,rlmin))
520 znotc=znotc+znot_rip
521# endif
522!
523!-----------------------------------------------------------------------
524! Compute bottom stress (m2/s2) components based on total roughnes.
525!-----------------------------------------------------------------------
526!
527 cff1=vonkar/log(zr(i,j)/znotc)
528 tau_c=cff1*cff1*umag(i,j)*umag(i,j)
529 tau_w=scf1*((znotc*fwave_bot)**scf2)*(ub(i,j)**scf3)
530 tau_cw=tau_c* &
531 & (1.0_r8+scf4*((tau_w/(tau_w+tau_c))**scf5))
532!! tau_cwb=SQRT((tau_cw+tau_w*COS(phiCW))**2+ &
533!! & (tau_w*SIN(phiCW))**2)
534!! tauCWmax(i,j)=tau_cwb
535# endif
536!
537!-----------------------------------------------------------------------
538! Compute effective bottom shear velocity (m/s) relevant for flow and
539! eddy-diffusivities/viscosity.
540!-----------------------------------------------------------------------
541!
542!! tauC(i,j)=tau_cw
543 taucw(i,j)=tau_cw
544 tauw(i,j)=tau_w
545 ELSE IF (ub(i,j).le.0.01_r8) THEN
546!
547! If current-only, tauCWmax=tauC(skin) for use in sediment model; tauC
548! is determined using roughness due to current ripples.
549! In this limit, tauCWmax=tauCW=tauC
550!
551 taucwmax(i,j)=tauc(i,j)
552 tauw(i,j)=0.0_r8
553!! tauW(i,j)=tau_cs
554!! tauCW(i,j)=tau_cs
555
556# ifdef MB_Z0RIP
557!! IF (tauC(i,j).gt.tau_up) THEN
558 IF (tau_cs(i,j).gt.tau_up) THEN
559 rhgt=0.0_r8
560 rlen=0.0_r8
561!! ELSE IF (tauC(i,j).lt.tau_cb) THEN
562 ELSE IF (tau_cs(i,j).lt.tau_cb) THEN
563 rlen=bottom(i,j,irlen)
564 rhgt=bottom(i,j,irhgt)
565 ELSE
566 rlen=1000.0_r8*d50 ! Yalin (1964)
567 rhgt=0.0308_r8*(rlen**1.19_r8)
568!! rhgt=100.0_r8*0.074_r8* &
569!! & (0.01_r8*rlen)**1.19_r8 ! Allen (1970)
570 END IF
571# ifdef MB_CALC_ZNOT
572 znotc=znotc+0.92_r8*rhgt*rhgt/(max(rlen,rlmin))
573# endif
574# endif
575 cff1=vonkar/log(zr(i,j)/znotc)
576 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
577!! tauc(i,j)=cff2*Umag(i,j)*Umag(i,j)
578 taucw(i,j)=cff2*umag(i,j)*umag(i,j)
579 END IF
580!
581!-----------------------------------------------------------------------
582! Load variables for output purposes.
583!-----------------------------------------------------------------------
584!
585 bottom(i,j,ibwav)=ab
586 bottom(i,j,irlen)=rlen
587 bottom(i,j,irhgt)=rhgt
588 bottom(i,j,izdef)=znot ! Zob(ng)
589 bottom(i,j,izapp)=znotc
590 END DO
591 END DO
592!
593!-----------------------------------------------------------------------
594! Compute kinematic bottom stress (m2/s2) components for flow due to
595! combined current and wind-induced waves.
596!-----------------------------------------------------------------------
597!
598 DO j=jstr,jend
599 DO i=istru,iend
600 anglec=ur_mb(i,j)/(0.5*(umag(i-1,j)+umag(i,j)))
601 bustr(i,j)=0.5_r8*(taucw(i-1,j)+taucw(i,j))*anglec
602# ifdef WET_DRY
603 cff2=0.75_r8*0.5_r8*(z_w(i-1,j,1)+z_w(i,j,1)- &
604 & z_w(i-1,j,0)-z_w(i,j,0))
605 bustr(i,j)=sign(1.0_r8,bustr(i,j))*min(abs(bustr(i,j)), &
606 & abs(u(i,j,1,nrhs))*cff2/dt(ng))
607# endif
608 END DO
609 END DO
610 DO j=jstrv,jend
611 DO i=istr,iend
612 anglec=vr_mb(i,j)/(0.5_r8*(umag(i,j-1)+umag(i,j)))
613 bvstr(i,j)=0.5_r8*(taucw(i,j-1)+taucw(i,j))*anglec
614# ifdef WET_DRY
615 cff2=0.75_r8*0.5_r8*(z_w(i,j-1,1)+z_w(i,j,1)- &
616 & z_w(i,j-1,0)-z_w(i,j,0))
617 bvstr(i,j)=sign(1.0_r8,bvstr(i,j))*min(abs(bvstr(i,j)), &
618 & abs(v(i,j,1,nrhs))*cff2/dt(ng))
619# endif
620 END DO
621 END DO
622 DO j=jstr,jend
623 DO i=istr,iend
624 anglec=ucur(i,j)/umag(i,j)
625 anglew=cos(1.5_r8*pi-dwave(i,j)-angler(i,j))
626 bustrc(i,j)=taucw(i,j)*anglec
627 bustrw(i,j)=tauw(i,j)*anglew
628 bustrcwmax(i,j)=taucwmax(i,j)*anglew
629 ubot(i,j)=ub(i,j)*anglew
630 ur(i,j)=ucur(i,j)
631!
632 anglec=vcur(i,j)/umag(i,j)
633 anglew=sin(1.5_r8*pi-dwave(i,j)-angler(i,j))
634 bvstrc(i,j)=taucw(i,j)*anglec
635 bvstrw(i,j)=tauw(i,j)*anglew
636 bvstrcwmax(i,j)=taucwmax(i,j)*anglew
637 vbot(i,j)=ub(i,j)*anglew
638 vr(i,j)=vcur(i,j)
639 END DO
640 END DO
641!
642! Apply periodic or gradient boundary conditions for output purposes.
643!
644 CALL bc_u2d_tile (ng, tile, &
645 & lbi, ubi, lbj, ubj, &
646 & bustr)
647 CALL bc_v2d_tile (ng, tile, &
648 & lbi, ubi, lbj, ubj, &
649 & bvstr)
650 CALL bc_r2d_tile (ng, tile, &
651 & lbi, ubi, lbj, ubj, &
652 & bustrc)
653 CALL bc_r2d_tile (ng, tile, &
654 & lbi, ubi, lbj, ubj, &
655 & bvstrc)
656 CALL bc_r2d_tile (ng, tile, &
657 & lbi, ubi, lbj, ubj, &
658 & bustrw)
659 CALL bc_r2d_tile (ng, tile, &
660 & lbi, ubi, lbj, ubj, &
661 & bvstrw)
662 CALL bc_u2d_tile (ng, tile, &
663 & lbi, ubi, lbj, ubj, &
664 & bustrcwmax)
665 CALL bc_r2d_tile (ng, tile, &
666 & lbi, ubi, lbj, ubj, &
667 & bvstrcwmax)
668 CALL bc_r2d_tile (ng, tile, &
669 & lbi, ubi, lbj, ubj, &
670 & ubot)
671 CALL bc_r2d_tile (ng, tile, &
672 & lbi, ubi, lbj, ubj, &
673 & vbot)
674 CALL bc_r2d_tile (ng, tile, &
675 & lbi, ubi, lbj, ubj, &
676 & ur)
677 CALL bc_r2d_tile (ng, tile, &
678 & lbi, ubi, lbj, ubj, &
679 & vr)
680 CALL bc_r2d_tile (ng, tile, &
681 & lbi, ubi, lbj, ubj, &
682 & bottom(:,:,ibwav))
683 CALL bc_r2d_tile (ng, tile, &
684 & lbi, ubi, lbj, ubj, &
685 & bottom(:,:,irhgt))
686 CALL bc_r2d_tile (ng, tile, &
687 & lbi, ubi, lbj, ubj, &
688 & bottom(:,:,irlen))
689 CALL bc_r2d_tile (ng, tile, &
690 & lbi, ubi, lbj, ubj, &
691 & bottom(:,:,izdef))
692 CALL bc_r2d_tile (ng, tile, &
693 & lbi, ubi, lbj, ubj, &
694 & bottom(:,:,izapp))
695
696#ifdef DISTRIBUTE
697 CALL mp_exchange2d (ng, tile, inlm, 4, &
698 & lbi, ubi, lbj, ubj, &
699 & nghostpoints, &
700 & ewperiodic(ng), nsperiodic(ng), &
701 & bustr, bvstr, bustrc, bvstrc)
702 CALL mp_exchange2d (ng, tile, inlm, 4, &
703 & lbi, ubi, lbj, ubj, &
704 & nghostpoints, &
705 & ewperiodic(ng), nsperiodic(ng), &
706 & bustrw, bvstrw, bustrcwmax, bvstrcwmax)
707 CALL mp_exchange2d (ng, tile, inlm, 4, &
708 & lbi, ubi, lbj, ubj, &
709 & nghostpoints, &
710 & ewperiodic(ng), nsperiodic(ng), &
711 & ubot, vbot, ur, vr)
712 CALL mp_exchange2d (ng, tile, inlm, 3, &
713 & lbi, ubi, lbj, ubj, &
714 & nghostpoints, &
715 & ewperiodic(ng), nsperiodic(ng), &
716 & bottom(:,:,ibwav), &
717 & bottom(:,:,irhgt), &
718 & bottom(:,:,irlen))
719 CALL mp_exchange2d (ng, tile, inlm, 2, &
720 & lbi, ubi, lbj, ubj, &
721 & nghostpoints, &
722 & ewperiodic(ng), nsperiodic(ng), &
723 & bottom(:,:,izdef), &
724 & bottom(:,:,izapp))
725#endif
726!
727 RETURN
bc_2d_mod
Definition bc_2d.F:2
bc_2d_mod::bc_r2d_tile
subroutine bc_r2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:44
bc_2d_mod::bc_v2d_tile
subroutine bc_v2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:345
bc_2d_mod::bc_u2d_tile
subroutine bc_u2d_tile(ng, tile, lbi, ubi, lbj, ubj, a)
Definition bc_2d.F:167
mod_param::nghostpoints
integer nghostpoints
Definition mod_param.F:710
mod_scalars::cdb_min
real(dp) cdb_min
Definition mod_scalars.F:772
mod_scalars::dt
real(dp), dimension(:), allocatable dt
Definition mod_scalars.F:269
mod_scalars::ewperiodic
logical, dimension(:), allocatable ewperiodic
Definition mod_scalars.F:510
mod_scalars::nsperiodic
logical, dimension(:), allocatable nsperiodic
Definition mod_scalars.F:511
mod_scalars::g
real(dp) g
Definition mod_scalars.F:464
mod_scalars::cdb_max
real(dp) cdb_max
Definition mod_scalars.F:773
mod_scalars::pi
real(dp), parameter pi
Definition mod_scalars.F:813
mod_sediment::iwsed
integer, parameter iwsed
Definition sediment_mod.h:126
mod_sediment::irlen
integer, parameter irlen
Definition sediment_mod.h:128
mod_sediment::irhgt
integer, parameter irhgt
Definition sediment_mod.h:129
mod_sediment::ibwav
integer, parameter ibwav
Definition sediment_mod.h:130
mod_sediment::idens
integer, parameter idens
Definition sediment_mod.h:125
mod_sediment::isd50
integer, parameter isd50
Definition sediment_mod.h:124
mod_sediment::izapp
integer, parameter izapp
Definition sediment_mod.h:132
mod_sediment::izdef
integer, parameter izdef
Definition sediment_mod.h:131
mod_sediment::itauc
integer, parameter itauc
Definition sediment_mod.h:127
mp_exchange_mod
Definition mp_exchange.F:2
mp_exchange_mod::mp_exchange2d
subroutine mp_exchange2d(ng, tile, model, nvar, lbi, ubi, lbj, ubj, nghost, ew_periodic, ns_periodic, a, b, c, d)
Definition mp_exchange.F:294

References bc_2d_mod::bc_r2d_tile(), bc_2d_mod::bc_u2d_tile(), bc_2d_mod::bc_v2d_tile(), mod_scalars::cdb_max, mod_scalars::cdb_min, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::g, mod_sediment::ibwav, mod_sediment::idens, mod_param::inlm, mod_sediment::irhgt, mod_sediment::irlen, mod_sediment::isd50, mod_sediment::itauc, mod_sediment::iwsed, mod_sediment::izapp, mod_sediment::izdef, mp_exchange_mod::mp_exchange2d(), mod_param::nghostpoints, mod_scalars::nsperiodic, mod_scalars::pi, and mod_scalars::vonkar.

Referenced by bblm().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sg_bbl_tile()

subroutine bbl_mod::sg_bbl_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) h,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) angler,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zobot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) hwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) uwave_rms,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) dwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pwave_bot,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) v,
real(r8), dimension(lbi:ubi,lbj:ubj,mbotp), intent(inout) bottom,
integer, dimension(lbi:ubi,lbj:ubj), intent(inout) iconv,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ubot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vbot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ur,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstr )
private

Definition at line 97 of file sg_bbl.h.

116!***********************************************************************
117!
118 USE mod_param
119 USE mod_scalars
120 USE mod_sediment
121!
122 USE bc_2d_mod
123#ifdef DISTRIBUTE
124 USE mp_exchange_mod, ONLY : mp_exchange2d
125#endif
126!
127! Imported variable declarations.
128!
129 integer, intent(in) :: ng, tile
130 integer, intent(in) :: LBi, UBi, LBj, UBj
131 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
132 integer, intent(in) :: nrhs
133!
134#ifdef ASSUMED_SHAPE
135 integer, intent(inout) :: Iconv(LBi:,LBj:)
136
137 real(r8), intent(in) :: h(LBi:,LBj:)
138 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
139 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
140 real(r8), intent(in) :: angler(LBi:,LBj:)
141 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
142# if defined SG_CALC_UB
143 real(r8), intent(in) :: Hwave(LBi:,LBj:)
144# else
145 real(r8), intent(in) :: Uwave_rms(LBi:,LBj:)
146# endif
147 real(r8), intent(in) :: Dwave(LBi:,LBj:)
148 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
149 real(r8), intent(in) :: rho(LBi:,LBj:,:)
150 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
151 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
152
153 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
154
155 real(r8), intent(out) :: Ubot(LBi:,LBj:)
156 real(r8), intent(out) :: Vbot(LBi:,LBj:)
157 real(r8), intent(out) :: Ur(LBi:,LBj:)
158 real(r8), intent(out) :: Vr(LBi:,LBj:)
159 real(r8), intent(out) :: bustrc(LBi:,LBj:)
160 real(r8), intent(out) :: bvstrc(LBi:,LBj:)
161 real(r8), intent(out) :: bustrw(LBi:,LBj:)
162 real(r8), intent(out) :: bvstrw(LBi:,LBj:)
163 real(r8), intent(out) :: bustrcwmax(LBi:,LBj:)
164 real(r8), intent(out) :: bvstrcwmax(LBi:,LBj:)
165 real(r8), intent(out) :: bustr(LBi:,LBj:)
166 real(r8), intent(out) :: bvstr(LBi:,LBj:)
167#else
168 integer, intent(inout) :: Iconv(LBi:UBi,LBj:UBj)
169
170 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
171 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
172 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
173 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
174 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
175# if defined SG_CALC_UB
176 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
177# else
178 real(r8), intent(in) :: Uwave_rms(LBi:UBi,LBj:UBj)
179# endif
180 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
181 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
182 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
183 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
184 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
185
186 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
187
188 real(r8), intent(out) :: Ubot(LBi:UBi,LBj:UBj)
189 real(r8), intent(out) :: Vbot(LBi:UBi,LBj:UBj)
190 real(r8), intent(out) :: Ur(LBi:UBi,LBj:UBj)
191 real(r8), intent(out) :: Vr(LBi:UBi,LBj:UBj)
192 real(r8), intent(out) :: bustrc(LBi:UBi,LBj:UBj)
193 real(r8), intent(out) :: bvstrc(LBi:UBi,LBj:UBj)
194 real(r8), intent(out) :: bustrw(LBi:UBi,LBj:UBj)
195 real(r8), intent(out) :: bvstrw(LBi:UBi,LBj:UBj)
196 real(r8), intent(out) :: bustrcwmax(LBi:UBi,LBj:UBj)
197 real(r8), intent(out) :: bvstrcwmax(LBi:UBi,LBj:UBj)
198 real(r8), intent(out) :: bustr(LBi:UBi,LBj:UBj)
199 real(r8), intent(out) :: bvstr(LBi:UBi,LBj:UBj)
200#endif
201!
202! Local variable declarations.
203!
204 logical :: ITERATE
205
206 integer :: Iter, i, j, k
207
208 real(r8), parameter :: eps = 1.0e-10_r8
209
210 real(r8) :: Fwave_bot, Kb, Kbh, KboKb0, Kb0, Kdelta, Ustr
211 real(r8) :: anglec, anglew
212 real(r8) :: cff, cff1, cff2, cff3, og, fac, fac1, fac2
213 real(r8) :: sg_ab, sg_abokb, sg_a1, sg_b1, sg_chi, sg_c1, sg_dd
214 real(r8) :: sg_epsilon, sg_eta, sg_fofa, sg_fofb, sg_fofc, sg_fwm
215 real(r8) :: sg_kbs, sg_lambda, sg_mu, sg_phicw, sg_ro, sg_row
216 real(r8) :: sg_shdnrm, sg_shld, sg_shldcr, sg_scf, sg_ss, sg_star
217 real(r8) :: sg_ub, sg_ubokur, sg_ubouc, sg_ubouwm, sg_ur
218 real(r8) :: sg_ustarc, sg_ustarcw, sg_ustarwm, sg_znot, sg_znotp
219 real(r8) :: sg_zr, sg_zrozn, sg_z1, sg_z1ozn, sg_z2, twopi, z1, z2
220
221 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ab
222 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauc
223 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauw
224 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taucwmax
225 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sg
226 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sg
227 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ub
228 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ucur
229 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Umag
230 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vcur
231 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr
232 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic
233 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phicw
234 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rheight
235 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rlength
236 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u100
237 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znot
238 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znotc
239
240#include "set_bounds.h"
241!
242!-----------------------------------------------------------------------
243! Initalize to default values.
244!-----------------------------------------------------------------------
245!
246 DO j=jstrv-1,jend
247 DO i=istru-1,iend
248 tauc(i,j)=0.0_r8
249 tauw(i,j)=0.0_r8
250 taucwmax(i,j)=0.0_r8
251 u100(i,j)=0.0_r8
252 rheight(i,j)=0.0_r8
253 rlength(i,j)=0.0_r8
254 znot(i,j)=zobot(i,j)
255 znotc(i,j)=0.0_r8
256 END DO
257 END DO
258!
259!-----------------------------------------------------------------------
260! Set currents above bed.
261!-----------------------------------------------------------------------
262!
263 DO j=jstrv-1,jend+1
264 DO i=istru-1,iend+1
265 zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
266 ur_sg(i,j)=u(i,j,1,nrhs)
267 vr_sg(i,j)=v(i,j,1,nrhs)
268#ifdef SG_LOGINT
269!
270! If current height is less than z1ur, interpolate logarithmically
271! to z1ur.
272!
273 IF (zr(i,j).lt.sg_z1min) THEN
274 DO k=2,n(ng)
275 z1=z_r(i,j,k-1)-z_w(i,j,0)
276 z2=z_r(i,j,k )-z_w(i,j,0)
277 IF ((z1.lt.sg_z1min).and.(sg_z1min.lt.z2)) THEN
278 fac=1.0_r8/log(z2/z1)
279 fac1=fac*log(z2/sg_z1min)
280 fac2=fac*log(sg_z1min/z1)
281 ur_sg(i,j)=fac1*u(i,j,k-1,nrhs)+fac2*u(i,j,k,nrhs)
282 vr_sg(i,j)=fac1*v(i,j,k-1,nrhs)+fac2*v(i,j,k,nrhs)
283 zr(i,j)=sg_z1min
284 END IF
285 END DO
286 END IF
287#endif
288 END DO
289 END DO
290!
291!-----------------------------------------------------------------------
292! Compute bed wave orbital velocity (m/s) and excursion amplitude
293! (m) from wind-induced waves. Use linear wave theory dispersion
294! relation for wave number.
295!-----------------------------------------------------------------------
296!
297 twopi=2.0_r8*pi
298 og=1.0_r8/g
299 DO j=jstrv-1,jend
300 DO i=istru-1,iend
301!
302! Compute first guess for wavenumber, Kb0. Use deep water (Kb0*h>1)
303! and shallow water (Kb0*H<1) approximations.
304!
305 fwave_bot=twopi/max(pwave_bot(i,j),0.05_r8)
306!
307! Compute bed wave orbital velocity and excursion amplitude.
308!
309#ifdef SG_CALC_UB
310 kb0=fwave_bot*fwave_bot*og
311 IF (kb0*h(i,j).ge.1.0_r8) THEN
312 kb=kb0
313 ELSE
314 kb=fwave_bot/sqrt(g*h(i,j))
315 END IF
316!
317! Compute bottom wave number via Newton-Raphson method.
318!
319 iterate=.true.
320 DO iter=1,sg_n
321 IF (iterate) THEN
322 kbh=kb*h(i,j)
323 kbokb0=kb/kb0
324 kdelta=(1.0_r8-kbokb0*tanh(kbh))/ &
325 & (1.0_r8+kbh*(kbokb0-1.0_r8/kbokb0))
326 iterate=abs(kb*kdelta) .ge. sg_tol
327 kb=kb*(1.0_r8+kdelta)
328 END IF
329 END DO
330 ab(i,j)=0.5_r8*hwave(i,j)/sinh(kb*h(i,j))+eps
331 ub(i,j)=fwave_bot*ab(i,j)+eps
332#else
333 ub(i,j)=abs(uwave_rms(i,j))+eps
334 ab(i,j)=ub(i,j)/fwave_bot+eps
335#endif
336!
337! Compute bottom current magnitude at RHO-points.
338!
339 ucur(i,j)=0.5_r8*(ur_sg(i,j)+ur_sg(i+1,j))
340 vcur(i,j)=0.5_r8*(vr_sg(i,j)+vr_sg(i,j+1))
341 umag(i,j)=sqrt(ucur(i,j)*ucur(i,j)+vcur(i,j)*vcur(i,j))+eps
342!
343! Compute angle between currents and waves (radians)
344!
345 IF (ucur(i,j).eq.0.0_r8) THEN
346 phic(i,j)=0.5_r8*pi*sign(1.0_r8,vcur(i,j))
347 ELSE
348 phic(i,j)=atan2(vcur(i,j),ucur(i,j))
349 ENDIF
350 phicw(i,j)=1.5_r8*pi-dwave(i,j)-phic(i,j)-angler(i,j)
351 END DO
352 END DO
353!
354!-----------------------------------------------------------------------
355! Set default logarithmic profile.
356!-----------------------------------------------------------------------
357!
358 DO j=jstrv-1,jend
359 DO i=istru-1,iend
360 IF (umag(i,j).gt.0.0_r8) THEN
361!! Ustr=MIN(sg_ustarcdef,Umag(i,j)*vonKar/ &
362!! & LOG(Zr(i,j)/ZoBot(i,j)))
363!! Tauc(i,j)=Ustr*Ustr
364 cff1=vonkar/log(zr(i,j)/zobot(i,j))
365 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
366 tauc(i,j)=cff2*umag(i,j)*umag(i,j)
367 END IF
368 END DO
369 END DO
370!
371!-----------------------------------------------------------------------
372! Wave-current interaction case.
373!-----------------------------------------------------------------------
374!
375 DO j=jstrv-1,jend
376 DO i=istru-1,iend
377 sg_dd=bottom(i,j,isd50)
378 sg_ss=bottom(i,j,idens)/(rho(i,j,1)+1000.0_r8)
379 sg_ab=ab(i,j)
380 sg_ub=ub(i,j)
381 sg_phicw=phicw(i,j)
382 sg_ur=umag(i,j)
383 sg_zr=zr(i,j)
384!
385! Compute hydraulic roughness "Znot" (m), ripple height "eta" (m),
386! and ripple length "lambda" (m).
387!
388#ifdef SG_CALC_ZNOT
389 sg_star=sg_dd/(4.0_r8*sg_nu)*sqrt((sg_ss-1.0_r8)*sg_g*sg_dd)
390!
391! Compute critical shield parameter based on grain diameter.
392! (sg_scf is a correction factor).
393!
394 sg_scf=1.0_r8
395 IF (sg_star.le.1.5_r8) THEN
396 sg_shldcr=sg_scf*0.0932_r8*sg_star**(-0.707_r8)
397 ELSE IF ((1.5_r8.lt.sg_star).and.(sg_star.lt.4.0_r8)) THEN
398 sg_shldcr=sg_scf*0.0848_r8*sg_star**(-0.473_r8)
399 ELSE IF ((4.0_r8.le.sg_star).and.(sg_star.lt.10.0_r8)) THEN
400 sg_shldcr=sg_scf*0.0680_r8*sg_star**(-0.314_r8)
401 ELSE IF ((10.0_r8.le.sg_star).and.(sg_star.lt.34.0_r8)) THEN
402 sg_shldcr=sg_scf*0.033_r8
403 ELSE IF ((34.0_r8.le.sg_star).and.(sg_star.lt.270.0_r8)) THEN
404 sg_shldcr=sg_scf*0.0134_r8*sg_star**(0.255_r8)
405 ELSE
406 sg_shldcr=sg_scf*0.056_r8
407 END IF
408!
409! Calculate skin friction shear stress based on Ole Madsen (1994)
410! empirical formula. Check initiation of sediment motion criteria,
411! to see if we compute sg_znot based on the wave-formed ripples.
412! If the skin friction calculation indicates that sediment is NOT
413! in motion, the ripple model is invalid and take the default value,
414! ZoBot.
415!
416 sg_abokb=sg_ab/sg_dd
417 IF (sg_abokb.le.100.0_r8) THEN
418 sg_fwm=exp(7.02_r8*sg_abokb**(-0.078_r8)-8.82_r8)
419 ELSE
420 sg_fwm=exp(5.61_r8*sg_abokb**(-0.109_r8)-7.30_r8)
421 END IF
422 sg_ustarwm=sqrt(0.5_r8*sg_fwm)*sg_ub
423 sg_shdnrm=(sg_ss-1.0_r8)*sg_dd*sg_g
424 sg_shld=sg_ustarwm*sg_ustarwm/sg_shdnrm
425 IF ((sg_shld/sg_shldcr).le.1.0_r8) THEN
426 sg_znot=zobot(i,j)
427 sg_eta=0.0_r8
428 sg_lambda=0.0_r8
429 ELSE
430!
431! Calculate ripple height and length and bottom roughness
432!
433 sg_chi=4.0_r8*sg_nu*sg_ub*sg_ub/ &
434 & (sg_dd*((sg_ss-1.0_r8)*sg_g*sg_dd)**1.5_r8)
435 IF (sg_chi.le.2.0_r8) THEN
436 sg_eta=sg_ab*0.30_r8*sg_chi**(-0.39_r8)
437 sg_lambda=sg_ab*1.96_r8*sg_chi**(-0.28_r8)
438 ELSE
439 sg_eta=sg_ab*0.45_r8*sg_chi**(-0.99_r8)
440 sg_lambda=sg_ab*2.71_r8*sg_chi**(-0.75_r8)
441 END IF
442 sg_kbs=sg_ab*0.0655_r8* &
443 & (sg_ub*sg_ub/((sg_ss-1.0_r8)*sg_g*sg_ab))**1.4_r8
444 sg_znot=(sg_dd+2.3_r8*sg_eta+sg_kbs)/30.0_r8
445 END IF
446#else
447 sg_znot=zobot(i,j)
448 sg_chi=4.0_r8*sg_nu*sg_ub*sg_ub/ &
449 & (sg_dd*((sg_ss-1.0_r8)*sg_g*sg_dd)**1.5_r8)
450 IF (sg_chi.le.2.0_r8) THEN
451 sg_eta=sg_ab*0.32_r8*sg_chi**(-0.34_r8)
452 sg_lambda=sg_ab*2.04_r8*sg_chi**(-0.23_r8)
453 ELSE
454 sg_eta=sg_ab*0.52_r8*sg_chi**(-1.01_r8)
455 sg_lambda=sg_ab*2.7_r8*sg_chi**(-0.78_r8)
456 END IF
457#endif
458 znot(i,j)=sg_znot
459 rheight(i,j)=sg_eta
460 rlength(i,j)=sg_lambda
461!
462! Compute only when nonzero currents and waves.
463!
464 sg_zrozn=sg_zr/sg_znot
465 IF ((sg_ur.gt.0.0_r8).and.(sg_ub.gt.0.0_r8).and. &
466 & (sg_zrozn.gt.1.0_r8)) THEN
467!
468! Compute bottom stress based on ripple roughness.
469!
470 sg_ubokur=sg_ub/(sg_kappa*sg_ur)
471 sg_row=sg_ab/sg_znot
472 sg_a1=1.0e-6_r8
473 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
474 & sg_a1, sg_mu, sg_epsilon, sg_ro, sg_fofa)
475 sg_abokb=sg_ab/(30.0_r8*sg_znot)
476 IF (sg_abokb.le.100.0_r8) THEN
477 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
478 ELSE
479 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
480 END IF
481 sg_ubouwm=sqrt(2.0_r8/sg_fwm)
482!
483! Determine the maximum ratio of wave over combined shear stresses,
484! sg_ubouwm (ub/ustarwm).
485!
486 CALL sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
487!
488! Set initial guess of the ratio of wave over shear stress, sg_c1
489! (ub/ustarc).
490!
491 sg_b1=sg_ubouwm
492 sg_fofb=-sg_fofa
493 sg_c1=0.5_r8*(sg_a1+sg_b1)
494 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
495 & sg_c1, sg_mu, sg_epsilon, sg_ro, sg_fofc)
496!
497! Solve PDE via bi-section method.
498!
499 iterate=.true.
500 DO iter=1,sg_n
501 IF (iterate) THEN
502 IF ((sg_fofb*sg_fofc).lt.0.0_r8) THEN
503 sg_a1=sg_c1
504 ELSE
505 sg_b1=sg_c1
506 END IF
507 sg_c1=0.5_r8*(sg_a1+sg_b1)
508 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
509 & sg_c1, sg_mu, sg_epsilon, sg_ro, &
510 & sg_fofc)
511 iterate=(sg_b1-sg_c1) .ge. sg_tol
512 IF (iterate) iconv(i,j)=iter
513 END IF
514 END DO
515 sg_ubouc=sg_c1
516!
517! Compute bottom shear stress magnitude (m/s).
518!
519 sg_ustarcw=sg_ub/sg_ubouc
520 sg_ustarwm=sg_mu*sg_ustarcw
521!! sg_ustarc=MIN(sg_ustarcdef,sg_epsilon*sg_ustarcw)
522!! sg_ustarc=MIN(SQRT(Tauc(i,j)),sg_epsilon*sg_ustarcw)
523 sg_ustarc=max(sqrt(tauc(i,j)),sg_epsilon*sg_ustarcw)
524 tauc(i,j)=sg_ustarc*sg_ustarc
525 tauw(i,j)=sg_ustarwm*sg_ustarwm
526 taucwmax(i,j)=sqrt((tauc(i,j)+ &
527 & tauw(i,j)*cos(phicw(i,j)))**2+ &
528 & (tauw(i,j)*sin(phicw(i,j)))**2)
529!
530! Compute apparent hydraulic roughness (m).
531!
532 IF (sg_epsilon.gt.0.0_r8) THEN
533 sg_z1=sg_alpha*sg_kappa*sg_ab/sg_ubouc
534 sg_z2=sg_z1/sg_epsilon
535 sg_z1ozn=sg_z1/sg_znot
536 znotc(i,j)=sg_z2* &
537 & exp(-(1.0_r8-sg_epsilon+ &
538 & sg_epsilon*log(sg_z1ozn)))
539!
540! Compute mean (m/s) current at 100 cm above the bottom.
541!
542 IF (sg_z100.gt.sg_z2) THEN
543 u100(i,j)=sg_ustarc* &
544 & (log(sg_z100/sg_z2)+1.0_r8-sg_epsilon+ &
545 & sg_epsilon*log(sg_z1ozn))/sg_kappa
546 ELSE IF ((sg_z100.le.sg_z2).and.(sg_zr.gt.sg_z1)) THEN
547 u100(i,j)=sg_ustarc*sg_epsilon* &
548 & (sg_z100/sg_z1-1.0_r8+log(sg_z1ozn))/sg_kappa
549 ELSE
550 u100(i,j)=sg_ustarc*sg_epsilon* &
551 & log(sg_z100/sg_znot)/sg_kappa
552 END IF
553 END IF
554 END IF
555 END DO
556 END DO
557!
558!-----------------------------------------------------------------------
559! Compute kinematic bottom stress components due current and wind-
560! induced waves.
561!-----------------------------------------------------------------------
562!
563 DO j=jstr,jend
564 DO i=istru,iend
565 anglec=ur_sg(i,j)/(0.5*(umag(i-1,j)+umag(i,j)))
566 bustr(i,j)=0.5_r8*(tauc(i-1,j)+tauc(i,j))*anglec
567# ifdef WET_DRY
568 cff2=0.75_r8*0.5_r8*(z_w(i-1,j,1)+z_w(i,j,1)- &
569 & z_w(i-1,j,0)-z_w(i,j,0))
570 bustr(i,j)=sign(1.0_r8,bustr(i,j))*min(abs(bustr(i,j)), &
571 & abs(u(i,j,1,nrhs))*cff2/dt(ng))
572# endif
573 END DO
574 END DO
575 DO j=jstrv,jend
576 DO i=istr,iend
577 anglec=vr_sg(i,j)/(0.5_r8*(umag(i,j-1)+umag(i,j)))
578 bvstr(i,j)=0.5_r8*(tauc(i,j-1)+tauc(i,j))*anglec
579# ifdef WET_DRY
580 cff2=0.75_r8*0.5_r8*(z_w(i,j-1,1)+z_w(i,j,1)- &
581 & z_w(i,j-1,0)-z_w(i,j,0))
582 bvstr(i,j)=sign(1.0_r8,bvstr(i,j))*min(abs(bvstr(i,j)), &
583 & abs(v(i,j,1,nrhs))*cff2/dt(ng))
584# endif
585 END DO
586 END DO
587 DO j=jstr,jend
588 DO i=istr,iend
589 anglec=ucur(i,j)/umag(i,j)
590 anglew=cos(1.5_r8*pi-dwave(i,j)-angler(i,j))
591 bustrc(i,j)=tauc(i,j)*anglec
592 bustrw(i,j)=tauw(i,j)*anglew
593 bustrcwmax(i,j)=taucwmax(i,j)*anglew
594 ubot(i,j)=ub(i,j)*anglew
595 ur(i,j)=ucur(i,j)
596!
597 anglec=vcur(i,j)/umag(i,j)
598 anglew=sin(1.5_r8*pi-dwave(i,j)-angler(i,j))
599 bvstrc(i,j)=tauc(i,j)*anglec
600 bvstrw(i,j)=tauw(i,j)*anglew
601 bvstrcwmax(i,j)=taucwmax(i,j)*anglew
602 vbot(i,j)=ub(i,j)*anglew
603 vr(i,j)=vcur(i,j)
604!
605 bottom(i,j,ibwav)=ab(i,j)
606 bottom(i,j,irhgt)=rheight(i,j)
607 bottom(i,j,irlen)=rlength(i,j)
608 bottom(i,j,izdef)=znot(i,j)
609 bottom(i,j,izapp)=znotc(i,j)
610 END DO
611 END DO
612!
613! Apply periodic or gradient boundary conditions for output
614! purposes only.
615!
616 CALL bc_u2d_tile (ng, tile, &
617 & lbi, ubi, lbj, ubj, &
618 & bustr)
619 CALL bc_v2d_tile (ng, tile, &
620 & lbi, ubi, lbj, ubj, &
621 & bvstr)
622 CALL bc_r2d_tile (ng, tile, &
623 & lbi, ubi, lbj, ubj, &
624 & bustrc)
625 CALL bc_r2d_tile (ng, tile, &
626 & lbi, ubi, lbj, ubj, &
627 & bvstrc)
628 CALL bc_r2d_tile (ng, tile, &
629 & lbi, ubi, lbj, ubj, &
630 & bustrw)
631 CALL bc_r2d_tile (ng, tile, &
632 & lbi, ubi, lbj, ubj, &
633 & bvstrw)
634 CALL bc_u2d_tile (ng, tile, &
635 & lbi, ubi, lbj, ubj, &
636 & bustrcwmax)
637 CALL bc_r2d_tile (ng, tile, &
638 & lbi, ubi, lbj, ubj, &
639 & bvstrcwmax)
640 CALL bc_r2d_tile (ng, tile, &
641 & lbi, ubi, lbj, ubj, &
642 & ubot)
643 CALL bc_r2d_tile (ng, tile, &
644 & lbi, ubi, lbj, ubj, &
645 & vbot)
646 CALL bc_r2d_tile (ng, tile, &
647 & lbi, ubi, lbj, ubj, &
648 & ur)
649 CALL bc_r2d_tile (ng, tile, &
650 & lbi, ubi, lbj, ubj, &
651 & vr)
652 CALL bc_r2d_tile (ng, tile, &
653 & lbi, ubi, lbj, ubj, &
654 & bottom(:,:,ibwav))
655 CALL bc_r2d_tile (ng, tile, &
656 & lbi, ubi, lbj, ubj, &
657 & bottom(:,:,irhgt))
658 CALL bc_r2d_tile (ng, tile, &
659 & lbi, ubi, lbj, ubj, &
660 & bottom(:,:,irlen))
661 CALL bc_r2d_tile (ng, tile, &
662 & lbi, ubi, lbj, ubj, &
663 & bottom(:,:,izdef))
664 CALL bc_r2d_tile (ng, tile, &
665 & lbi, ubi, lbj, ubj, &
666 & bottom(:,:,izapp))
667#ifdef DISTRIBUTE
668 CALL mp_exchange2d (ng, tile, inlm, 4, &
669 & lbi, ubi, lbj, ubj, &
670 & nghostpoints, &
671 & ewperiodic(ng), nsperiodic(ng), &
672 & bustr, bvstr, bustrc, bvstrc)
673 CALL mp_exchange2d (ng, tile, inlm, 4, &
674 & lbi, ubi, lbj, ubj, &
675 & nghostpoints, &
676 & ewperiodic(ng), nsperiodic(ng), &
677 & bustrw, bvstrw, bustrcwmax, bvstrcwmax)
678 CALL mp_exchange2d (ng, tile, inlm, 4, &
679 & lbi, ubi, lbj, ubj, &
680 & nghostpoints, &
681 & ewperiodic(ng), nsperiodic(ng), &
682 & ubot, vbot, ur, vr)
683 CALL mp_exchange2d (ng, tile, inlm, 3, &
684 & lbi, ubi, lbj, ubj, &
685 & nghostpoints, &
686 & ewperiodic(ng), nsperiodic(ng), &
687 & bottom(:,:,ibwav), &
688 & bottom(:,:,irhgt), &
689 & bottom(:,:,irlen))
690 CALL mp_exchange2d (ng, tile, inlm, 4, &
691 & lbi, ubi, lbj, ubj, &
692 & nghostpoints, &
693 & ewperiodic(ng), nsperiodic(ng), &
694 & bottom(:,:,izdef), &
695 & bottom(:,:,izapp))
696#endif
697!
698 RETURN
mod_param::n
integer, dimension(:), allocatable n
Definition mod_param.F:479
mod_scalars::sg_kappa
real(dp) sg_kappa
Definition mod_scalars.F:1541
mod_scalars::sg_g
real(dp) sg_g
Definition mod_scalars.F:1540
mod_scalars::sg_z100
real(dp) sg_z100
Definition mod_scalars.F:1545
mod_scalars::sg_n
integer, parameter sg_n
Definition mod_scalars.F:1534
mod_scalars::sg_tol
real(dp) sg_tol
Definition mod_scalars.F:1543
mod_scalars::sg_alpha
real(dp) sg_alpha
Definition mod_scalars.F:1539
mod_scalars::sg_nu
real(dp) sg_nu
Definition mod_scalars.F:1542
mod_scalars::sg_z1min
real(dp) sg_z1min
Definition mod_scalars.F:1546

References bc_2d_mod::bc_r2d_tile(), bc_2d_mod::bc_u2d_tile(), bc_2d_mod::bc_v2d_tile(), mod_scalars::cdb_max, mod_scalars::cdb_min, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::g, mod_sediment::ibwav, mod_sediment::idens, mod_param::inlm, mod_sediment::irhgt, mod_sediment::irlen, mod_sediment::isd50, mod_sediment::izapp, mod_sediment::izdef, mp_exchange_mod::mp_exchange2d(), mod_param::nghostpoints, mod_scalars::nsperiodic, mod_scalars::pi, mod_scalars::sg_alpha, sg_bstress(), mod_scalars::sg_g, mod_scalars::sg_kappa, mod_scalars::sg_n, mod_scalars::sg_nu, sg_purewave(), mod_scalars::sg_tol, mod_scalars::sg_z100, mod_scalars::sg_z1min, and mod_scalars::vonkar.

Here is the call graph for this function:

◆ sg_bstress()

subroutine bbl_mod::sg_bstress ( real(r8), intent(in) sg_row,
real(r8), intent(in) sg_zrozn,
real(r8), intent(in) sg_phicw,
real(r8), intent(in) sg_ubokur,
real(r8), intent(inout) sg_ubouc,
real(r8), intent(out) sg_mu,
real(r8), intent(out) sg_epsilon,
real(r8), intent(out) sg_ro,
real(r8), intent(out) sg_fofx )
private

Definition at line 701 of file sg_bbl.h.

704!
705!=======================================================================
706! !
707! This routine computes bottom stresses via bottom boundary layer !
708! formulation of Styles and Glenn (1999). !
709! !
710! On Input: !
711! !
712! sg_row Ratio of wave excursion amplitude over roughness. !
713! sg_zrozn Ratio of height of current over roughness. !
714! sg_phiwc Angle between wave and currents (radians). !
715! sg_ubokur Ratio of wave over current velocity: !
716! ub/(vonKar*ur) !
717! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
718! (ub/ustarc), first guess. !
719! !
720! On Output: !
721! !
722! sg_ubouc Ratio of bed wave orbital over bottom shear stress !
723! (ub/ustarc), iterated value. !
724! sg_mu Ratio between wave and current bottom shear !
725! stresses (ustarwm/ustarc). !
726! sg_epsilon Ratio between combined (wave and current) and !
727! current bottom shear stresses (ustarc/ustarcw). !
728! sg_ro Internal friction Rossby number: !
729! ustarc/(omega*znot) !
730! sg_fofx Root of PDE used for convergence. !
731! !
732!=======================================================================
733!
734 USE mod_param
735 USE mod_scalars
736!
737! Imported variable declarations.
738!
739 real(r8), intent(in) :: sg_row, sg_zrozn, sg_phicw, sg_ubokur
740
741 real(r8), intent(inout) :: sg_ubouc
742
743 real(r8), intent(out) :: sg_mu, sg_epsilon, sg_ro, sg_fofx
744!
745! Local variable declarations.
746!
747 logical :: ITERATE
748
749 integer :: Iter
750
751 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_cosphi
752 real(r8) :: sg_eps2, sg_kei, sg_keip, sg_ker, sg_kerp, sg_mu2
753 real(r8) :: sg_phi, sg_ror, sg_x, sg_z2p, sg_znotp, sg_zroz1
754 real(r8) :: sg_zroz2, sg_z1ozn, sg_z2ozn
755
756 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
757 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
758 complex(c8) :: sg_ll, sg_nn
759!
760!-----------------------------------------------------------------------
761! Compute bottom stresses.
762!-----------------------------------------------------------------------
763!
764! Compute nondimensional bottom wave shear, phi. Iterate to make
765! sure that there is an upper limit in "ubouc". It usually requires
766! only one pass.
767!
768 iterate=.true.
769 DO iter=1,sg_n
770 IF (iterate) THEN
771 sg_ro=sg_row/sg_ubouc
772 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
773 IF ((sg_z1p/sg_znotp).gt.1.0_r8) THEN
774 sg_x=2.0_r8*sqrt(sg_znotp)
775 IF (sg_x.le.8.0_r8) THEN
776 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
777 & sg_berp, sg_beip, sg_kerp, sg_keip)
778 ELSE
779 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
780 & sg_kerp, sg_keip, sg_berp, sg_beip)
781 END IF
782 cff=1.0_r8/sqrt(sg_znotp)
783 sg_bnot =cmplx(sg_ber,sg_bei)
784 sg_knot =cmplx(sg_ker,sg_kei)
785 sg_bnotp=cmplx(sg_berp,sg_beip)*cff
786 sg_knotp=cmplx(sg_kerp,sg_keip)*cff
787!
788 sg_x=2.0_r8*sqrt(sg_z1p)
789 IF (sg_x.le.8.0_r8) THEN
790 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
791 & sg_berp, sg_beip, sg_kerp, sg_keip)
792 ELSE
793 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
794 & sg_kerp, sg_keip, sg_berp, sg_beip)
795 END IF
796 cff=1.0_r8/sqrt(sg_z1p)
797 sg_b1 =cmplx(sg_ber,sg_bei)
798 sg_k1 =cmplx(sg_ker,sg_kei)
799 sg_b1p=cmplx(sg_berp,sg_beip)*cff
800 sg_k1p=cmplx(sg_kerp,sg_keip)*cff
801!
802 sg_ll=sg_mp*sg_b1+sg_b1p
803 sg_nn=sg_mp*sg_k1+sg_k1p
804 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
805 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
806 sg_gammai=-sg_kappa*sg_znotp*sg_argi
807 sg_phi=cabs(sg_gammai)
808 ELSE
809 sg_gammai=-sg_kappa*sg_z1p*sg_mp
810 sg_phi=cabs(sg_gammai)
811 END IF
812!
813 IF (sg_ubouc.gt.(1.0_r8/sg_phi)) THEN
814 sg_ubouc=1.0_r8/sg_phi
815 ELSE
816 iterate=.false.
817 END IF
818 END IF
819 END DO
820!
821! Compute ratio of wave over current bottom shear stresses.
822!
823 sg_mu=sqrt(sg_ubouc*sg_phi)
824!
825! Compute ratio of current over combined bottom shear stresses.
826!
827 IF (sg_mu.eq.1.0_r8) THEN
828 sg_epsilon=0.0_r8
829 ELSE
830 sg_mu2=sg_mu*sg_mu
831 sg_cosphi=abs(cos(sg_phicw))
832 sg_eps2=-sg_mu2*sg_cosphi+ &
833 & sqrt(1.0_r8+sg_mu2*sg_mu2*(sg_cosphi*sg_cosphi-1.0_r8))
834 sg_epsilon=sqrt(sg_eps2)
835 END IF
836!
837! Determine root of PDE used for convergence.
838!
839 IF (sg_epsilon.ne.0.0_r8) THEN
840 sg_z2p=sg_z1p/sg_epsilon
841 sg_ror=sg_ro/sg_zrozn
842 sg_zroz1=1.0_r8/(sg_alpha*sg_kappa*sg_ror)
843 sg_zroz2=sg_epsilon*sg_zroz1
844 sg_z1ozn=sg_alpha*sg_kappa*sg_ro
845 sg_z2ozn=sg_z1ozn/sg_epsilon
846!
847 IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
848 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
849 & (log(sg_zroz2)+1.0_r8-sg_epsilon+ &
850 & sg_epsilon*log(sg_z1ozn))
851 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
852 & (sg_z1ozn.gt.1.0_r8)) THEN
853 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
854 & (sg_zroz1-1.0_r8+log(sg_z1ozn))
855 ELSE IF ((sg_zroz1.le.1.0_r8).and.(sg_z1ozn.gt.1.0_r8)) THEN
856 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
857 & log(sg_zrozn)
858 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z1ozn.le.1.0_r8).and. &
859 & (sg_z2ozn.gt.1.0_r8)) THEN
860 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon* &
861 & (log(sg_zroz2)+1.0_r8-1.0_r8/sg_z2ozn)
862 ELSE IF ((sg_zroz2.le.1.0_r8).and.(sg_zroz1.gt.1.0_r8).and. &
863 & (sg_z1ozn.le.1.0_r8).and.(sg_z2ozn.gt.1.0_r8)) THEN
864 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*sg_epsilon* &
865 & (sg_zroz1-1.0_r8/sg_z1ozn)
866 ELSE IF ((sg_zroz2.gt.1.0_r8).and.(sg_z2ozn.le.1.0_r8)) THEN
867 sg_fofx=-sg_ubouc+sg_ubokur*sg_epsilon*log(sg_zrozn)
868 END IF
869 END IF
870!
871 RETURN
mod_scalars::sg_mp
complex(c8) sg_mp
Definition mod_scalars.F:1551
mod_scalars::sg_z1p
real(dp) sg_z1p
Definition mod_scalars.F:1547

References mod_scalars::sg_alpha, mod_scalars::sg_kappa, sg_kelvin8m(), sg_kelvin8p(), mod_scalars::sg_mp, mod_scalars::sg_n, and mod_scalars::sg_z1p.

Referenced by sg_bbl_tile(), and ssw_bbl_tile().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ sg_kelvin8m()

subroutine bbl_mod::sg_kelvin8m ( real(r8), intent(in) x,
real(r8), intent(out) ber,
real(r8), intent(out) bei,
real(r8), intent(out) ker,
real(r8), intent(out) kei,
real(r8), intent(out) berp,
real(r8), intent(out) beip,
real(r8), intent(out) kerp,
real(r8), intent(out) keip )
private

Definition at line 979 of file sg_bbl.h.

981!
982!=======================================================================
983! !
984! This rotuine computes the Kelvin functions for arguments less !
985! than eight. !
986! !
987!=======================================================================
988!
989 USE mod_scalars
990!
991! Imported variable declarations.
992!
993 real(r8), intent(in) :: x
994 real(r8), intent(out) :: ber, bei, ker, kei
995 real(r8), intent(out) :: berp, beip, kerp, keip
996!
997! Local variable declarations.
998!
999 integer :: i
1000
1001 real(r8) :: cff, xhalf
1002
1003 real(r8), dimension(28) :: xp
1004!
1005!-----------------------------------------------------------------------
1006! Compute Kelvin functions.
1007!-----------------------------------------------------------------------
1008!
1009 cff=0.125_r8*x
1010 xp(1)=cff
1011 DO i=2,28
1012 xp(i)=xp(i-1)*cff
1013 END DO
1014 xhalf=0.5_r8*x
1015!
1016 ber=1.0_r8- &
1017 & 64.0_r8*xp(4)+113.77777774_r8*xp(8)- &
1018 & 32.36345652_r8*xp(12)+2.64191397_r8*xp(16)- &
1019 & 0.08349609_r8*xp(20)+0.00122552_r8*xp(24)- &
1020 & 0.00000901_r8*xp(28)
1021 bei=16.0_r8*xp(2)-113.77777774_r8*xp(6)+ &
1022 & 72.81777742*xp(10)-10.56765779_r8*xp(14)+ &
1023 & 0.52185615_r8*xp(18)-0.01103667_r8*xp(22)+ &
1024 & 0.00011346*xp(26)
1025!
1026 ker=-ber*log(xhalf)+0.25_r8*pi*bei- &
1027 & 0.57721566_r8-59.05819744*xp(4)+ &
1028 & 171.36272133_r8*xp(8)-60.60977451_r8*xp(12)+ &
1029 & 5.65539121_r8*xp(16)-0.19636347_r8*xp(20)+ &
1030 & 0.00309699_r8*xp(24)-0.00002458_r8*xp(28)
1031 kei=-bei*log(xhalf)-0.25_r8*pi*ber+ &
1032 & 6.76454936_r8*xp(2)-142.91827687_r8*xp(6)+ &
1033 & 124.23569650_r8*xp(10)-21.30060904_r8*xp(14)+ &
1034 & 1.17509064_r8*xp(18)-0.02695875_r8*xp(22)+ &
1035 & 0.00029532_r8*xp(26)
1036!
1037 berp=x*(-4.0_r8*xp(2)+14.22222222_r8*xp(6)- &
1038 & 6.06814810_r8*xp(10)+0.66047849_r8*xp(14)- &
1039 & 0.02609253_r8*xp(18)+0.00045957_r8*xp(22)- &
1040 & 0.00000394_r8*xp(26))
1041 beip=x*(0.5_r8-10.66666666_r8*xp(4)+11.37777772_r8*xp(8)- &
1042 & 2.31167514_r8*xp(12)+0.14677204_r8*xp(16)- &
1043 & 0.00379386_r8*xp(20)+0.00004609_r8*xp(24))
1044!
1045 kerp=-berp*log(xhalf)-ber/x+0.25*pi*beip+ &
1046 & x*(-3.69113734_r8*xp(2)+21.42034017_r8*xp(6)- &
1047 & 11.36433272_r8*xp(10)+1.41384780_r8*xp(14)- &
1048 & 0.06136358_r8*xp(18)+0.00116137_r8*xp(22)- &
1049 & 0.00001075*xp(26))
1050 keip=-beip*log(xhalf)-bei/x-0.25_r8*pi*berp+ &
1051 & x*(0.21139217_r8-13.39858846_r8*xp(4)+ &
1052 & 19.41182758_r8*xp(8)-4.65950823_r8*xp(12)+ &
1053 & 0.33049424_r8*xp(16)-0.00926707_r8*xp(20)+ &
1054 & 0.00011997_r8*xp(24))
1055!
1056 RETURN

References mod_scalars::pi.

Referenced by sg_bstress(), sg_purewave(), and ssw_bbl_tile().

Here is the caller graph for this function:

◆ sg_kelvin8p()

subroutine bbl_mod::sg_kelvin8p ( real(r8), intent(in) x,
real(r8), intent(out) ker,
real(r8), intent(out) kei,
real(r8), intent(out) ber,
real(r8), intent(out) bei,
real(r8), intent(out) kerp,
real(r8), intent(out) keip,
real(r8), intent(out) berp,
real(r8), intent(out) beip )
private

Definition at line 1059 of file sg_bbl.h.

1061!
1062!=======================================================================
1063! !
1064! This rotuine computes the Kelvin functions for arguments greater !
1065! than eight. !
1066! !
1067!=======================================================================
1068!
1069 USE mod_scalars
1070!
1071! Imported variable declarations.
1072!
1073 real(r8), intent(in) :: x
1074 real(r8), intent(out) :: ker, kei, ber, bei
1075 real(r8), intent(out) :: kerp, keip, berp, beip
1076!
1077! Local variable declarations.
1078!
1079 integer :: i
1080
1081 real(r8) :: cff, xhalf
1082
1083 real(r8), dimension(6) :: xm, xp
1084
1085 complex(c8) :: argm, argp, fofx, gofx, phim, phip, thetam, thetap
1086!
1087!-----------------------------------------------------------------------
1088! Compute Kelvin functions.
1089!-----------------------------------------------------------------------
1090!
1091 cff=8.0_r8/x
1092 xp(1)=cff
1093 xm(1)=-cff
1094 DO i=2,6
1095 xp(i)=xp(i-1)*cff
1096 xm(i)=-xm(i-1)*cff
1097 END DO
1098!
1099 thetap=cmplx(0.0_r8,-0.3926991_r8)+ &
1100 & cmplx(0.0110486_r8,-0.0110485_r8)*xp(1)+ &
1101 & cmplx(0.0_r8,-0.0009765_r8)*xp(2)+ &
1102 & cmplx(-0.0000906_r8,-0.0000901_r8)*xp(3)+ &
1103 & cmplx(-0.0000252_r8,0.0_r8)*xp(4)+ &
1104 & cmplx(-0.0000034_r8,0.0000051_r8)*xp(5)+ &
1105 & cmplx(0.0000006,0.0000019)*xp(6)
1106 thetam=cmplx(0.0_r8,-0.3926991_r8)+ &
1107 & cmplx(0.0110486_r8,-0.0110485_r8)*xm(1)+ &
1108 & cmplx(0.0_r8,-0.0009765_r8)*xm(2)+ &
1109 & cmplx(-0.0000906_r8,-0.0000901_r8)*xm(3)+ &
1110 & cmplx(-0.0000252_r8,0.0_r8)*xm(4)+ &
1111 & cmplx(-0.0000034_r8,0.0000051_r8)*xm(5)+ &
1112 & cmplx(0.0000006_r8,0.0000019_r8)*xm(6)
1113!
1114 phip=cmplx(0.7071068_r8,0.7071068_r8)+ &
1115 & cmplx(-0.0625001_r8,-0.0000001_r8)*xp(1)+ &
1116 & cmplx(-0.0013813_r8,0.0013811_r8)*xp(2)+ &
1117 & cmplx(0.0000005_r8,0.0002452_r8)*xp(3)+ &
1118 & cmplx(0.0000346_r8,0.0000338_r8)*xp(4)+ &
1119 & cmplx(0.0000117_r8,-0.0000024_r8)*xp(5)+ &
1120 & cmplx(0.0000016_r8,-0.0000032_r8)*xp(6)
1121 phim=cmplx(0.7071068_r8,0.7071068_r8)+ &
1122 & cmplx(-0.0625001_r8,-0.0000001_r8)*xm(1)+ &
1123 & cmplx(-0.0013813_r8,0.0013811_r8)*xm(2)+ &
1124 & cmplx(0.0000005_r8,0.0002452_r8)*xm(3)+ &
1125 & cmplx(0.0000346_r8,0.0000338_r8)*xm(4)+ &
1126 & cmplx(0.0000117_r8,-0.0000024_r8)*xm(5)+ &
1127 & cmplx(0.0000016_r8,-0.0000032_r8)*xm(6)
1128!
1129 cff=x/sqrt(2.0_r8)
1130 argm=-cff*cmplx(1.0_r8,1.0_r8)+thetam
1131 fofx=sqrt(pi/(2.0_r8*x))*cexp(argm)
1132 ker=real(fofx)
1133 kei=aimag(fofx)
1134!
1135 argp=cff*cmplx(1.0_r8,1.0_r8)+thetap
1136 gofx=1.0_r8/sqrt(2.0_r8*pi*x)*cexp(argp)
1137 ber=real(gofx)-kei/pi
1138 bei=aimag(gofx)+ker/pi
1139!
1140 kerp=real(-fofx*phim)
1141 keip=aimag(-fofx*phim)
1142!
1143 berp=real(gofx*phip)-keip/pi
1144 beip=aimag(gofx*phip)+kerp/pi
1145!
1146 RETURN

References mod_scalars::pi.

Referenced by sg_bstress(), sg_purewave(), and ssw_bbl_tile().

Here is the caller graph for this function:

◆ sg_purewave()

subroutine bbl_mod::sg_purewave ( real(r8), intent(in) sg_row,
real(r8), intent(inout) sg_ubouwm,
real(r8), intent(out) sg_znotp,
real(r8), intent(out) sg_ro )
private

Definition at line 874 of file sg_bbl.h.

875!
876!=======================================================================
877! !
878! This routine determines the maximum ratio of waves over combined !
879! bottom shear stress. !
880! !
881! On Input: !
882! !
883! sg_row Ratio of wave excursion amplitude over roughness. !
884! sg_ubouwm Maximum ratio of waves over combined bottom shear !
885! stress. !
886! !
887! On Output: !
888! !
889! sg_ubouwm Maximum ratio of waves over combined bottom shear !
890! stress. !
891! sg_znotp Ratio of hydraulic roughness over scaled height !
892! of bottom boundary layer. !
893! sg_ro Internal friction Rossby number: !
894! ustarc/(omega*znot) !
895! !
896!=======================================================================
897!
898 USE mod_param
899 USE mod_scalars
900!
901! Imported variable declarations.
902!
903 real(r8), intent(in) :: sg_row
904
905 real(r8), intent(inout) :: sg_ubouwm
906
907 real(r8), intent(out) :: sg_znotp, sg_ro
908!
909! Local variable declarations.
910!
911 integer :: Iter
912
913 real(r8) :: cff, sg_bei, sg_beip, sg_ber, sg_berp, sg_kei
914 real(r8) :: sg_keip, sg_ker, sg_kerp, sg_phi, sg_ubouwmn, sg_x
915
916 complex(c8) :: sg_argi, sg_bnot, sg_bnotp, sg_b1, sg_b1p
917 complex(c8) :: sg_gammai, sg_knot, sg_knotp, sg_k1, sg_k1p
918 complex(c8) :: sg_ll, sg_nn
919!
920!-----------------------------------------------------------------------
921! Compute wind-induced wave stress.
922!-----------------------------------------------------------------------
923!
924 DO iter=1,sg_n
925 sg_ro=sg_row/sg_ubouwm
926 sg_znotp=1.0_r8/(sg_kappa*sg_ro)
927 IF (sg_z1p/sg_znotp.gt.1.0_r8) THEN
928 sg_x=2.0_r8*sqrt(sg_znotp)
929 IF (sg_x.le.8.0_r8) THEN
930 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
931 & sg_berp, sg_beip, sg_kerp, sg_keip)
932 ELSE
933 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
934 & sg_kerp, sg_keip, sg_berp, sg_beip)
935 END IF
936 cff=1.0_r8/sqrt(sg_znotp)
937 sg_bnot =cmplx(sg_ber,sg_bei)
938 sg_knot =cmplx(sg_ker,sg_kei)
939 sg_bnotp=cmplx(sg_berp,sg_beip)*cff
940 sg_knotp=cmplx(sg_kerp,sg_keip)*cff
941!
942 sg_x=2.0*sqrt(sg_z1p)
943 IF (sg_x.le.8.0_r8) THEN
944 CALL sg_kelvin8m (sg_x, sg_ber, sg_bei, sg_ker, sg_kei, &
945 & sg_berp, sg_beip, sg_kerp, sg_keip)
946 ELSE
947 CALL sg_kelvin8p (sg_x, sg_ker, sg_kei, sg_ber, sg_bei, &
948 & sg_kerp, sg_keip, sg_berp, sg_beip)
949 END IF
950 cff=1.0_r8/sqrt(sg_z1p)
951 sg_b1 =cmplx(sg_ber,sg_bei)
952 sg_k1 =cmplx(sg_ker,sg_kei)
953 sg_b1p=cmplx(sg_berp,sg_beip)*cff
954 sg_k1p=cmplx(sg_kerp,sg_keip)*cff
955!
956 sg_ll=sg_mp*sg_b1+sg_b1p
957 sg_nn=sg_mp*sg_k1+sg_k1p
958 sg_argi=sg_bnotp*sg_nn/(sg_bnot*sg_nn-sg_knot*sg_ll)+ &
959 & sg_knotp*sg_ll/(sg_knot*sg_ll-sg_bnot*sg_nn)
960 sg_gammai=-sg_kappa*sg_znotp*sg_argi
961 sg_phi=cabs(sg_gammai)
962 ELSE
963 sg_gammai=-sg_kappa*sg_z1p*sg_mp
964 sg_phi=cabs(sg_gammai)
965 END IF
966!
967 sg_ubouwmn=1.0_r8/sg_phi
968 IF (abs((sg_ubouwmn-sg_ubouwm)/sg_ubouwmn).le.sg_tol) THEN
969 sg_ubouwm=sg_ubouwmn
970 RETURN
971 ELSE
972 sg_ubouwm=sg_ubouwmn
973 END IF
974 END DO
975!
976 RETURN

References mod_scalars::sg_kappa, sg_kelvin8m(), sg_kelvin8p(), mod_scalars::sg_mp, mod_scalars::sg_n, mod_scalars::sg_tol, and mod_scalars::sg_z1p.

Referenced by sg_bbl_tile(), and ssw_bbl_tile().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ ssw_bbl_tile()

subroutine bbl_mod::ssw_bbl_tile ( integer, intent(in) ng,
integer, intent(in) tile,
integer, intent(in) lbi,
integer, intent(in) ubi,
integer, intent(in) lbj,
integer, intent(in) ubj,
integer, intent(in) imins,
integer, intent(in) imaxs,
integer, intent(in) jmins,
integer, intent(in) jmaxs,
integer, intent(in) nrhs,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) h,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) z_r,
real(r8), dimension(lbi:ubi,lbj:ubj,0:n(ng)), intent(in) z_w,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) angler,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) zobot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) hwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) uwave_rms,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) dwave,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(in) pwave_bot,
real(r8), dimension(lbi:ubi,lbj:ubj,1:nst), intent(in) bedldu,
real(r8), dimension(lbi:ubi,lbj:ubj,1:nst), intent(in) bedldv,
real(r8), dimension(lbi:ubi,lbj:ubj,mbotp), intent(inout) bottom,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) rho,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) u,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng),2), intent(in) v,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) u_stokes,
real(r8), dimension(lbi:ubi,lbj:ubj,n(ng)), intent(in) v_stokes,
real(r8), dimension(lbi:ubi,lbj:ubj,3), intent(in) zeta,
integer, dimension(lbi:ubi,lbj:ubj), intent(inout) iconv,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ubot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vbot,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) ur,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) vr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrc,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrw,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstrcwmax,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bustr,
real(r8), dimension(lbi:ubi,lbj:ubj), intent(out) bvstr )
private

Definition at line 113 of file ssw_bbl.h.

141!***********************************************************************
142!
143 USE mod_param
144 USE mod_parallel
145 USE mod_iounits
146 USE mod_scalars
147 USE mod_sediment
148!
149 USE bc_2d_mod
150#ifdef DISTRIBUTE
151 USE mp_exchange_mod, ONLY : mp_exchange2d
152#endif
153
154!
155! Imported variable declarations.
156!
157 integer, intent(in) :: ng, tile
158 integer, intent(in) :: LBi, UBi, LBj, UBj
159 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
160 integer, intent(in) :: nrhs
161
162#ifdef ASSUMED_SHAPE
163 integer, intent(inout) :: Iconv(LBi:,LBj:)
164
165 real(r8), intent(in) :: h(LBi:,LBj:)
166 real(r8), intent(in) :: z_r(LBi:,LBj:,:)
167 real(r8), intent(in) :: z_w(LBi:,LBj:,0:)
168 real(r8), intent(in) :: angler(LBi:,LBj:)
169 real(r8), intent(in) :: ZoBot(LBi:,LBj:)
170# if defined SSW_CALC_UB
171 real(r8), intent(in) :: Hwave(LBi:,LBj:)
172# else
173 real(r8), intent(in) :: Uwave_rms(LBi:,LBj:)
174# endif
175 real(r8), intent(in) :: Dwave(LBi:,LBj:)
176 real(r8), intent(in) :: Pwave_bot(LBi:,LBj:)
177# ifdef BEDLOAD
178 real(r8), intent(in) :: bedldu(LBi:,LBj:,:)
179 real(r8), intent(in) :: bedldv(LBi:,LBj:,:)
180# endif
181 real(r8), intent(inout) :: bottom(LBi:,LBj:,:)
182 real(r8), intent(in) :: rho(LBi:,LBj:,:)
183 real(r8), intent(in) :: u(LBi:,LBj:,:,:)
184 real(r8), intent(in) :: v(LBi:,LBj:,:,:)
185# if defined SSW_LOGINT_STOKES
186 real(r8), intent(in) :: u_stokes(LBi:,LBj:,:)
187 real(r8), intent(in) :: v_stokes(LBi:,LBj:,:)
188# endif
189# if defined SSW_CALC_UB
190 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
191# endif
192 real(r8), intent(out) :: Ubot(LBi:,LBj:)
193 real(r8), intent(out) :: Vbot(LBi:,LBj:)
194 real(r8), intent(out) :: Ur(LBi:,LBj:)
195 real(r8), intent(out) :: Vr(LBi:,LBj:)
196 real(r8), intent(out) :: bustrc(LBi:,LBj:)
197 real(r8), intent(out) :: bvstrc(LBi:,LBj:)
198 real(r8), intent(out) :: bustrw(LBi:,LBj:)
199 real(r8), intent(out) :: bvstrw(LBi:,LBj:)
200 real(r8), intent(out) :: bustrcwmax(LBi:,LBj:)
201 real(r8), intent(out) :: bvstrcwmax(LBi:,LBj:)
202 real(r8), intent(out) :: bustr(LBi:,LBj:)
203 real(r8), intent(out) :: bvstr(LBi:,LBj:)
204#else
205 integer, intent(inout) :: Iconv(LBi:UBi,LBj:UBj)
206
207 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
208 real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
209 real(r8), intent(in) :: z_w(LBi:UBi,LBj:UBj,0:N(ng))
210 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
211 real(r8), intent(in) :: ZoBot(LBi:UBi,LBj:UBj)
212# if defined SSW_CALC_UB
213 real(r8), intent(in) :: Hwave(LBi:UBi,LBj:UBj)
214# else
215 real(r8), intent(in) :: Uwave_rms(LBi:UBi,LBj:UBj)
216# endif
217 real(r8), intent(in) :: Dwave(LBi:UBi,LBj:UBj)
218 real(r8), intent(in) :: Pwave_bot(LBi:UBi,LBj:UBj)
219# ifdef BEDLOAD
220 real(r8), intent(in) :: bedldu(LBi:UBi,LBj:UBj,1:NST)
221 real(r8), intent(in) :: bedldv(LBi:UBi,LBj:UBj,1:NST)
222# endif
223 real(r8), intent(inout) :: bottom(LBi:UBi,LBj:UBj,MBOTP)
224 real(r8), intent(in) :: rho(LBi:UBi,LBj:UBj,N(ng))
225 real(r8), intent(in) :: u(LBi:UBi,LBj:UBj,N(ng),2)
226 real(r8), intent(in) :: v(LBi:UBi,LBj:UBj,N(ng),2)
227# if defined SSW_LOGINT_STOKES
228 real(r8), intent(in) :: u_stokes(LBi:UBi,LBj:UBj,N(ng))
229 real(r8), intent(in) :: v_stokes(LBi:UBi,LBj:UBj,N(ng))
230# endif
231# if defined SSW_CALC_UB
232 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
233# endif
234 real(r8), intent(out) :: Ubot(LBi:UBi,LBj:UBj)
235 real(r8), intent(out) :: Vbot(LBi:UBi,LBj:UBj)
236 real(r8), intent(out) :: Ur(LBi:UBi,LBj:UBj)
237 real(r8), intent(out) :: Vr(LBi:UBi,LBj:UBj)
238 real(r8), intent(out) :: bustrc(LBi:UBi,LBj:UBj)
239 real(r8), intent(out) :: bvstrc(LBi:UBi,LBj:UBj)
240 real(r8), intent(out) :: bustrw(LBi:UBi,LBj:UBj)
241 real(r8), intent(out) :: bvstrw(LBi:UBi,LBj:UBj)
242 real(r8), intent(out) :: bustrcwmax(LBi:UBi,LBj:UBj)
243 real(r8), intent(out) :: bvstrcwmax(LBi:UBi,LBj:UBj)
244 real(r8), intent(out) :: bustr(LBi:UBi,LBj:UBj)
245 real(r8), intent(out) :: bvstr(LBi:UBi,LBj:UBj)
246#endif
247!
248! Local variable declarations.
249!
250 logical :: ITERATE
251
252 integer :: Iter, i, j, k
253
254 real(r8), parameter :: eps = 1.0e-10_r8
255
256 real(r8) :: Kbh, Kbh2, Kdh
257 real(r8) :: taucr, wsedr, tstar, coef_st
258 real(r8) :: coef_b1, coef_b2, coef_b3, d0
259 real(r8) :: dolam, dolam1, doeta1, doeta2, fdo_etaano
260 real(r8) :: lamorb, lamanorb
261 real(r8) :: m_ubr, m_wr, m_ucr, m_zr, m_phicw, m_kb
262 real(r8) :: m_ustrc, m_ustrwm, m_ustrr, m_fwc, m_zoa, m_dwc
263 real(r8) :: zo, Dstp
264 real(r8) :: Kb, Kdelta, Ustr
265 real(r8) :: anglec, anglew
266 real(r8) :: cff, cff1, cff2, cff3, og, fac, fac1, fac2
267 real(r8) :: sg_ab, sg_abokb, sg_a1, sg_b1, sg_chi, sg_c1, d50
268 real(r8) :: sg_epsilon, ssw_eta, sg_fofa, sg_fofb, sg_fofc, sg_fwm
269 real(r8) :: sg_kbs, ssw_lambda, sg_mu, sg_phicw, sg_ro, sg_row
270 real(r8) :: sg_shdnrm, sg_shld, sg_shldcr, sg_scf, rhos, sg_star
271 real(r8) :: sg_ub, sg_ubokur, sg_ubouc, sg_ubouwm, sg_ur
272 real(r8) :: sg_ustarc, sg_ustarcw, sg_ustarwm, sg_znot, sg_znotp
273 real(r8) :: sg_zr, sg_zrozn, sg_z1, sg_z1ozn, sg_z2, z1, z2
274 real(r8) :: zoMIN, zoMAX
275 real(r8) :: coef_fd
276
277 real(r8), parameter :: twopi=2.0_r8*pi
278!
279 real(r8), parameter :: absolute_zoMIN = 5.0d-5 ! in Harris-Wiberg
280!! real(r8), parameter :: absolute_zoMIN = 5.0d-8 ! in Harris-Wiberg
281 real(r8), parameter :: Cd_fd = 0.5_r8
282
283 real(r8), parameter :: K1 = 0.6666666666_r8 ! Coefficients for
284 real(r8), parameter :: K2 = 0.3555555555_r8 ! explicit
285 real(r8), parameter :: K3 = 0.1608465608_r8 ! wavenumber
286 real(r8), parameter :: K4 = 0.0632098765_r8 ! calculation
287 real(r8), parameter :: K5 = 0.0217540484_r8 ! (Dean and
288 real(r8), parameter :: K6 = 0.0065407983_r8 ! Dalrymple, 1991)
289
290 real(r8), parameter :: coef_a1=0.095_r8 ! Coefficients for
291 real(r8), parameter :: coef_a2=0.442_r8 ! ripple predictor
292 real(r8), parameter :: coef_a3=2.280_r8 ! (Wiberg-Harris)
293
294#if defined GM82_RIPRUF
295 real(r8), parameter :: ar = 27.7_r8/30.0_r8 ! Grant-Madsen (1982)
296#elif defined N92_RIPRUF
297 real(r8), parameter :: ar = 0.267_r8 ! Nielsen (1992)
298#elif defined R88_RIPRUF
299 real(r8), parameter :: ar = 0.533_r8 ! Raudkivi (1988)
300#else
301 no ripple roughness coeff. chosen
302#endif
303
304 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ab
305 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Fwave_bot
306 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauc
307 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Tauw
308 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Taucwmax
309 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sg
310 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sg
311 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ub
312 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ucur
313 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Umag
314 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vcur
315 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr
316 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic
317 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phicw
318 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rheight
319 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: rlength
320 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: u100
321 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znot
322 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: znotc
323 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoN
324 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoST
325 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoBF
326 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoDEF
327 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: zoBIO
328 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: thck_wbl
329#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
330 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ksd_wbl
331 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: ustrc_wbl
332#endif
333#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
334 defined bedload_vandera_direct_udelta
335 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: udelta_wbl
336 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Zr_wbl
337 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: phic_sgwbl
338#endif
339!
340 real(r8), dimension(1:N(ng)) :: Urz, Vrz
341#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
342 defined bedload_vandera_direct_udelta
343 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Ur_sgwbl
344 real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: Vr_sgwbl
345 real(r8) :: Ucur_sgwbl, Vcur_sgwbl
346#endif
347!
348#include "set_bounds.h"
349!
350!-----------------------------------------------------------------------
351! Set currents above the bed.
352!-----------------------------------------------------------------------
353!
354 DO j=jstrv-1,jend
355 DO i=istru-1,iend
356!
357! Calculate bottom cell thickness
358!
359 zr(i,j)=z_r(i,j,1)-z_w(i,j,0)
360!
361#if defined SSW_LOGINT
362!
363 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
364!
365# if defined SSW_LOGINT_DIRECT
366!
367! Cap minimum Zr 0.9*depth and user input.
368!
369 cff1=min(0.9_r8*dstp, max(zr(i,j), sg_zwbl(ng)))
370# elif defined SSW_LOGINT_WBL
371!
372! Cap minimum Zr 0.9*depth and computed wave bound layer
373!
374 cff1=min(0.98_r8*dstp,max(zr(i,j), bottom(i,j,idtbl)*1.1_r8))
375# else
376!
377! Use the original coded ssw_logint formulation
378!
379 cff1=sg_z1min
380# endif
381!
382! If using the logarithmic interpolation
383!
384 DO k=1,n(ng)
385 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
386 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
387# ifdef SSW_LOGINT_STOKES
388 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
389 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
390# endif
391 END DO
392 CALL log_interp( n(ng), dstp, cff1, &
393 & urz, vrz, &
394 & z_r(i,j,:), z_w(i,j,:), &
395 & bottom(i,j,isd50), bottom(i,j,izapp), &
396 & zr(i,j), &
397 & ur_sg(i,j), vr_sg(i,j) )
398#else
399!
400! Regular method to get reference velocity for Madsen
401! without using log interp by the bottom cell.
402!
403 ur_sg(i,j)=0.5_r8*(u(i,j,1,nrhs)+u(i+1,j,1,nrhs))
404 vr_sg(i,j)=0.5_r8*(v(i,j,1,nrhs)+v(i,j+1,1,nrhs))
405# ifdef SSW_LOGINT_STOKES
406 ur_sg(i,j)=ur_sg(i,j)+0.5_r8*(u_stokes(i,j,1)+u_stokes(i+1,j,1))
407 vr_sg(i,j)=vr_sg(i,j)+0.5_r8*(v_stokes(i,j,1)+v_stokes(i,j+1,1))
408# endif
409#endif
410 END DO
411 END DO
412!
413!-----------------------------------------------------------------------
414! Compute bottom stresses.
415!-----------------------------------------------------------------------
416!
417 DO j=jstrv-1,jend
418 DO i=istru-1,iend
419!
420! Set bed wave orbital velocity and excursion amplitude. Use data
421! from wave models (SWAN) or use Dean and Dalrymple (1991) 6th-degree
422! polynomial to approximate wave number on shoaling water.
423
424 fwave_bot(i,j)=twopi/max(pwave_bot(i,j),1.0_r8)
425#ifdef SSW_CALC_UB
426 kdh=(h(i,j)+zeta(i,j,nrhs))*fwave_bot(i,j)**2/g
427 kbh2=kdh*kdh+ &
428 & kdh/(1.0_r8+kdh*(k1+kdh*(k2+kdh*(k3+kdh*(k4+ &
429 & kdh*(k5+k6*kdh))))))
430 kbh=sqrt(kbh2)
431 ab(i,j)=0.5_r8*hwave(i,j)/sinh(kbh)+eps
432 ub(i,j)=fwave_bot(i,j)*ab(i,j)+eps
433#else
434 ub(i,j)=max(uwave_rms(i,j),0.0_r8)+eps
435 ab(i,j)=ub(i,j)/fwave_bot(i,j)+eps
436#endif
437!
438! Compute bottom current magnitude at RHO-points.
439!
440 ucur(i,j)=ur_sg(i,j)
441 vcur(i,j)=vr_sg(i,j)
442!
443 umag(i,j)=sqrt(ucur(i,j)*ucur(i,j)+vcur(i,j)*vcur(i,j)+eps)
444!
445! Compute angle between currents and waves (radians)
446!
447 IF (ucur(i,j).eq.0.0_r8) THEN
448 phic(i,j)=0.5_r8*pi*sign(1.0_r8,vcur(i,j))
449 ELSE
450 phic(i,j)=atan2(vcur(i,j),ucur(i,j))
451 ENDIF
452 phicw(i,j)=1.5_r8*pi-dwave(i,j)-phic(i,j)-angler(i,j)
453 END DO
454 END DO
455!
456! Loop over RHO points.
457!
458 DO j=jstrv-1,jend
459 DO i=istru-1,iend
460!
461! Load sediment properties, stresses, and roughness from previous time
462! step (stresses are in m2 s-2).
463!
464 d50=bottom(i,j,isd50)
465 rhos=bottom(i,j,idens)/(rho(i,j,1)+1000.0_r8)
466 wsedr=bottom(i,j,iwsed)
467 taucr=bottom(i,j,itauc)
468 tauc(i,j)=sqrt(bustrc(i,j)**2+bvstrc(i,j)**2)
469 tauw(i,j)=sqrt(bustrw(i,j)**2+bvstrw(i,j)**2)
470 taucwmax(i,j)=sqrt( bustrcwmax(i,j)**2+bvstrcwmax(i,j)**2)
471!
472 rheight(i,j)=bottom(i,j,irhgt)
473 rlength(i,j)=bottom(i,j,irlen)
474 zomax=0.9_r8*zr(i,j)
475 zomin=max(absolute_zomin,2.5_r8*d50/30.0_r8)
476!
477! Initialize arrays.
478!
479 zon(i,j)=min(max(2.5_r8*d50/30.0_r8, zomin ),zomax)
480 zost(i,j)=0.0_r8
481 zobf(i,j)=0.0_r8
482 zobio(i,j)=0.0_r8
483!! zoDEF(i,j)=bottom(i,j,izdef)
484 zodef(i,j)=zobot(i,j)
485
486#ifdef SSW_CALC_ZNOT
487!
488! Calculate components of roughness and sum: zo = zoN + zoST + zoBF
489! Determine whether sediment is in motion. Use Shields criterion to
490! determine if sediment is mobile.
491!
492 tstar=taucwmax(i,j)/(taucr+eps)
493 IF (tstar.lt.1.0_r8) THEN ! no motion
494 zost(i,j)=0.0_r8
495 zobf(i,j)=ar*rheight(i,j)**2/(rlength(i,j)+eps)
496 ELSE
497!
498! Threshold of motion exceeded - calculate new zoST and zoBF
499! Calculate saltation roughness according to Wiberg & Rubin (1989)
500! (Eqn. 11 in Harris & Wiberg, 2001)
501! (d50 is in m, but this formula needs cm)
502!
503 coef_st=0.0204_r8*log(100.0_r8*d50+eps)**2+ &
504 & 0.0220_r8*log(100.0_r8*d50+eps)+0.0709_r8
505 zost(i,j)=0.056_r8*d50*0.68_r8*tstar/ &
506 & (1.0_r8+coef_st*tstar)
507 IF (zost(i,j).lt.0.0_r8) THEN
508 IF (master) THEN
509 WRITE (stdout,'(/,a)') &
510 & 'Warning zoST < 0: tstar, d50, coef_st:'
511 WRITE (stdout,*) tstar, d50, coef_st
512 END IF
513 END IF
514!
515! Calculate ripple height and wavelength.
516! Use Malarkey and Davies (2003) explict version of Wiberg and Harris.
517!
518 coef_b1=1.0_r8/coef_a1
519 coef_b2=0.5_r8*(1.0_r8 + coef_a2)*coef_b1
520 coef_b3=coef_b2**2-coef_a3*coef_b1
521 d0=2.0_r8*ab(i,j)
522 IF ((d0/d50).gt.13000.0_r8) THEN ! sheet flow
523 rheight(i,j)=0.0_r8
524 rlength(i,j)=535.0_r8*d50 ! does not matter since
525 ELSE ! rheight=0
526 dolam1=d0/(535.0_r8*d50)
527 doeta1=exp(coef_b2-sqrt(coef_b3-coef_b1*log(dolam1)))
528 lamorb=0.62_r8*d0
529 lamanorb=535.0_r8*d50
530 IF (doeta1.lt.20.0_r8) THEN
531 dolam=1.0_r8/0.62_r8
532 ELSE IF (doeta1.gt.100.0_r8) THEN
533 dolam=dolam1
534 ELSE
535 fdo_etaano=-log(lamorb/lamanorb)* &
536 & log(0.01_r8*doeta1)/log(5.0_r8)
537 dolam=dolam1*exp(-fdo_etaano)
538 END IF
539 doeta2=exp(coef_b2-sqrt(coef_b3-coef_b1*log(dolam)))
540 rheight(i,j)=d0/doeta2
541 rlength(i,j)=d0/dolam
542 END IF
543!
544! Value of ar can range from 0.3 to 3 (Soulsby, 1997, p. 124)
545!
546 zobf(i,j)=ar*rheight(i,j)**2/rlength(i,j)
547 END IF
548 zo=zon(i,j)
549# ifdef SSW_ZOBL
550 zo=zo+zost(i,j)
551# endif
552# ifdef SSW_ZORIP
553 zo=zo+zobf(i,j)
554# endif
555# ifdef SSW_ZOBIO
556 zo=zo+zobio(i,j)
557# endif
558#else
559 IF (zodef(i,j).lt.absolute_zomin) THEN
560 zodef(i,j)=absolute_zomin
561 IF (master) THEN
562 WRITE (stdout,*) ' Warning: default zo < 0.05 mm,', &
563 & ' replaced with: ', zodef
564 END IF
565 END IF
566 zo=zodef(i,j)
567#endif
568!
569! Compute stresses.
570!
571! Default stress calcs for pure currents
572!
573 zo=min(max(zo,zomin),zomax)
574!
575 cff1=vonkar/log(zr(i,j)/zo)
576 cff2=min(cdb_max,max(cdb_min,cff1*cff1))
577 tauc(i,j)=cff2*umag(i,j)*umag(i,j)
578 tauw(i,j)=0.0_r8
579 taucwmax(i,j)=tauc(i,j)
580 znot(i,j)=zo
581 znotc(i,j)=zo
582#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
583 ksd_wbl(i,j)=zo
584 ustrc_wbl(i,j)=sqrt(tauc(i,j)+eps)
585#endif
586!
587 IF ((umag(i,j).le.eps).and.(ub(i,j).ge.eps)) THEN
588!
589! Pure waves - use wave friction factor approach from Madsen
590! (1994, eqns 32-33).
591!
592 sg_abokb=ab(i,j)/(30.0_r8*zo)
593 sg_fwm=0.3_r8
594 IF ((sg_abokb.gt.0.2_r8).and.(sg_abokb.le.100.0_r8)) THEN
595 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
596 ELSE IF (sg_abokb.gt.100.0_r8)THEN
597 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
598 END IF
599 tauc(i,j)= 0.0_r8
600 tauw(i,j)= 0.5_r8*sg_fwm*ub(i,j)*ub(i,j)
601 taucwmax(i,j)=tauw(i,j)
602 znot(i,j)=zo
603 znotc(i,j)=zo
604 ELSE IF ((umag(i,j).gt.eps).and.(ub(i,j).ge.eps).and. &
605 & ((zr(i,j)/zo).le.1.0_r8)) THEN
606!
607! Waves and currents, but zr <= zo.
608!
609 IF (master) THEN
610 WRITE (stdout,*) ' Warning: w-c calcs ignored because', &
611 & ' zr <= zo'
612 END IF
613 ELSE IF ((umag(i,j).gt.eps).and.(ub(i,j).ge.eps).and. &
614 & ((zr(i,j)/zo).gt.1.0_r8)) THEN
615!
616! Waves and currents, zr > zo.
617!
618#if defined SGWC
619 sg_zrozn=zr(i,j)/zo
620 sg_ubokur=ub(i,j)/(sg_kappa*umag(i,j))
621 sg_row=ab(i,j)/zo
622 sg_a1=1.0d-6
623 sg_phicw=phicw(i,j)
624 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
625 & sg_a1, sg_mu, sg_epsilon, sg_ro, sg_fofa)
626 sg_abokb=ab(i,j)/(30.0_r8*zo)
627 IF (sg_abokb.le.100.0_r8) THEN
628 sg_fwm=exp(-8.82_r8+7.02_r8*sg_abokb**(-0.078_r8))
629 ELSE
630 sg_fwm=exp(-7.30_r8+5.61_r8*sg_abokb**(-0.109_r8))
631 END IF
632 sg_ubouwm=sqrt(2.0_r8/sg_fwm)
633!
634! Determine the maximum ratio of wave over combined shear stresses,
635! sg_ubouwm (ub/ustarwm).
636!
637 CALL sg_purewave (sg_row, sg_ubouwm, sg_znotp, sg_ro)
638!
639! Set initial guess of the ratio of wave over shear stress, sg_c1
640! (ub/ustarc).
641!
642 sg_b1=sg_ubouwm
643 sg_fofb=-sg_fofa
644 sg_c1=0.5_r8*(sg_a1+sg_b1)
645 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
646 & sg_c1, sg_mu, sg_epsilon, sg_ro, sg_fofc)
647!
648! Solve PDE via bi-section method.
649!
650 iterate=.true.
651 DO iter=1,sg_n
652 IF (iterate) THEN
653 IF ((sg_fofb*sg_fofc).lt.0.0_r8) THEN
654 sg_a1=sg_c1
655 ELSE
656 sg_b1=sg_c1
657 END IF
658 sg_c1=0.5_r8*(sg_a1+sg_b1)
659 CALL sg_bstress (sg_row, sg_zrozn, sg_phicw, sg_ubokur, &
660 & sg_c1, sg_mu, sg_epsilon, sg_ro, &
661 & sg_fofc)
662 iterate=(sg_b1-sg_c1) .ge. sg_tol
663 IF (iterate) iconv(i,j)=iter
664 END IF
665 END DO
666 sg_ubouc=sg_c1
667!
668! Compute bottom shear stress magnitude (m/s).
669!
670 sg_ustarcw=ub(i,j)/sg_ubouc
671 sg_ustarwm=sg_mu*sg_ustarcw
672!! sg_ustarc=MIN(sg_ustarcdef,sg_epsilon*sg_ustarcw) !original
673 sg_ustarc=max(sqrt(tauc(i,j)),sg_epsilon*sg_ustarcw)
674 tauc(i,j)=sg_ustarc*sg_ustarc
675 tauw(i,j)=sg_ustarwm*sg_ustarwm
676 taucwmax(i,j)=sqrt((tauc(i,j)+ &
677 & tauw(i,j)*cos(phicw(i,j)))**2+ &
678 & (tauw(i,j)*sin(phicw(i,j)))**2)
679!
680! Compute apparent hydraulic roughness (m).
681!
682 IF (sg_epsilon.gt.0.0_r8) THEN
683 sg_z1=sg_alpha*sg_kappa*ab(i,j)/sg_ubouc
684 sg_z2=sg_z1/sg_epsilon
685 sg_z1ozn=sg_z1/zo
686 znotc(i,j)=sg_z2* &
687 & exp(-(1.0_r8-sg_epsilon+ &
688 & sg_epsilon*log(sg_z1ozn)))
689!
690! Compute mean (m/s) current at 100 cm above the bottom.
691!
692 IF (sg_z100.gt.sg_z2) THEN
693 u100(i,j)=sg_ustarc* &
694 & (log(sg_z100/sg_z2)+1.0_r8-sg_epsilon+ &
695 & sg_epsilon*log(sg_z1ozn))/sg_kappa
696 ELSE IF ((sg_z100.le.sg_z2).and.(zr(i,j).gt.sg_z1)) THEN
697 u100(i,j)=sg_ustarc*sg_epsilon* &
698 & (sg_z100/sg_z1-1.0_r8+log(sg_z1ozn))/sg_kappa
699 ELSE
700 u100(i,j)=sg_ustarc*sg_epsilon* &
701 & log(sg_z100/zo)/sg_kappa
702 END IF
703 END IF
704#elif defined M94WC
705 m_ubr=ub(i,j)
706 m_wr=fwave_bot(i,j)
707 m_ucr=umag(i,j)
708 m_zr=zr(i,j)
709 m_phicw=phicw(i,j)
710 m_kb=30.0_r8*zo
711 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
712 CALL madsen94 (m_ubr, m_wr, m_ucr, &
713 & m_zr, m_phicw, m_kb, dstp, &
714 & m_ustrc, m_ustrwm, m_ustrr, m_fwc, m_zoa, &
715 & m_dwc)
716 tauc(i,j)=m_ustrc*m_ustrc
717 tauw(i,j)=m_ustrwm*m_ustrwm
718 taucwmax(i,j)=m_ustrr*m_ustrr
719 znotc(i,j)=min( m_zoa, zomax )
720 u100(i,j)=(m_ustrc/vonkar)*log(1.0_r8/m_zoa)
721 thck_wbl(i,j)=m_dwc
722#endif
723#if defined SSW_FORM_DRAG_COR
724 IF (rheight(i,j).gt.(zon(i,j)+zost(i,j))) THEN
725 coef_fd=0.5_r8*cd_fd*(rheight(i,j)/rlength(i,j))* &
726 & (1.0_r8/(vonkar*vonkar))* &
727 & (log(rheight(i,j)/ &
728 & (zon(i,j)+zost(i,j)))-1.0_r8)**2
729 taucwmax(i,j)=taucwmax(i,j)/(1.0_r8+coef_fd)
730 taucwmax(i,j)=taucwmax(i,j)*(1.0_r8+8.0_r8* &
731 & rheight(i,j)/rlength(i,j))
732 END IF
733#endif
734#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
735 ksd_wbl(i,j)=m_zoa
736 ustrc_wbl(i,j)=m_ustrc
737#endif
738 END IF
739 END DO
740 END DO
741!
742#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
743!
744! Find the near-bottom current velocity(udelta) at a given elevation.
745! Use the Madsen output of current shear stress, apparent roughness
746! to get udelta.
747! Find the angle at that near-bottom current velocity.
748!
749 DO j=jstr,jend
750 DO i=istr,iend
751!
752 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
753!
754! Use wave boundary layer (wbl) thickness based on Madsen to get
755! near bottom current velocity.
756!
757 cff=min( 0.98_r8*dstp, thck_wbl(i,j) )
758!
759! Make sure that wbl is under total depth and greater than
760! apparent roughness.
761!
762 cff1=max(cff, 1.1_r8*ksd_wbl(i,j))
763 cff2=log(cff1/ksd_wbl(i,j))
764 udelta_wbl(i,j)=(ustrc_wbl(i,j)/vonkar)*cff2
765!
766 DO k=1,n(ng)
767 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
768 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
769# ifdef SSW_LOGINT_STOKES
770 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
771 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
772# endif
773 END DO
774 CALL log_interp( n(ng), dstp, cff1, &
775 & urz, vrz, &
776 & z_r(i,j,:), z_w(i,j,:), &
777 & bottom(i,j,isd50), bottom(i,j,izapp), &
778 & zr_wbl(i,j), &
779 & ur_sgwbl(i,j), vr_sgwbl(i,j) )
780!
781! Compute bottom current magnitude at RHO-points.
782!
783 ucur_sgwbl=ur_sgwbl(i,j)
784 vcur_sgwbl=vr_sgwbl(i,j)
785!
786 IF (ucur_sgwbl.eq.0.0_r8) THEN
787 phic_sgwbl(i,j)=0.5_r8*pi*sign(1.0_r8,vcur_sgwbl)
788 ELSE
789 phic_sgwbl(i,j)=atan2(vcur_sgwbl,ucur_sgwbl)
790 ENDIF
791 END DO
792 END DO
793#endif
794#if defined BEDLOAD_VANDERA_DIRECT_UDELTA
795!
796! Find the near-bottom current velocity directly at a given
797! elevation (doesnot require Madsen output)
798! Find the angle at that near-bottom current velocity.
799!
800 DO j=jstr,jend
801 DO i=istr,iend
802 dstp=z_r(i,j,n(ng))-z_w(i,j,0)
803 cff=min( 0.98_r8*dstp, sg_zwbl(ng) )
804 DO k=1,n(ng)
805 urz(k)=0.5_r8*(u(i,j,k,nrhs)+u(i+1,j,k,nrhs))
806 vrz(k)=0.5_r8*(v(i,j,k,nrhs)+v(i,j+1,k,nrhs))
807# ifdef SSW_LOGINT_STOKES
808 urz(k)=urz(k)+0.5_r8*(u_stokes(i,j,k)+u_stokes(i+1,j,k))
809 vrz(k)=vrz(k)+0.5_r8*(v_stokes(i,j,k)+v_stokes(i,j+1,k))
810# endif
811 END DO
812 CALL log_interp( n(ng), dstp, cff, &
813 & urz, vrz, &
814 & z_r(i,j,:), z_w(i,j,:), &
815 & bottom(i,j,isd50), bottom(i,j,izapp), &
816 & zr_wbl(i,j), &
817 & ur_sgwbl(i,j), vr_sgwbl(i,j) )
818!
819! Compute bottom current magnitude at RHO-points.
820!
821 ucur_sgwbl=ur_sgwbl(i,j)
822 vcur_sgwbl=vr_sgwbl(i,j)
823 udelta_wbl(i,j)=sqrt(ur_sgwbl(i,j)*ur_sgwbl(i,j)+ &
824 & vr_sgwbl(i,j)*vr_sgwbl(i,j)+eps)
825!
826 IF (ucur_sgwbl.eq.0.0_r8) THEN
827 phic_sgwbl(i,j)=0.5_r8*pi*sign(1.0_r8,vcur_sgwbl)
828 ELSE
829 phic_sgwbl(i,j)=atan2(vcur_sgwbl,ucur_sgwbl)
830 ENDIF
831 END DO
832 END DO
833!
834#endif
835!
836!-----------------------------------------------------------------------
837! Compute kinematic bottom stress components due current and wind-
838! induced waves.
839!-----------------------------------------------------------------------
840!
841 DO j=jstr,jend
842 DO i=istru,iend
843 anglec=0.5_r8*(ur_sg(i,j)+ur_sg(i-1,j))/ &
844 & (0.5_r8*(umag(i-1,j)+umag(i,j)))
845 bustr(i,j)=0.5_r8*(tauc(i-1,j)+tauc(i,j))*anglec
846#ifdef WET_DRY
847 cff2=0.75_r8*0.5_r8*(z_w(i-1,j,1)+z_w(i,j,1)- &
848 & z_w(i-1,j,0)-z_w(i,j,0))
849 bustr(i,j)=sign(1.0_r8,bustr(i,j))*min(abs(bustr(i,j)), &
850 & abs(u(i,j,1,nrhs))*cff2/dt(ng))
851#endif
852 END DO
853 END DO
854 DO j=jstrv,jend
855 DO i=istr,iend
856 anglec=0.5_r8*(vr_sg(i,j)+vr_sg(i,j-1))/ &
857 & (0.5_r8*(umag(i,j-1)+umag(i,j)))
858 bvstr(i,j)=0.5_r8*(tauc(i,j-1)+tauc(i,j))*anglec
859#ifdef WET_DRY
860 cff2=0.75_r8*0.5_r8*(z_w(i,j-1,1)+z_w(i,j,1)- &
861 & z_w(i,j-1,0)-z_w(i,j,0))
862 bvstr(i,j)=sign(1.0_r8,bvstr(i,j))*min(abs(bvstr(i,j)), &
863 & abs(v(i,j,1,nrhs))*cff2/dt(ng))
864#endif
865 END DO
866 END DO
867 DO j=jstr,jend
868 DO i=istr,iend
869 anglec=ucur(i,j)/umag(i,j)
870 anglew=cos(1.5_r8*pi-dwave(i,j)-angler(i,j))
871 bustrc(i,j)=tauc(i,j)*anglec
872 bustrw(i,j)=tauw(i,j)*anglew
873 bustrcwmax(i,j)=taucwmax(i,j)*anglew
874 ubot(i,j)=ub(i,j)*anglew
875 ur(i,j)=ucur(i,j)
876!
877 anglec=vcur(i,j)/umag(i,j)
878 anglew=sin(1.5_r8*pi-dwave(i,j)-angler(i,j))
879 bvstrc(i,j)=tauc(i,j)*anglec
880 bvstrw(i,j)=tauw(i,j)*anglew
881 bvstrcwmax(i,j)=taucwmax(i,j)*anglew
882 vbot(i,j)=ub(i,j)*anglew
883 vr(i,j)=vcur(i,j)
884!
885 bottom(i,j,irlen)=rlength(i,j)
886 bottom(i,j,irhgt)=rheight(i,j)
887 bottom(i,j,ibwav)=ab(i,j)
888 bottom(i,j,izdef)=zodef(i,j)
889 bottom(i,j,izapp)=znotc(i,j)
890 bottom(i,j,iznik)=zon(i,j)
891 bottom(i,j,izbio)=zobio(i,j)
892 bottom(i,j,izbfm)=zobf(i,j)
893 bottom(i,j,izbld)=zost(i,j)
894 bottom(i,j,izwbl)=znot(i,j)
895 bottom(i,j,idtbl)=thck_wbl(i,j)
896#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
897 bottom(i,j,idksd)=ksd_wbl(i,j)
898 bottom(i,j,idusc)=ustrc_wbl(i,j)
899#endif
900#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
901 defined bedload_vandera_direct_udelta
902 bottom(i,j,idubl)=udelta_wbl(i,j)
903 bottom(i,j,idzrw)=zr_wbl(i,j)
904 bottom(i,j,idpcx)=phic_sgwbl(i,j)
905#else
906 bottom(i,j,idpcx)=phic(i,j)
907 bottom(i,j,idpwc)=phicw(i,j)
908#endif
909 END DO
910 END DO
911!
912! Apply periodic or gradient boundary conditions for output
913! purposes only.
914!
915 CALL bc_u2d_tile (ng, tile, &
916 & lbi, ubi, lbj, ubj, &
917 & bustr)
918 CALL bc_v2d_tile (ng, tile, &
919 & lbi, ubi, lbj, ubj, &
920 & bvstr)
921 CALL bc_r2d_tile (ng, tile, &
922 & lbi, ubi, lbj, ubj, &
923 & bustrc)
924 CALL bc_r2d_tile (ng, tile, &
925 & lbi, ubi, lbj, ubj, &
926 & bvstrc)
927 CALL bc_r2d_tile (ng, tile, &
928 & lbi, ubi, lbj, ubj, &
929 & bustrw)
930 CALL bc_r2d_tile (ng, tile, &
931 & lbi, ubi, lbj, ubj, &
932 & bvstrw)
933 CALL bc_r2d_tile (ng, tile, &
934 & lbi, ubi, lbj, ubj, &
935 & bustrcwmax)
936 CALL bc_r2d_tile (ng, tile, &
937 & lbi, ubi, lbj, ubj, &
938 & bvstrcwmax)
939 CALL bc_r2d_tile (ng, tile, &
940 & lbi, ubi, lbj, ubj, &
941 & ubot)
942 CALL bc_r2d_tile (ng, tile, &
943 & lbi, ubi, lbj, ubj, &
944 & vbot)
945 CALL bc_r2d_tile (ng, tile, &
946 & lbi, ubi, lbj, ubj, &
947 & ur)
948 CALL bc_r2d_tile (ng, tile, &
949 & lbi, ubi, lbj, ubj, &
950 & vr)
951 CALL bc_r2d_tile (ng, tile, &
952 & lbi, ubi, lbj, ubj, &
953 & bottom(:,:,irlen))
954 CALL bc_r2d_tile (ng, tile, &
955 & lbi, ubi, lbj, ubj, &
956 & bottom(:,:,irhgt))
957 CALL bc_r2d_tile (ng, tile, &
958 & lbi, ubi, lbj, ubj, &
959 & bottom(:,:,ibwav))
960 CALL bc_r2d_tile (ng, tile, &
961 & lbi, ubi, lbj, ubj, &
962 & bottom(:,:,izdef))
963 CALL bc_r2d_tile (ng, tile, &
964 & lbi, ubi, lbj, ubj, &
965 & bottom(:,:,izapp))
966 CALL bc_r2d_tile (ng, tile, &
967 & lbi, ubi, lbj, ubj, &
968 & bottom(:,:,iznik))
969 CALL bc_r2d_tile (ng, tile, &
970 & lbi, ubi, lbj, ubj, &
971 & bottom(:,:,izbio))
972 CALL bc_r2d_tile (ng, tile, &
973 & lbi, ubi, lbj, ubj, &
974 & bottom(:,:,izbfm))
975 CALL bc_r2d_tile (ng, tile, &
976 & lbi, ubi, lbj, ubj, &
977 & bottom(:,:,izbld))
978 CALL bc_r2d_tile (ng, tile, &
979 & lbi, ubi, lbj, ubj, &
980 & bottom(:,:,izwbl))
981#if defined BEDLOAD_VANDERA_MADSEN_UDELTA
982 CALL bc_r2d_tile (ng, tile, &
983 & lbi, ubi, lbj, ubj, &
984 & bottom(:,:,idtbl))
985 CALL bc_r2d_tile (ng, tile, &
986 & lbi, ubi, lbj, ubj, &
987 & bottom(:,:,idksd))
988 CALL bc_r2d_tile (ng, tile, &
989 & lbi, ubi, lbj, ubj, &
990 & bottom(:,:,idusc))
991#endif
992#if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
993 defined bedload_vandera_direct_udelta
994 CALL bc_r2d_tile (ng, tile, &
995 & lbi, ubi, lbj, ubj, &
996 & bottom(:,:,idubl))
997 CALL bc_r2d_tile (ng, tile, &
998 & lbi, ubi, lbj, ubj, &
999 & bottom(:,:,idzrw))
1000#else
1001 CALL bc_r2d_tile (ng, tile, &
1002 & lbi, ubi, lbj, ubj, &
1003 & bottom(:,:,idpwc))
1004#endif
1005 CALL bc_r2d_tile (ng, tile, &
1006 & lbi, ubi, lbj, ubj, &
1007 & bottom(:,:,idpcx))
1008#ifdef DISTRIBUTE
1009 CALL mp_exchange2d (ng, tile, inlm, 4, &
1010 & lbi, ubi, lbj, ubj, &
1011 & nghostpoints, &
1012 & ewperiodic(ng), nsperiodic(ng), &
1013 & bustr, bvstr, bustrc, bvstrc)
1014 CALL mp_exchange2d (ng, tile, inlm, 4, &
1015 & lbi, ubi, lbj, ubj, &
1016 & nghostpoints, &
1017 & ewperiodic(ng), nsperiodic(ng), &
1018 & bustrw, bvstrw, bustrcwmax, bvstrcwmax)
1019 CALL mp_exchange2d (ng, tile, inlm, 4, &
1020 & lbi, ubi, lbj, ubj, &
1021 & nghostpoints, &
1022 & ewperiodic(ng), nsperiodic(ng), &
1023 & ubot, vbot, ur, vr)
1024 CALL mp_exchange2d (ng, tile, inlm, 3, &
1025 & lbi, ubi, lbj, ubj, &
1026 & nghostpoints, &
1027 & ewperiodic(ng), nsperiodic(ng), &
1028 & bottom(:,:,irlen), &
1029 & bottom(:,:,irhgt), &
1030 & bottom(:,:,ibwav))
1031 CALL mp_exchange2d (ng, tile, inlm, 3, &
1032 & lbi, ubi, lbj, ubj, &
1033 & nghostpoints, &
1034 & ewperiodic(ng), nsperiodic(ng), &
1035 & bottom(:,:,izdef), &
1036 & bottom(:,:,izapp), &
1037 & bottom(:,:,iznik))
1038 CALL mp_exchange2d (ng, tile, inlm, 4, &
1039 & lbi, ubi, lbj, ubj, &
1040 & nghostpoints, &
1041 & ewperiodic(ng), nsperiodic(ng), &
1042 & bottom(:,:,izbio), &
1043 & bottom(:,:,izbfm), &
1044 & bottom(:,:,izbld), &
1045 & bottom(:,:,izwbl))
1046# if defined BEDLOAD_VANDERA_MADSEN_UDELTA
1047 CALL mp_exchange2d (ng, tile, inlm, 3, &
1048 & lbi, ubi, lbj, ubj, &
1049 & nghostpoints, &
1050 & ewperiodic(ng), nsperiodic(ng), &
1051 & bottom(:,:,idtbl), &
1052 & bottom(:,:,idksd), &
1053 & bottom(:,:,idusc))
1054# endif
1055# if defined BEDLOAD_VANDERA_MADSEN_UDELTA || \
1056 defined bedload_vandera_direct_udelta
1057 CALL mp_exchange2d (ng, tile, inlm, 2, &
1058 & lbi, ubi, lbj, ubj, &
1059 & nghostpoints, &
1060 & ewperiodic(ng), nsperiodic(ng), &
1061 & bottom(:,:,idzrw), &
1062 & bottom(:,:,idubl))
1063# else
1064 CALL mp_exchange2d (ng, tile, inlm, 1, &
1065 & lbi, ubi, lbj, ubj, &
1066 & nghostpoints, &
1067 & ewperiodic(ng), nsperiodic(ng), &
1068 & bottom(:,:,idpwc))
1069# endif
1070 CALL mp_exchange2d (ng, tile, inlm, 1, &
1071 & lbi, ubi, lbj, ubj, &
1072 & nghostpoints, &
1073 & ewperiodic(ng), nsperiodic(ng), &
1074 & bottom(:,:,idpcx))
1075#endif
1076!
1077 RETURN
mod_iounits
Definition mod_iounits.F:2
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_sediment::iznik
integer, parameter iznik
Definition sediment_mod.h:133
mod_sediment::izbld
integer, parameter izbld
Definition sediment_mod.h:136
mod_sediment::izbio
integer, parameter izbio
Definition sediment_mod.h:134
mod_sediment::izwbl
integer, parameter izwbl
Definition sediment_mod.h:137
mod_sediment::izbfm
integer, parameter izbfm
Definition sediment_mod.h:135

References bc_2d_mod::bc_r2d_tile(), bc_2d_mod::bc_u2d_tile(), bc_2d_mod::bc_v2d_tile(), mod_scalars::cdb_max, mod_scalars::cdb_min, mod_scalars::dt, mod_scalars::ewperiodic, mod_scalars::g, mod_sediment::ibwav, mod_sediment::idens, mod_param::inlm, mod_sediment::irhgt, mod_sediment::irlen, mod_sediment::isd50, mod_sediment::itauc, mod_sediment::iwsed, mod_sediment::izapp, mod_sediment::izbfm, mod_sediment::izbio, mod_sediment::izbld, mod_sediment::izdef, mod_sediment::iznik, mod_sediment::izwbl, log_interp(), madsen94(), mod_parallel::master, mp_exchange_mod::mp_exchange2d(), mod_param::nghostpoints, mod_scalars::nsperiodic, mod_scalars::pi, mod_scalars::sg_alpha, sg_bstress(), mod_scalars::sg_kappa, sg_kelvin8m(), sg_kelvin8p(), mod_scalars::sg_mp, mod_scalars::sg_n, sg_purewave(), mod_scalars::sg_tol, mod_scalars::sg_z100, mod_scalars::sg_z1min, mod_scalars::sg_z1p, mod_iounits::stdout, and mod_scalars::vonkar.

Here is the call graph for this function: