Opened 5 years ago

Last modified 4 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 5 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

Corrected code supplied

comment:2 Changed 5 years ago by pesei

  • Version none deleted

comment:3 Changed 5 years ago by hasod

  • Status changed from assigned to accepted

comment:4 Changed 5 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 5 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 4 years ago by hasod

  • Resolution set to fixed
  • Status changed from accepted to closed

implemented as suggested

comment:7 Changed 4 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/)

Note: See TracTickets for help on using tickets.
hosted by ZAMG