Skip to content

Commit

Permalink
Add notes for WWM vortex formulation.
Browse files Browse the repository at this point in the history
  • Loading branch information
@josephzhang8
josephzhang8 committed Mar 19, 2024
1 parent ef4a584 commit d325ed3
Show file tree
Hide file tree
Showing 3 changed files with 66 additions and 53 deletions.
6 changes: 3 additions & 3 deletions src/Utility/Cluster_files/run_kuro
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
#!/bin/tcsh
#SBATCH --job-name=R01a
##64 cores/nodes
##64 cores/nodes. Use sinfo to see system status
#SBATCH -N 13 --ntasks-per-node=64
##Max wall clock time is XX days
##Max wall clock time is 2 days
#SBATCH -t 24:00:00

srun ./pschism_i24 6
srun ./pschism_KURO_BLD_STANDALONE_TVD-VL 6
111 changes: 62 additions & 49 deletions src/WWMIII/wwm_coupl_selfe.F90
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
!Note: most arrays in this file are from SCHISM directly (too account for
!Note: most arrays in this file are from SCHISM directly (to account for
!quads)
#include "wwm_functions.h"
#ifdef SCHISM
Expand Down Expand Up @@ -87,7 +87,7 @@ SUBROUTINE RADIATION_STRESS_SCHISM
! Total water depth at sides
HTOT = MAX((DEP(isidenode(1,IS)) + DEP(isidenode(2,IS)))/2.0D0,hmin_radstress)

! Wave forces
! Wave forces (HTOT>0)
WWAVE_FORCE(1,:,IS) = WWAVE_FORCE(1,:,IS) - (DSXX3D(1,:,IS) + DSXY3D(2,:,IS)) / HTOT
WWAVE_FORCE(2,:,IS) = WWAVE_FORCE(2,:,IS) - (DSXY3D(1,:,IS) + DSYY3D(2,:,IS)) / HTOT
END DO !IS
Expand Down Expand Up @@ -123,7 +123,7 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
IF(idry(IP) == 1) CYCLE

! Total water depth at the node
D_loc = MAX( DEP(IP) , hmin_radstress )
D_loc = MAX( DEP(IP) , hmin_radstress ) !>0

! Initialization of the local Stokes drift and J variables
USTOKES_loc = 0.D0; VSTOKES_loc = 0.D0;
Expand All @@ -135,6 +135,10 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
Uint = 0.D0; Vint = 0.D0; Urint = 0.D0; Vrint = 0.D0
k_loc = MIN(KDMAX/DEP(IP),WK(IS,IP))
kD_loc = MIN(KDMAX,WK(IS,IP)*D_loc)
IF(kD_loc == 0) THEN
WRITE(errmsg,*)'WWM: kD_loc=0'
CALL parallel_abort(errmsg)
END IF

! Loop on the directions
DO ID = 1, MDC
Expand Down Expand Up @@ -162,6 +166,10 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
! influence Wst, which is computed from the continuity equation for waves only

IF (IROLLER == 1) THEN
IF(CROLP(IP)== 0) THEN
WRITE(errmsg,*)'WWM: CROLP(IP)=0'
CALL parallel_abort(errmsg)
END IF
Urint = 2.D0*COS(DROLP(IP))*EROL2(IP)/(CROLP(IP)*D_loc)
Vrint = 2.D0*SIN(DROLP(IP))*EROL2(IP)/(CROLP(IP)*D_loc)

Expand Down Expand Up @@ -217,7 +225,8 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
END DO
END DO !nsa

!... Compute bottom Stokes drift z-vel. at elements
!... Compute bottom Stokes drift w-vel. at elements
! Used only first 3 nodes of quad
STOKES_WVEL_ELEM = 0.D0
DO IE = 1,nea
IF(idry_e(IE) == 1) CYCLE
Expand All @@ -240,41 +249,41 @@ SUBROUTINE STOKES_STRESS_INTEGRAL_SCHISM
STOKES_WVEL_ELEM(kbe(IE),IE) = -dhdx*ubar - dhdy*vbar
END DO !nea

!... Compute bottom Stokes z-vel. at nodes
!... Compute bottom Stokes w-vel. at nodes
STOKES_WVEL = 0.D0
DO IP = 1,np !residents only
IF(idry(IP) == 1) CYCLE

!Bottom Stokes z-vel.
!Bottom Stokes w-vel.
tmp0 = 0.D0
DO j = 1,nne(IP)
ie = indel(j,IP)
IF(idry_e(ie)==0) THEN
STOKES_WVEL(kbp(IP),IP) = STOKES_WVEL(kbp(IP),IP) + STOKES_WVEL_ELEM(kbe(ie),ie)*area(ie)
END IF
tmp0 = tmp0 + area(ie)
tmp0 = tmp0 + area(ie) !>0
END DO !j
STOKES_WVEL(kbp(IP),IP) = STOKES_WVEL(kbp(IP),IP)/tmp0
END DO !np

!... Compute horizontal gradient of Stokes x and y-vel. (to compute Stokes z-vel.)
!... Compute horizontal gradient of Stokes x and y-vel. (to compute Stokes w-vel.)
ws_tmp1 = 0.D0; ws_tmp2 = 0.D0
CALL hgrad_nodes(2,0,NVRT,MNP,nsa,STOKES_HVEL(1,:,:),dr_dxy_loc)
ws_tmp1(:,:) = dr_dxy_loc(1,:,:) !valid only in resident
CALL hgrad_nodes(2,0,NVRT,MNP,nsa,STOKES_HVEL(2,:,:),dr_dxy_loc)
ws_tmp2(:,:) = dr_dxy_loc(2,:,:)

!... Compute Stokes z-vel. at all levels: STOKES_WVEL_SIDE(NVRT,nsa)
!... Compute Stokes w-vel. at all levels: STOKES_WVEL_SIDE(NVRT,nsa)
STOKES_WVEL_SIDE = 0.D0
DO IS = 1,ns !residents only
IF(idry_s(IS) == 1) CYCLE
n1 = isidenode(1,IS)
n2 = isidenode(2,IS)

!Bottom Stokes z-vel.
!Bottom Stokes w-vel.
STOKES_WVEL_SIDE(kbs(IS),IS) = (STOKES_WVEL(max(kbs(IS),kbp(n1)),n1) + STOKES_WVEL(max(kbs(IS),kbp(n2)),n2))/2.D0

!Stokes z-vel. at all levels
!Stokes w-vel. at all levels
DO k = kbs(IS)+1, NVRT
ztmp = zs(k,IS) - zs(k-1,IS)
STOKES_WVEL_SIDE(k,IS) = STOKES_WVEL_SIDE(k-1,IS) &
Expand Down Expand Up @@ -404,14 +413,14 @@ END SUBROUTINE COMPUTE_CONSERVATIVE_VF_TERMS_SCHISM

!**********************************************************************
!* This routine is used with RADFLAG=VOR (3D vortex formulation, after Bennis et al., 2011)
!* => Computation of the non-conservative terms due to depth-induced breaking (term Fb from Eq. (11) and (12))
!* => Computation of the non-conservative terms due to depth-induced breaking (term Fd from Eq. (11) and (12))
!* March 2022 : update LRU team
! Accounts for depth-induced breaking, roller (if turned on) and whitecapping contribution
!**********************************************************************
SUBROUTINE COMPUTE_BREAKING_VF_TERMS_SCHISM
USE DATAPOOL
USE schism_glbl, ONLY: hmin_radstress, kbs, ns, isbs, dps, h0, out_wwm, &
&zs,nsa,idry_s,isidenode
&zs,nsa,idry_s,isidenode,errmsg
USE schism_msgp
IMPLICIT NONE

Expand All @@ -430,7 +439,7 @@ SUBROUTINE COMPUTE_BREAKING_VF_TERMS_SCHISM
n1 = isidenode(1,IS); n2 = isidenode(2,IS)
eta_tmp = (eta2(n1) + eta2(n2))/2.D0
!htot = max(h0,dps(IS)+eta_tmp,hmin_radstress) ! KM
htot = max(h0,dps(IS)+eta_tmp)
htot = max(h0,dps(IS)+eta_tmp) !>0

IF(kbs(IS)+1 == NVRT) THEN !2D

Expand Down Expand Up @@ -469,9 +478,9 @@ SUBROUTINE COMPUTE_BREAKING_VF_TERMS_SCHISM
IF (ZPROF_BREAK == 1) swild_3D(k) = 1.D0
! Hyperbolic distribution
IF (ZPROF_BREAK == 2) swild_3D(k) = cosh((dps(IS)+zs(k,IS))/(0.2D0*tmp0))
IF (ZPROF_BREAK == 3) swild_3D(k) = 1.D0 - dtanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**2.D0)
IF (ZPROF_BREAK == 4) swild_3D(k) = 1.D0 - dtanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**4.D0)
IF (ZPROF_BREAK == 5) swild_3D(k) = 1.D0 - dtanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**8.D0)
IF (ZPROF_BREAK == 3) swild_3D(k) = 1.D0 - tanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**2.D0)
IF (ZPROF_BREAK == 4) swild_3D(k) = 1.D0 - tanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**4.D0)
IF (ZPROF_BREAK == 5) swild_3D(k) = 1.D0 - tanh(((eta_tmp-zs(k,IS))/(0.5D0*tmp0))**8.D0)
! All in the two surface layers
IF (ZPROF_BREAK == 6 .AND. k .GE. NVRT-1) swild_3D(k)=1.D0
END DO !k
Expand All @@ -486,6 +495,10 @@ SUBROUTINE COMPUTE_BREAKING_VF_TERMS_SCHISM
DO k = kbs(IS), NVRT-1
sum_3D = sum_3D + (swild_3D(k+1) + swild_3D(k))/2.D0*(zs(k+1,IS) - zs(k,IS))
END DO !NVRT-1
IF(sum_3D==0) THEN
WRITE(errmsg,*)'WWM: sum_3D=0'
CALL parallel_abort(errmsg)
END IF

