Particle loss at poles

[alaedera]

Running FLEXPART 9.1 in domainfilling global mode results in loosing particles at the poles.

[hasod]

Corrected code supplied

[pesei]

Is this related to ticket:39 (Interpolation close to the poles)?

[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

[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 ---

[hasod]

implemented as suggested

[pesei]

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/)