Opened 7 years ago
Last modified 6 years ago
#54 reopened Defect
Particle loss at poles
Reported by: | alaedera | Owned by: | hasod |
---|---|---|---|
Priority: | major | Milestone: | FLEXPART 9.2 |
Component: | FP physics/numerics | Version: | |
Keywords: | Cc: |
Description
Running FLEXPART 9.1 in domainfilling global mode results in loosing particles at the poles.
Change History (7)
comment:1 Changed 7 years ago by hasod
- Component changed from FP coding/compilation to FP physics/numerics
- Owner changed from somebody to hasod
- Status changed from new to assigned
comment:2 Changed 7 years ago by pesei
- Version none deleted
Is this related to ticket:39 (Interpolation close to the poles)?
comment:3 Changed 7 years ago by hasod
- Status changed from assigned to accepted
comment:4 Changed 7 years ago by alaedera
@pesei
I'm not aware about the discussion results of the "FLEXPART Developer Workshop February 2012" leading to the ticket #39 (Interpolation close to the poles) as I started my PhD later on.
Anyway, the particle loss has been a consequence of missing handling of particles passing the poles and some .ge. versus .gt. checks in the advance.f90 routine.
I hope this information might help.
Best regards,
Alex
comment:5 Changed 7 years ago by hasod
In file advance.f90, lines 711 change the code as shown below (block BUGFIX)
! If global data are available, use cyclic boundary condition !************************************************************ if (xglobal) then if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1) if (xt.lt.0.) xt=xt+real(nxmin1) if (xt.le.eps) xt=eps if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps endif ! --- BUGFIX BEGIN --- if ( yt.lt.0. ) then xt=mod(xt+180.,360.) yt=-yt else if ( yt.gt.real(nymin1) ) then xt=mod(xt+180.,360.) yt=2*real(nymin1)-yt endif ! Check position: If trajectory outside model domain, terminate it !***************************************************************** if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. & (yt.gt.real(nymin1))) then !BUGFIX: changed .ge. to .gt. nstop=3 return endif ! --- BUGFIX END ---
Very similar changes are required in advance.f90, lines 881 and below (section indicated with BUGFIX):
! If global data are available, use cyclic boundary condition !************************************************************ if (xglobal) then if (xt.ge.real(nxmin1)) xt=xt-real(nxmin1) if (xt.lt.0.) xt=xt+real(nxmin1) if (xt.le.eps) xt=eps if (abs(xt-real(nxmin1)).le.eps) xt=real(nxmin1)-eps endif ! --- BUGFIX BEGIN --- if ( yt.lt.0. ) then xt=mod(xt+180.,360.) yt=-yt else if ( yt.gt.real(nymin1) ) then xt=mod(xt+180.,360.) yt=2*real(nymin1)-yt endif ! Check position: If trajectory outside model domain, terminate it !***************************************************************** if ((xt.lt.0.).or.(xt.ge.real(nxmin1)).or.(yt.lt.0.).or. & (yt.gt.real(nymin1))) then !BUGFIX: changed .ge. to .gt. nstop=3 endif ! --- BUGFIX END ---
comment:6 Changed 6 years ago by hasod
- Resolution set to fixed
- Status changed from accepted to closed
implemented as suggested
comment:7 Changed 6 years ago by pesei
- Resolution fixed deleted
- Status changed from closed to reopened
I think this fix is working properly only for global domains.
As I would like to understand better the exact condition where the problem appears, I'd ask alaedera to post details or to please contact me directly (http://homepage.univie.ac.at/petra.seibert/)
Corrected code supplied