DO k = kbs(IS), NVRT

Expand Down Expand Up @@ -524,14 +537,14 @@ END SUBROUTINE COMPUTE_BREAKING_VF_TERMS_SCHISM

!**********************************************************************
!* This routine is used with RADFLAG=VOR (3D vortex formulation, after Bennis et al., 2011)
!* => Computation of the non-conservative terms due to bottom friction (see Uchiyama et al., 2010)
!* => Computation of the non-conservative terms (Fb) due to bottom friction (see Uchiyama et al., 2010)
!* TO DO : pass the vertical distribution in option similar to breaking wave force
!**********************************************************************
SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
! MP
USE DATAPOOL
USE schism_glbl, ONLY: hmin_radstress, kbs, ns, isbs, dps, h0, out_wwm,&
&zs, idry_s, isidenode, nchi, rough_p, iwbl, delta_wbl
&zs,idry_s,isidenode,nchi,rough_p,iwbl,delta_wbl,errmsg
USE schism_msgp
IMPLICIT NONE

Expand All @@ -549,7 +562,7 @@ SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
n1 = isidenode(1,IS); n2 = isidenode(2,IS)
eta_tmp = (eta2(n1) + eta2(n2))/2.D0
!htot = max(h0,dps(IS)+eta_tmp,hmin_radstress) ! KM
htot = max(h0,dps(IS)+eta_tmp)
htot = max(h0,dps(IS)+eta_tmp) !>0

