From 3696fa08aeafc0638eb95cb1372b2e12994192b6 Mon Sep 17 00:00:00 2001
From: lucieb92 <lucie.bakels@univie.ac.at>
Date: Wed, 1 May 2024 17:19:57 +0200
Subject: [PATCH] Set maximum for settling velocity: it cannot be larger than
the height of the particle to avoid settling reflection
Former-commit-id: 25284cb5001652da5d55b14d3098613812b45e50
---
src/advance_mod.f90 | 17 ++++++++++++++++-
1 file changed, 16 insertions(+), 1 deletion(-)
diff --git a/src/advance_mod.f90 b/src/advance_mod.f90
index ee8963b0..98e94423 100644
--- a/src/advance_mod.f90
+++ b/src/advance_mod.f90
@@ -410,10 +410,16 @@ subroutine adv_above_pbl(itime,itimec,dxsave,dysave,ux,vy,tropop,nrand,ipart)
call get_settling(xts,yts,zts,nsp,part(ipart)%settling)
#ifdef ETA
call update_zeta_to_z(itime,ipart)
+ if (part(ipart)%z+part(ipart)%settling*dt.lt.0) then
+ part(ipart)%settling=-part(ipart)%z/dt
+ endif
call w_to_weta(itime,dt,part(ipart)%xlon,part(ipart)%ylat, &
part(ipart)%z,part(ipart)%zeta,part(ipart)%settling,weta_settling)
- weta=weta+weta_settling
+ weta=weta+weta_settling
#else
+ if (part(ipart)%z+part(ipart)%settling*dt.lt.0) then
+ part(ipart)%settling=-part(ipart)%z/dt
+ endif
w=w+part(ipart)%settling
#endif
end if
@@ -586,6 +592,9 @@ subroutine adv_in_pbl(itime,itimec, dxsave,dysave,dawsave,dcwsave, abovePBL, &
endif
if (density(nsp).gt.0.) then
call get_settling(xts,yts,zts,nsp,part(ipart)%settling) !bugfix
+ if (part(ipart)%z+part(ipart)%settling*dt.lt.0) then
+ part(ipart)%settling=-part(ipart)%z/dt
+ endif
w=w+part(ipart)%settling
end if
end if
@@ -712,6 +721,9 @@ subroutine petterssen_corr(itime,ipart)
call update_z_to_zeta(itime+part(ipart)%idt,ipart)
zts=real(part(ipart)%z)
call get_settling(xts,yts,zts,nsp,part(ipart)%settling) !bugfix
+ if (part(ipart)%z+part(ipart)%settling*part(ipart)%idt.lt.0) then
+ part(ipart)%settling=-part(ipart)%z/part(ipart)%idt
+ endif
call w_to_weta( &
itime+part(ipart)%idt, real(part(ipart)%idt), part(ipart)%xlon, &
part(ipart)%ylat, part(ipart)%z, part(ipart)%zeta, &
@@ -721,6 +733,9 @@ subroutine petterssen_corr(itime,ipart)
!real(part(ipart)%zeta-part(ipart)%zeta_prev)/real(part(ipart)%idt*ldirect)
#else
call get_settling(xts,yts,zts,nsp,part(ipart)%settling)
+ if (part(ipart)%z+part(ipart)%settling*part(ipart)%idt.lt.0) then
+ part(ipart)%settling=-part(ipart)%z/part(ipart)%idt
+ endif
w=w+part(ipart)%settling
#endif
end if
--
GitLab