Skip to content

Commit 6f81cda

Browse files
LeynseLeynse
committed
- Already do the bugfix of PR #283 by switching to srcig based on E and redistribute by /E * ee
1 parent 111f782 commit 6f81cda

1 file changed

Lines changed: 89 additions & 36 deletions

File tree

‎source/src/snapwave/snapwave_solver.f90‎

Lines changed: 89 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -487,7 +487,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
487487
!
488488
! Actual determining of source term:
489489
!
490-
call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
490+
call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
491491
!
492492
! inout: alphaig_local, srcig_local - eeprev, eeprev_ig, cgprev, beta_local
493493
! in: the rest
@@ -575,7 +575,7 @@ subroutine solve_energy_balance2Dstat(x,y,dhdx, dhdy, no_nodes,inner, &
575575
!
576576
! Actual determining of source term - every first sweep of iteration
577577
!
578-
call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
578+
call determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
579579
!
580580
endif
581581
!
@@ -1076,7 +1076,15 @@ subroutine baldock (rho, g, alfa, gamma, depth, H, T, iexp, Dw, Hmax)
10761076
!
10771077
end subroutine baldock
10781078

1079-
subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, cg_ig, nwav, depth, zb, H, ee, ee_ig, eeprev, eeprev_ig, cgprev, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
1079+
1080+
subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, ds, prev, dtheta, cg_ig, nwav, depth, zb, H, ee, ee_ig, ig_opt, alphaigfac, alphaig_local, beta_local, srcig_local)
1081+
!
1082+
! Determining of IG source term as defined in Leijnse et al. 2024
1083+
!
1084+
! inout: alphaig_local, srcig_local, beta_local
1085+
! in: the rest
1086+
!
1087+
! NOTE - This is based on the energy in the previous SnapWave timestep 'ee' and 'ee_ig', and waveheight 'H', which should therefore be made available.
10801088
!
10811089
implicit none
10821090
!
@@ -1095,6 +1103,7 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d
10951103
real*4, dimension(ntheta,no_nodes), intent(in) :: ee_ig ! energy density infragravity waves
10961104
integer, intent(in) :: ig_opt ! option of IG wave settings (1 = default = conservative shoaling based dSxx and Baldock breaking)
10971105
real*4, intent(in) :: alphaigfac ! Multiplication factor for IG shoaling source/sink term, default = 1.0
1106+
real*4, intent(in) :: dtheta ! directional resolution
10981107
!
10991108
! Inout variables
11001109
real*4, dimension(:,:), intent(inout) :: alphaig_local ! Local infragravity wave shoaling parameter alpha
@@ -1107,20 +1116,60 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d
11071116
integer :: itheta ! directional counter
11081117
integer :: k ! counters (k is grid index)
11091118
integer :: k1,k2 ! upwind counters (k is grid index)
1110-
real*4 :: gam ! local gamma (Hinc / depth ratio)
1111-
real*4, dimension(ntheta,no_nodes) :: depthprev ! water depth at upwind intersection point
1112-
real*4, dimension(ntheta,no_nodes) :: Sxx ! Radiation Stress
1113-
real*4, dimension(:), allocatable :: Sxxprev ! radiation stress at upwind intersection point
1114-
real*4, dimension(:), allocatable :: Hprev ! Incident wave height at upwind intersection point
1119+
real*4 :: gam ! local gamma (Hinc / depth ratio)
1120+
real*4, dimension(ntheta,no_nodes) :: depthprev ! water depth at upwind intersection point
1121+
real*4, dimension(no_nodes) :: Sxx ! radiation Stress
1122+
real*4, dimension(ntheta) :: Sxxprev ! radiation stress at upwind point
1123+
real*4, dimension(ntheta) :: Hprev ! wave height at upwind point
1124+
real*4, dimension(ntheta) :: cgprev ! group velocity at upwind point
1125+
real*4, dimension(ntheta) :: Eprev ! Mean incident wave energy at upwind intersection point
1126+
real*4, dimension(ntheta) :: Eprev_ig ! Mean infragravity wave energy at upwind intersection point
1127+
real*4, dimension(no_nodes) :: E_local ! mean wave energy waves - just local
1128+
real*4, dimension(no_nodes) :: E_ig_local ! mean wave energy infragravity waves - just local
11151129
real*4 :: dSxx ! difference in Radiation stress
1116-
real*4 :: Sxx_cons ! conservative estimate of radiation stress using conservative shoaling
1117-
!
1118-
! Allocate internal variables
1119-
allocate(Sxxprev(ntheta))
1120-
allocate(Hprev(ntheta))
1130+
real*4 :: Sxx_cons ! conservative estimate of radiation stress using conservative shoaling
1131+
!
1132+
! Set internal variables
11211133
!
11221134
Sxx = 0.0
1135+
Hprev = 0.0
1136+
Eprev = 0.0
1137+
Eprev_ig = 0.0
1138+
!
1139+
E_local = 0.0
1140+
E_ig_local = 0.0
1141+
Sxx = 0.0
1142+
! Pre-compute Sxx for all nodes
1143+
!
1144+
!$omp parallel do schedule(static)
1145+
do k = 1, no_nodes
1146+
!
1147+
if (inner(k)) then !TODO: check whether should be on only 'inner' or not
1148+
!
1149+
! Update E (not saved from previous timestep)
1150+
!
1151+
E_local(k) = sum(ee(:,k)) * dtheta
1152+
!
1153+
! Update E_ig (not saved from previous timestep)
1154+
!
1155+
E_ig_local(k) = sum(ee_ig(:, k)) * dtheta
1156+
!
1157+
endif
1158+
!
1159+
Sxx(k) = ((2.0 * max(0.0, min(1.0, nwav(k)))) - 0.5) * E_local(k)
1160+
!
1161+
enddo
1162+
!$omp end parallel do
1163+
!
1164+
! Main loop: compute IG source/sink term per node.
1165+
! All writes target the column (itheta, k), so the loop is data-independent across k.
1166+
! Per-k scratch arrays (cgprev, Eprev, Eprev_ig, Sxxprev, Hprev) are
1167+
! listed as private so each thread gets its own copy on the stack.
11231168
!
1169+
!$omp parallel do &
1170+
!$omp& private(itheta, k1, k2, gam, dSxx, Sxx_cons, &
1171+
!$omp& cgprev, Eprev, Eprev_ig, Sxxprev, Hprev) &
1172+
!$omp& schedule(static)
11241173
do k = 1, no_nodes
11251174
!
11261175
if (inner(k)) then
@@ -1152,10 +1201,10 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d
11521201
Sxx(itheta,k1) = ((2.0 * max(0.0,min(1.0,nwav(k1)))) - 0.5) * ee(itheta, k1) ! limit so value of nwav is between 0 and 1
11531202
Sxx(itheta,k2) = ((2.0 * max(0.0,min(1.0,nwav(k2)))) - 0.5) * ee(itheta, k2) ! limit so value of nwav is between 0 and 1
11541203
!
1155-
Sxxprev(itheta) = w(1, itheta, k)*Sxx(itheta,k1) + w(2, itheta, k)*Sxx(itheta,k2)
1204+
Sxxprev(itheta) = w(1, itheta, k) * Sxx(k1) + w(2, itheta, k) * Sxx(k2)
11561205
!
1157-
eeprev(itheta) = w(1, itheta, k)*ee(itheta, k1) + w(2, itheta, k)*ee(itheta, k2)
1158-
eeprev_ig(itheta) = w(1, itheta, k)*ee_ig(itheta, k1) + w(2, itheta, k)*ee_ig(itheta, k2)
1206+
Eprev(itheta) = w(1, itheta, k) * E_local(k1) + w(2, itheta, k) * E_local(k2)
1207+
Eprev_ig(itheta) = w(1, itheta, k) * E_ig_local(k1) + w(2, itheta, k) * E_ig_local(k2)
11591208
!
11601209
Hprev(itheta) = w(1, itheta, k)*H(k1) + w(2, itheta, k)*H(k2)
11611210
!
@@ -1179,26 +1228,30 @@ subroutine determine_infragravity_source_sink_term(inner, no_nodes, ntheta, w, d
11791228
srcig_local(itheta, k) = 0.0 !Avoid big jumps in dSxx that can happen if a upwind point is a boundary point with Hinc=0
11801229
!
11811230
else
1182-
!
1183-
if (ig_opt == 1) then ! Option using conservative shoaling for dSxx/dx
1184-
!
1185-
! Calculate Sxx based on conservative shoaling of upwind point's energy:
1186-
! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5)
1187-
Sxx_cons = eeprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0,min(1.0,nwav(k)))) - 0.5)
1188-
! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite
1189-
!
1190-
dSxx = Sxx_cons - Sxxprev(itheta)
1191-
!
1192-
elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx
1193-
!
1194-
dSxx = Sxx(itheta,k) - Sxxprev(itheta)
1195-
!
1196-
endif
1197-
!
1198-
dSxx = max(dSxx, 0.0)
1199-
!
1200-
srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(eeprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k)
1201-
!
1231+
!
1232+
if (ig_opt == 1) then ! Option using conservative shoaling for dSxx/dx
1233+
!
1234+
! Calculate Sxx based on conservative shoaling of upwind point's energy:
1235+
! Sxx_cons = E(i-1) * Cg(i-1) / Cg * (2 * n(i) - 0.5)
1236+
!
1237+
Sxx_cons = Eprev(itheta) * cgprev(itheta) / cg_ig(k) * ((2.0 * max(0.0, min(1.0, nwav(k)))) - 0.5)
1238+
!
1239+
! Note - limit so value of nwav is between 0 and 1, and Sxx therefore doesn't become NaN for nwav=Infinite
1240+
!
1241+
dSxx = Sxx_cons - Sxxprev(itheta)
1242+
!
1243+
elseif (ig_opt == 2) then ! Option taking actual difference for dSxx/dx
1244+
!
1245+
dSxx = Sxx(k) - Sxxprev(itheta)
1246+
!
1247+
endif
1248+
!
1249+
dSxx = max(dSxx, 0.0)
1250+
!
1251+
! Base on E_prev_ig instead of eeprev_ig(itheta) > no bins but total energy
1252+
! NOTE - already here multiplied with ee(itheta,k), for direct inclusion in 'R'-term
1253+
srcig_local(itheta, k) = alphaigfac * alphaig_local(itheta,k) * sqrt(Eprev_ig(itheta)) * cgprev(itheta) / depthprev(itheta,k) * dSxx / ds(itheta, k) /max(E_local(k), 1.0e-6) * ee(itheta,k)
1254+
!
12021255
endif
12031256
!
12041257
else ! TL: option to add future parameterisations here for e.g. coral reef type coasts

0 commit comments

Comments
 (0)