IF(kbs(IS)+1 == NVRT) THEN !2D

Expand All @@ -568,6 +581,7 @@ SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
! we note 1/kwd = tmp0
tmp0 = (delta_wbl(n1) + delta_wbl(n2))/2.D0
IF(tmp0 .LT. SMALL) CYCLE
!tmp0>0

! Vertical distribution function of qdm
swild_3D = 0.D0
Expand All @@ -577,9 +591,9 @@ SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
! Hyperbolic distribution - Type of profile 1
!swild_3D(k) = cosh((eta_tmp-zs(k,IS))/tmp0)
! Hyperbolic distribution - Type of profile 2
swild_3D(k) = 1.D0 - dtanh(((dps(IS)+zs(k,IS))/tmp0)**2.D0)
!swild_3D(k) = 1.D0 - dtanh(((dps(IS)+zs(k,IS))/tmp0)**4.D0)
!swild_3D(k) = 1.D0 - dtanh(((dps(IS)+zs(k,IS))/tmp0)**8.D0)
swild_3D(k) = 1.D0 - tanh(((dps(IS)+zs(k,IS))/tmp0)**2.D0)
!swild_3D(k) = 1.D0 - tanh(((dps(IS)+zs(k,IS))/tmp0)**4.D0)
!swild_3D(k) = 1.D0 - tanh(((dps(IS)+zs(k,IS))/tmp0)**8.D0)
END DO !k

! In shallow depths, we make the vertical profile tend to a vertical-uniform one
Expand All @@ -591,13 +605,17 @@ SUBROUTINE COMPUTE_STREAMING_VF_TERMS_SCHISM
DO k = kbs(IS), NVRT-1
sum_3D = sum_3D + (swild_3D(k+1) + swild_3D(k))/2.D0*(zs(k+1,IS) - zs(k,IS))
END DO !NVRT-1
IF(sum_3D==0) THEN
WRITE(errmsg,*)'WWM: sum_3D=0'
CALL parallel_abort(errmsg)
END IF

DO k = kbs(IS), NVRT
Fws_x_loc = - swild_3D(k)*(SBF(1,n1) + SBF(1,n2))/2.D0/sum_3D
Fws_y_loc = - swild_3D(k)*(SBF(2,n1) + SBF(2,n2))/2.D0/sum_3D
! Saving wave streaming
WWAVE_FORCE(1,k,IS) = WWAVE_FORCE(1,k,IS) + Fws_x_loc
WWAVE_FORCE(2,k,IS) = WWAVE_FORCE(2,k,IS) + Fws_y_loc
WWAVE_FORCE(1,k,IS) = WWAVE_FORCE(1,k,IS) + Fws_x_loc
WWAVE_FORCE(2,k,IS) = WWAVE_FORCE(2,k,IS) + Fws_y_loc
END DO
END IF !2D/3D

Expand Down Expand Up @@ -639,8 +657,8 @@ SUBROUTINE COMPUTE_VEGDISS_VF_TERMS_SCHISM
n1 = isidenode(1,IS); n2 = isidenode(2,IS)
eta_tmp = (eta2(n1) + eta2(n2))/2.D0
!htot = max(h0,dps(IS)+eta_tmp,hmin_radstress)
htot = max(h0,dps(IS)+eta_tmp)
htot = max(h0,dps(IS)+eta_tmp) !>0

!vegetation height at sides
SAV_H_tmp = (VLTH(n1)+VLTH(n2))/2.D0

Expand All @@ -653,19 +671,18 @@ SUBROUTINE COMPUTE_VEGDISS_VF_TERMS_SCHISM
WWAVE_FORCE(2,:,IS) = WWAVE_FORCE(2,:,IS) + Fveg_y_loc

ELSE !3D

tmp0 = 0.D0
topveg_k = NVRT
IF (SAV_H_tmp >= htot) THEN !emergent vegetation
topveg_k = NVRT
ELSE
ELSE !submerged vegetation
DO k = kbs(IS), NVRT-1
tmp0 = tmp0 + zs(k+1,IS) - zs(k,IS)
IF (tmp0 >= SAV_H_tmp) THEN !submerged vegetation
IF (tmp0 >= SAV_H_tmp) THEN
topveg_k = k+1
EXIT
EXIT
ENDIF
END DO
END DO !k
ENDIF

! Vertical distribution function of qdm
Expand All @@ -682,8 +699,9 @@ SUBROUTINE COMPUTE_VEGDISS_VF_TERMS_SCHISM
DO k = kbs(IS), topveg_k-1
sum_3D = sum_3D + (swild_3D(k+1) + swild_3D(k))/2.D0*(zs(k+1,IS) - zs(k,IS))
END DO !topveg_k-1
IF(sum_3D == 0) CALL parallel_abort('Vertical profile in wave vegetation force: integral=0')
IF(sum_3D==0) CALL parallel_abort('WWM: Vertical profile in wave vegetation force: integral=0')

!'
DO k = kbs(IS), topveg_k
Fveg_x_loc = - swild_3D(k)*(SVEG(1,n1) + SVEG(1,n2))/2.D0/sum_3D
Fveg_y_loc = - swild_3D(k)*(SVEG(2,n1) + SVEG(2,n2))/2.D0/sum_3D
Expand Down Expand Up @@ -742,7 +760,7 @@ SUBROUTINE COMPUTE_INTRAWAVE_VEG_FORCE
!Initialization
WWAVE_FORCE_VEG_NL_TOT = 0.D0

Uorbi(:) = 0.D0
Uorbi(:) = 0.D0
etaw(:) = 0.D0

! Check IF dry segment or open bnd segment
Expand All @@ -752,14 +770,14 @@ SUBROUTINE COMPUTE_INTRAWAVE_VEG_FORCE
n1 = isidenode(1,IS); n2 = isidenode(2,IS)

eta_tmp = (eta2(n1) + eta2(n2))/2.D0
htot = max(h0,dps(IS)+eta_tmp)
htot = max(h0,dps(IS)+eta_tmp) !>0

!Vegetation characteristics at sides
SAV_CD_tmp = (VCD(n1) + VCD(n2))/2.D0
SAV_BV_tmp = (VDM(n1) + VDM(n2))/2.D0
SAV_NV_tmp = (VNV(n1) + VNV(n2))/2.D0
SAV_H_tmp = (VLTH(n1) + VLTH(n2))/2.D0

!Check IF this is a vegetated SIDE
IF (SAV_CD_tmp*SAV_BV_tmp*SAV_NV_tmp*SAV_H_tmp == 0.D0) CYCLE

Expand All @@ -774,30 +792,26 @@ SUBROUTINE COMPUTE_INTRAWAVE_VEG_FORCE
w_dir_tmp = DEG*PI/180.d0

IF (Hs_tmp <= 0.00D0) CYCLE
IF (TM10_tmp <= 0.00D0) CYCLE
IF (TM10_tmp <= 0.00D0) CYCLE

!Bed slope at sides
tanbeta_x_tmp = (tanbeta_x(n1) + tanbeta_x(n2))/2.D0
tanbeta_y_tmp = (tanbeta_y(n1) + tanbeta_y(n2))/2.D0

! Computation of wave orbital velocities and water surface elevations over a wave cycle in case of non linear waves.
! Both of the proposed methods compute near bed orbital velocity. As effects of coastal vegetation on waves usually
! occur in shallow water conditions, near bed orbital velocity is used as the depth-averaged value to a first approximation.
!The vegetation force is then empirically distributed over the vegetation height. The implementation of a more accurate
! method computing vertical profiles of orbital velocity is under development (See Fenton, 1988).

! Choose ONE of the following methods : WAVE_ASYMMETRY_ELFRINK_VEG or WAVE_ASYMMETRY_RF.

! The routine WAVE_ASYMMETRY_ELFRINK_VEG computes a time series of wave orbital velocity
! and water surface elevation over a wave cycle based on Elfrink et al. (2006).

!CALL WAVE_ASYMMETRY_ELFRINK_VEG(Hs_tmp,TM10_tmp,w_dir_tmp,htot,tanbeta_x_tmp,tanbeta_y_tmp,&
!ech+1,Uorbi,Ucrest,Utrough,T_crest,T_trough,etaw)

! The routine WAVE_ASYMMETRY_RF computes a time series of wave orbital velocity
! and water surface elevation over a wave cycle based on Rienecker and Fenton (1981) and Ruessink et al. (2012).

CALL WAVE_ASYMMETRY_RF(Hrms_tmp,TM10_tmp,klm_tmp,htot,ech,Uorbi,etaw)
CALL WAVE_ASYMMETRY_RF(Hrms_tmp,TM10_tmp,klm_tmp,htot,ech,Uorbi,etaw)

!Fv,w = 1/Trep int 0 -> Trep 0.5*rhow*Cd*bv*Nv*h'v*uw*abs(uw) dt
dt = TM10_tmp/(ech-1)
Expand All @@ -822,12 +836,12 @@ SUBROUTINE COMPUTE_INTRAWAVE_VEG_FORCE
topveg_k = NVRT
IF (SAV_H_tmp >= htot) THEN !emergent vegetation
topveg_k = NVRT
ELSE
ELSE !submerged vegetation
DO k = kbs(IS), NVRT-1
tmp0 = tmp0 + zs(k+1,IS) - zs(k,IS)
IF (tmp0 >= SAV_H_tmp) THEN !submerged vegetation
IF (tmp0 >= SAV_H_tmp) THEN
topveg_k = k+1
EXIT
EXIT
ENDIF
END DO
ENDIF
Expand All @@ -844,20 +858,19 @@ SUBROUTINE COMPUTE_INTRAWAVE_VEG_FORCE
sum_3D = sum_3D + (swild_3D(k+1) + swild_3D(k))/2.D0*(zs(k+1,IS) - zs(k,IS))
END DO !topveg_k-1

IF(sum_3D == 0) CALL parallel_abort('Vertical profile in intrawave vegetation force: integral=0')
IF(sum_3D==0) CALL parallel_abort('Vertical profile in intrawave vegetation force: integral=0')
!'

DO k = kbs(IS), topveg_k
WWAVE_FORCE_VEG_NL_x = - swild_3D(k)*WWAVE_FORCE_VEG_NL_TOT*COS(w_dir_tmp)/sum_3D
WWAVE_FORCE_VEG_NL_y = - swild_3D(k)*WWAVE_FORCE_VEG_NL_TOT*SIN(w_dir_tmp)/sum_3D

WWAVE_FORCE(1,k,IS) = WWAVE_FORCE(1,k,IS) + WWAVE_FORCE_VEG_NL_x
WWAVE_FORCE(2,k,IS) = WWAVE_FORCE(2,k,IS) + WWAVE_FORCE_VEG_NL_y
END DO
ENDIF !2D/3D
END DO

CALL exchange_s3d_2(WWAVE_FORCE)

END SUBROUTINE

!**********************************************************************
Expand Down
2 changes: 1 addition & 1 deletion src/WWMIII/wwm_datapl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ MODULE DATAPOOL
!diagnosis
& NE_GLOBAL => ne_global, &! Global number of elements
& NP_GLOBAL => np_global, &! Global number of nodes
& INETMP => elnode_wwm, & ! Element connection table of the augmented domain?
& INETMP => elnode_wwm, & ! Element connection table of the augmented domain (after splitting quads)
& iplg, & ! node local to global mapping
& ipgl, & ! node global to local mapping
! & ielg, & ! element local to global mapping
Expand Down

0 comments on commit d325ed3

Please sign in to comment.