Modules split into source files - test version

src/adv_adv_above_pbl.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine adv_above_pbl(itime,itimec,dxsave,dysave,ux,vy,tropop,nrand,ipart)
+
+  implicit none
+  integer, intent(in) ::          &
+    itime,                        & ! time index
+    ipart                           ! particle index
+  integer, intent(inout) ::       &
+    itimec,                       & ! next timestep
+    nrand                           ! random number used for turbulence
+  real, intent(in) ::             &
+    tropop                          ! height of troposphere
+  real, intent(inout) ::          &
+    ux,vy,                        & ! random turbulent velocities above PBL
+    dxsave,dysave                   ! accumulated displacement in long and lat
+  real ::                         &
+    dt,                           & ! real(ldt)
+    xts,yts,zts,ztseta,           & ! local 'real' copy of the particle position
+    weta_settling,                & ! settling velocity in eta coordinates
+    wp                              ! random turbulence velocities
+  integer ::                      &
+    insp,nsp                        ! loop variables for number of species
+
+  zts=real(part(ipart)%z)
+  ztseta=real(part(ipart)%zeta)
+  xts=real(part(ipart)%xlon)
+  yts=real(part(ipart)%ylat)
+  if (lsettling) part(ipart)%settling=0.
+
+  call interpol_wind(itime,xts,yts,zts,ztseta,ipart)
+
+  ! Compute everything for above the PBL
+
+  ! Assume constant, uncorrelated, turbulent perturbations
+  ! In the stratosphere, use a small vertical diffusivity d_strat,
+  ! in the troposphere, use a larger horizontal diffusivity d_trop.
+  ! Turbulent velocity scales are determined based on sqrt(d_trop/dt)
+  !******************************************************************
+
+  part(ipart)%idt=abs(lsynctime-itimec+itime)
+  dt=real(part(ipart)%idt)
+
+  if (.not.turboff) then
+    call turb_stratosph(dt,nrand,ux,vy,wp,tropop,zts)
+  else
+    !sec switch off turbulence
+    ux=0.0
+    vy=0.0
+    wp=0.0
+  endif
+
+  ! If particle represents only a single species, add gravitational settling
+  ! velocity. The settling velocity is zero for gases
+  !*************************************************************************
+  ! Does not work in eta coordinates yet
+  if (mdomainfill.eq.0) then
+    if (lsettling) then
+      if ((ipin.ne.3).and.(ipin.ne.4)) then
+        do insp=1,nspec
+          nsp=insp
+          if (xmass(part(ipart)%npoint,nsp).gt.eps3) exit
+        end do
+      else
+        nsp=1
+      endif
+      ! LB change to eta coords?
+      if (density(nsp).gt.0.) then
+        call get_settling(itime,xts,yts,zts,nsp,part(ipart)%settling)
+        select case (wind_coord_type)
+          case ('ETA')
+            call update_zeta_to_z(itime,ipart)
+            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
+          case ('METER')
+            w=w+part(ipart)%settling
+          case default
+            w=w+part(ipart)%settling
+        end select
+      end if
+    endif
+  end if
+
+  ! Calculate position at time step itime+lsynctime
+  !************************************************
+  dxsave=dxsave+(u+ux)*dt
+  dysave=dysave+(v+vy)*dt
+
+  select case (wind_coord_type)
+
+    case ('ETA')
+
+      if (wp.ne.0.) then
+        call update_zeta_to_z(itime,ipart)
+        call update_z(ipart,wp*dt*real(ldirect))
+        if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+          ! if particle below ground -> reflection
+        call update_z_to_zeta(itime,ipart)
+      endif
+      call update_zeta(ipart,weta*dt*real(ldirect))
+      if (part(ipart)%zeta.ge.1.) call set_zeta(ipart,1.-(part(ipart)%zeta-1.))
+      if (part(ipart)%zeta.eq.1.) call update_zeta(ipart,-eps_eta)
+
+    case ('METER')
+
+      call update_z(ipart,(w+wp)*dt*real(ldirect))
+      if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+
+    case default
+
+      call update_z(ipart,(w+wp)*dt*real(ldirect))
+      if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+
+  end select
+
+end subroutine adv_above_pbl

src/adv_adv_in_pbl.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine adv_in_pbl(itime,itimec, dxsave,dysave,dawsave,dcwsave, abovePBL,  &
+  nrand,ipart,thread)
+
+  use drydepo_mod, only: drydepo_probability
+
+  implicit none
+
+  logical, intent(inout) ::       &
+    abovePBL
+  ! flag will be set to 'true' if computation needs to be completed above PBL
+  integer, intent(in) ::          &
+    itime,                        & ! time index
+    ipart,                        & ! particle index
+    thread                          ! number of the omp thread
+  real, intent(inout) ::          &
+    dxsave,dysave,                & ! accumulated displacement in long and lat
+    dawsave,dcwsave               ! accumulated displacement in wind directions
+  integer, intent(inout) ::       &
+    itimec,                       & ! next timestep
+    nrand                           ! random number used for turbulence
+  real ::                         &
+    dt,                           & ! real(ldt)
+    xts,yts,zts,ztseta,           & ! local 'real' copy of the particle position
+    rhoa,                         & ! air density, used in CBL
+    rhograd                     ! vertical gradient of air density, used in CBL
+  integer ::                      &
+    loop,                         & ! loop variable for time in the PBL
+    nsp,insp                        ! loop variable for species
+  real :: vdepo(maxspec)  ! deposition velocities for all species
+
+  eps=nxmax/3.e5
+  if (lsettling) part(ipart)%settling=0.
+
+  ! BEGIN TIME LOOP
+  !================
+  ! For wind_coord_type=ETA:
+  ! Within this loop, only METER coordinates are used, and the new z value will
+  ! be updated to ETA coordinates at the end
+  !****************************************************************************
+
+  call update_zeta_to_z(itime,ipart)
+
+  loop=0
+  pbl_loop: do
+
+    loop=loop+1
+    if (method.eq.1) then
+      part(ipart)%idt=min(part(ipart)%idt,abs(lsynctime-itimec+itime))
+      itimec=itimec+part(ipart)%idt*ldirect
+    else
+      part(ipart)%idt=abs(lsynctime)
+      itimec=itime+lsynctime
+    endif
+    dt=real(part(ipart)%idt)
+    xts=real(part(ipart)%xlon)
+    yts=real(part(ipart)%ylat)
+    zts=real(part(ipart)%z)
+
+    zeta=zts/h
+    if (loop.eq.1) then ! Temporal interpolation only for the first iteration
+
+      if (ngrid.le.0) then
+        xts=real(part(ipart)%xlon)
+        yts=real(part(ipart)%ylat)
+        call interpol_pbl(itime,xts,yts,zts,real(part(ipart)%zeta))
+      else
+        call interpol_pbl(itime,xtn,ytn,zts,real(part(ipart)%zeta))
+      endif
+
+    else
+
+      ! Determine the level below the current position for u,v,rho
+      !***********************************************************
+      call find_z_level_meters(zts)
+
+      ! If one of the levels necessary is not yet available,
+      ! calculate it
+      !*****************************************************
+      call interpol_pbl_misslev()
+
+    endif
+
+  ! Vertical interpolation of u,v,w,rho and drhodz
+  !***********************************************
+
+  ! Vertical distance to the level below and above current position
+  ! both in terms of (u,v) and (w) fields
+  !****************************************************************
+
+    call interpol_pbl_short(zts,rhoa,rhograd) ! Vertical interpolation
+
+  ! Compute the turbulent disturbances
+  ! Determine the sigmas and the timescales
+  !****************************************
+
+    if (.not.turboff) then
+      call turb_pbl(ipart,nrand,dt,zts,rhoa,rhograd,thread)
+      ! Note: zts and nrand get updated
+
+      ! Determine time step for next integration
+      !*****************************************
+      if (turbswitch) then
+        part(ipart)%idt = int( &
+          min( tlw, &
+               h/max( 2.*abs(part(ipart)%turbvel%w*sigw), 1.e-5 ), &
+               0.5/abs(dsigwdz) &
+             ) *ctl)
+      else
+        part(ipart)%idt = int( &
+          min( tlw, &
+               h/max( 2.*abs(part(ipart)%turbvel%w), 1.e-5) &
+              ) *ctl)
+      endif
+    else
+      part(ipart)%turbvel%u=0.
+      part(ipart)%turbvel%v=0.
+      part(ipart)%turbvel%w=0.
+    endif
+
+    part(ipart)%idt=max(part(ipart)%idt,mintime)
+
+
+  ! If particle represents only a single species, add gravitational settling
+  ! velocity. The settling velocity is zero for gases, or if particle
+  ! represents more than one species
+  !*************************************************************************
+
+    if (mdomainfill.eq.0) then
+      if (lsettling) then
+        if ((ipin.ne.3).and.(ipin.ne.4)) then
+          do insp=1,nspec
+            nsp=insp
+            if (xmass(part(ipart)%npoint,nsp).gt.eps3) exit
+          end do
+        else
+          nsp=1
+        endif
+        if (density(nsp).gt.0.) then
+          call get_settling(itime,xts,yts,zts,nsp,part(ipart)%settling)  !bugfix
+          w=w+part(ipart)%settling
+        end if
+      end if
+    endif
+
+  ! Horizontal displacements during time step dt are small real values compared
+  ! to the position; adding the two, would result in large numerical errors.
+  ! Thus, displacements are accumulated during lsynctime and are added to the
+  ! position at the end
+  !****************************************************************************
+
+    dxsave=dxsave+u*dt
+    dysave=dysave+v*dt
+    dawsave=dawsave+part(ipart)%turbvel%u*dt
+    dcwsave=dcwsave+part(ipart)%turbvel%v*dt
+    ! How can I change the w to w(eta) efficiently?
+
+    select case (wind_coord_type)
+
+      case ('ETA')
+
+        call update_z(ipart,w*dt*real(ldirect))
+        zts=real(part(ipart)%z)
+        ! HSO/AL: Particle managed to go over highest level -> interpolation
+        ! error in goto 700
+        !          alias interpol_wind (division by zero)
+        if (zts.ge.height(nz)) call set_z(ipart,height(nz)-100.*eps)
+         ! Manually for z instead
+
+      case ('METER')
+
+        call update_z(ipart,w*dt*real(ldirect))
+        call pushpartdown(ipart)
+
+    end select
+    zts=real(part(ipart)%z)
+
+    if (zts.gt.h) then
+      call update_z_to_zeta(itime,ipart)
+      if (itimec.ne.itime+lsynctime) abovePBL=.true.
+        ! complete the current interval above PBL
+      return
+    endif
+
+  ! Determine probability of deposition
+  !************************************
+    call drydepo_probability(part(ipart)%prob,dt,zts,vdepo)
+
+    if (zts.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+      ! if particle below ground -> reflection
+
+    if (itimec.eq.(itime+lsynctime)) then
+      ! Convert z position that changed by turbulent motions to eta coords
+      call update_z_to_zeta(itime,ipart)
+      return  ! finished
+    endif
+
+  end do pbl_loop
+
+  call update_z_to_zeta(itime,ipart)
+
+end subroutine adv_in_pbl

src/adv_advance.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine advance(itime,ipart,thread)
+
+  !*****************************************************************************
+  !                                                                            *
+  !  Calculation of turbulent particle trajectories utilizing a                *
+  !  zero-acceleration scheme, which is corrected by a numerically more        *
+  !  accurate Petterssen scheme whenever possible.                             *
+  !                                                                            *
+  !  Particle positions are read in, incremented, and returned to the calling  *
+  !  program.                                                                  *
+  !                                                                            *
+  !  In different regions of the atmosphere (PBL vs. free troposphere),        *
+  !  different parameters are needed for advection, parameterizing turbulent   *
+  !  velocities, etc. For efficiency, different interpolation routines have    *
+  !  been written for these different cases, with the disadvantage that there  *
+  !  exist several routines doing almost the same. They all share the          *
+  !  included file 'interpol_mod'. The following                               *
+  !  interpolation routines are used:                                          *
+  !                                                                            *
+  !  interpol_all(_nest)     interpolates everything (called inside the PBL)  *
+  !  interpol_misslev(_nest) if a particle moves vertically in the PBL,       *
+  !                           additional parameters are interpolated if it     *
+  !                           crosses a model level                            *
+  !  interpol_wind(_nest)    interpolates the wind and determines the         *
+  !                           standard deviation of the wind (called outside   *
+  !                           PBL) also interpolates potential vorticity       *
+  !  interpol_wind_short(_nest) only interpolates the wind (needed for the    *
+  !                           Petterssen scheme)                               *
+  !  interpol_vdep(_nest)    interpolates deposition velocities               *
+  !                                                                            *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     16 December 1997                                                       *
+  !                                                                            *
+  !  Changes:                                                                  *
+  !                                                                            *
+  !  8 April 2000: Deep convection parameterization                            *
+  !                                                                            *
+  !  May 2002: Petterssen scheme introduced                                    *
+  !                                                                            *
+  !  2021, L. Bakels:                                                          *
+  !         - Separated PBL and above PBL computations in different            *
+  !           subroutines                                                      *
+  !         - Moved all turbulence computations to turbulence_mod.f90          *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! icbt               1 if particle not transferred to forbidden state,       *
+  !                    else -1                                                 *
+  ! dawsave            accumulated displacement in along-wind direction        *
+  ! dcwsave            accumulated displacement in cross-wind direction        *
+  ! dxsave             accumulated displacement in longitude                   *
+  ! dysave             accumulated displacement in latitude                    *
+  ! h [m]              Mixing height                                           *
+  ! lwindinterv [s]    time interval between two wind fields                   *
+  ! itime [s]          time at which this subroutine is entered                *
+  ! itimec [s]         actual time, which is incremented in this subroutine    *
+  ! href [m]           height for which dry deposition velocity is calculated  *
+  ! ladvance [s]       Total integration time period                           *
+  ! ldirect            1 forward, -1 backward                                  *
+  ! ldt [s]            Time step for the next integration                      *
+  ! lsynctime [s]      Synchronisation interval of FLEXPART                    *
+  ! ngrid              index which grid is to be used                          *
+  ! nrand              index for a variable to be picked from rannumb          *
+  ! nstop              if > 1 particle has left domain and must be stopped     *
+  ! prob               probability of absorption due to dry deposition         *
+  ! rannumb(maxrand)   normally distributed random variables                   *
+  ! rhoa               air density                                             *
+  ! rhograd            vertical gradient of the air density                    *
+  ! up,vp,wp           random velocities due to turbulence (along wind, cross  *
+  !                    wind, vertical wind                                     *
+  ! usig,vsig,wsig     mesoscale wind fluctuations                             *
+  ! xt,yt,zt           Particle position                                       *
+  !                                                                            *
+  !*****************************************************************************
+
+  ! openmp change
+  use omp_lib, only: OMP_GET_THREAD_NUM
+  ! openmp change end
+
+  implicit none
+  integer, intent(in) ::          &
+    itime,                        & ! time index
+    ipart,                        & ! particle index
+    thread                          ! OMP thread
+  integer ::                      &
+    itimec,                       &
+    i,j,k,                        & ! loop variables
+    nrand,                        & ! random number used for turbulence
+    memindnext,                   & ! seems useless
+    ngr,                          & ! temporary new grid index of moved particle
+    nsp                             ! loop variables for number of species
+  real ::                         &
+    ux,vy,                        & ! random turbulent velocities above PBL
+    weta_settling,                & ! Settling velocity in eta coordinates
+    tropop,                       & ! height of troposphere
+    dxsave,dysave,                & ! accumulated displacement in long and lat
+    dawsave,dcwsave                ! accumulated displacement in wind directions
+  logical ::                      &
+    abovePBL
+    ! flag will be set to 'true' if computation needs to be completed above PBL
+
+  eps=nxmax/3.e5
+
+  part(ipart)%nstop=.false.
+  do i=1,nmixz
+    indzindicator(i)=.true.
+  end do
+
+  if (DRYDEP) then    ! reset probability for deposition
+    depoindicator=.true.
+    part(ipart)%prob=0.
+  endif
+
+  if (lsettling) part(ipart)%settling=0.
+
+  !if (ipart.eq.1) write(*,*) 'Mass: ', part(ipart)%mass(:), itime
+  dxsave=0.           ! reset position displacements
+  dysave=0.           ! due to mean wind
+  dawsave=0.          ! and turbulent wind
+  dcwsave=0.
+
+  itimec=itime
+
+  nrand=int(ran3(iseed1(thread),thread)*real(maxrand-1))+1
+
+  ! Determine whether lat/long grid or polarstereographic projection
+  ! is to be used
+  ! Furthermore, determine which nesting level to be used
+  !*****************************************************************
+  call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)
+
+  !***************************
+  ! Interpolate necessary data
+  !***************************
+
+  if (abs(itime-memtime(1)).lt.abs(itime-memtime(2))) then
+    memindnext=1
+  else
+    memindnext=2
+  endif
+
+  ! Convert z(eta) to z(m) for the turbulence scheme, w(m/s)
+  ! is computed in verttransform_ecmwf.f90
+
+  call update_zeta_to_z(itime,ipart)
+
+  ! Determine nested grid coordinates
+  ! Determine the lower left corner and its distance to the current position
+  ! Calculate variables for time interpolation
+  !*******************************************
+  call init_interpol(itime, &
+    real(part(ipart)%xlon),real(part(ipart)%ylat),&
+    real(part(ipart)%z),   real(part(ipart)%zeta))
+
+  ! Compute maximum mixing height around particle position
+  !*******************************************************
+
+  ! Compute height of troposphere and PBL at x-y location of particle
+  call interpol_htropo_hmix(tropop,h)
+  zeta=real(part(ipart)%z)/h
+
+  !*************************************************************
+  ! If particle is in the PBL, interpolate once and then make a
+  ! time loop until end of interval is reached
+  !*************************************************************
+  ! In the PBL we use meters instead of eta coordinates for vertical transport
+
+  abovePBL=.true.
+  if (zeta.le.1.) then
+
+    abovePBL=.false.
+    call adv_in_pbl(itime,itimec,&
+      dxsave,dysave,dawsave,dcwsave,abovePBL,nrand,ipart,thread)
+    if (wind_coord_type.eq.'ETA' .and. lsettling) then
+      call w_to_weta(itime,real(part(ipart)%idt),part(ipart)%xlon, &
+        part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta, &
+        part(ipart)%settling,weta_settling)
+      weta=weta+weta_settling
+    endif
+
+  endif
+
+  !**********************************************************
+  ! For all particles that are outside the PBL, make a single
+  ! time step. Only horizontal turbulent disturbances are
+  ! calculated. Vertical disturbances are reset.
+  !**********************************************************
+
+  ! Interpolate the wind
+  !*********************
+
+  if (abovePBL) call adv_above_pbl(itime,itimec,dxsave,dysave, &
+    ux,vy,tropop,nrand,ipart)
+    ! Above PBL computation
+
+  !****************************************************************
+  ! Add mesoscale random disturbances
+  ! This is done only once for the whole lsynctime interval to save
+  ! computation time
+  !****************************************************************
+
+  ! Mesoscale wind velocity fluctuations are obtained by scaling
+  ! with the standard deviation of the grid-scale winds surrounding
+  ! the particle location, multiplied by a factor fturbmeso.
+  ! The autocorrelation time constant is taken as half the
+  ! time interval between wind fields
+  !****************************************************************
+
+  if (.not. turboff) then
+    ! mesoscale turbulence is found to give issues, so turned off
+    if (lmesoscale_turb) then
+      call interpol_mesoscale(itime, &
+        real(part(ipart)%xlon),real(part(ipart)%ylat), &
+        real(part(ipart)%z),real(part(ipart)%zeta))
+      call turb_mesoscale(nrand,dxsave,dysave,ipart, &
+        usig,vsig,wsig,wsigeta,eps_eta)
+    endif
+
+    !*************************************************************
+    ! Transform along and cross wind components to xy coordinates,
+    ! add them to u and v, transform u,v to grid units/second
+    ! and calculate new position
+    !*************************************************************
+
+    call windalign(dxsave,dysave,dawsave,dcwsave,ux,vy)
+    dxsave=dxsave+ux
+     ! comment by MC: comment this line to stop particles horizontally for tests
+    dysave=dysave+vy
+  endif
+
+  call update_xy(dxsave,dysave,ipart)
+  if (part(ipart)%nstop) return
+
+  ! If particle above highest model level, set it back into the domain
+  !*******************************************************************
+  call pushpartdown(ipart)
+
+  !************************************************************************
+  ! Now we could finish, as this was done in FLEXPART versions up to 4.0.
+  ! However, truncation errors of the advection can be significantly
+  ! reduced by doing one iteration of the Petterssen scheme, if this is
+  ! possible.
+  ! Note that this is applied only to the grid-scale winds, not to
+  ! the turbulent winds.
+  !************************************************************************
+
+  ! The Petterssen scheme can only applied with long time steps (only then u
+  ! is the "old" wind as required by the scheme); otherwise do nothing
+  !*************************************************************************
+
+  if (part(ipart)%idt .ne. abs(lsynctime)) return
+
+ ! The Petterssen scheme can only be applied if the ending time of the time step
+ ! (itime+ldt*ldirect) is still between the two wind fields held in memory;
+ ! otherwise do nothing
+ !******************************************************************************
+
+  if (abs(itime+part(ipart)%idt*ldirect).gt.abs(memtime(2))) return
+
+  ! Apply it also only if starting and ending point of current time step are on
+  ! the same grid; otherwise do nothing
+  !*****************************************************************************
+  ! ngr = ngrid
+  ! call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)
+
+  if (nglobal .and. real(part(ipart)%ylat).gt.switchnorthg) then
+    ngr=-1
+  else if (sglobal.and. real(part(ipart)%ylat).lt.switchsouthg) then
+    ngr=-2
+  else
+    ngr=0
+    ! Temporary fix for nested layer edges: replaced eps with dxn and dyn (LB)
+    do j=numbnests,1,-1
+      if (real(part(ipart)%xlon).gt.xln(j)+dxn(j) .and. &
+          real(part(ipart)%xlon).lt.xrn(j)-dxn(j) .and. &
+          real(part(ipart)%ylat).gt.yln(j)+dyn(j) .and. &
+          real(part(ipart)%ylat).lt.yrn(j)-dyn(j)) then
+        ngr=j
+        exit
+      endif
+    end do
+  endif
+
+  if (ngr.ne.ngrid) return
+
+  call petterson_corr(itime,ipart)
+
+end subroutine advance

src/adv_petterson_corr.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine petterson_corr(itime,ipart)
+
+  implicit none
+
+  integer, intent(in) ::          &
+    itime,                        & ! time index
+    ipart                           ! particle index
+  integer ::                      &
+    nsp,insp                        ! loop variables for number of species
+  real ::                         &
+    xts,yts,zts,ztseta,           & ! local 'real' copy of the particle position
+    uold,vold,wold,woldeta,       &
+    weta_settling
+  real(kind=dp) ::                &
+    ztemp                           ! temporarily storing z position
+
+
+  xts=real(part(ipart)%xlon)
+  yts=real(part(ipart)%ylat)
+  zts=real(part(ipart)%z)
+  ztseta=real(part(ipart)%zeta)
+  if (lsettling) part(ipart)%settling=0.
+
+  ! Memorize the old wind
+  !**********************
+
+  uold=u
+  vold=v
+
+  select case (wind_coord_type)
+    case ('ETA')
+      woldeta=weta
+    case ('METER')
+      wold=w
+    case default
+      wold=w
+  end select
+
+  ! Interpolate wind at new position and time
+  !******************************************
+  call interpol_wind_short(itime+part(ipart)%idt*ldirect,xts,yts,zts,ztseta)
+
+  if (mdomainfill.eq.0) then
+    if (lsettling) then
+      if ((ipin.ne.3).and.(ipin.ne.4)) then
+        do insp=1,nspec
+          nsp=insp
+          if (xmass(part(ipart)%npoint,nsp).gt.eps3) exit
+        end do
+      else
+        nsp=1
+      endif
+      if (density(nsp).gt.0.) then
+        select case (wind_coord_type)
+
+          case ('ETA')
+
+            call update_zeta_to_z(itime+part(ipart)%idt,ipart)
+            call update_z_to_zeta(itime+part(ipart)%idt,ipart)
+            zts=real(part(ipart)%z)
+            call get_settling( &
+             itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling) !bugfix
+            call w_to_weta( &
+              itime+part(ipart)%idt, real(part(ipart)%idt), part(ipart)%xlon, &
+              part(ipart)%ylat, part(ipart)%z, part(ipart)%zeta, &
+              part(ipart)%settling, weta_settling)
+            weta=weta+weta_settling
+           !woldeta=
+     !real(part(ipart)%zeta-part(ipart)%zeta_prev)/real(part(ipart)%idt*ldirect)
+
+          case ('METER')
+
+            call get_settling( &
+              itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
+            w=w+part(ipart)%settling
+
+          case default
+
+            call get_settling( &
+              itime+part(ipart)%idt,xts,yts,zts,nsp,part(ipart)%settling)
+            w=w+part(ipart)%settling
+
+        end select
+      end if
+    endif
+  end if
+
+  ! Determine the difference vector between new and old wind
+  ! (use half of it to correct position according to Petterssen)
+  !*************************************************************
+
+  u=(u-uold)*0.5
+  v=(v-vold)*0.5
+
+  select case (wind_coord_type)
+
+    case ('ETA')
+
+      weta=(weta-woldeta)/2.
+      call update_zeta(ipart,weta*real(part(ipart)%idt*ldirect))
+      if (part(ipart)%zeta.ge.1.) call set_zeta(ipart,1.-(part(ipart)%zeta-1.))
+      if (part(ipart)%zeta.eq.1.) call update_zeta(ipart,-eps_eta)
+
+    case ('METER')
+
+      w=(w-wold)/2.
+      call update_z(ipart,w*real(part(ipart)%idt*ldirect))
+      if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))          ! if particle below ground -> reflection
+
+    case default
+
+      w=(w-wold)/2.
+      call update_z(ipart,w*real(part(ipart)%idt*ldirect))
+      if (part(ipart)%z.lt.0.) call set_z(ipart,min(h-eps2,-1.*part(ipart)%z))
+
+  end select
+
+  ! Finally, correct the old position
+  !**********************************
+
+  call update_xy(u*part(ipart)%idt,v*part(ipart)%idt,ipart)
+
+  ! If particle above highest model level, set it back into the domain
+  !*******************************************************************
+
+  call pushpartdown(ipart)
+
+end subroutine petterson_corr

src/adv_pushpartdown.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine pushpartdown(ipart)
+
+  implicit none
+
+  integer, intent(in) ::          &
+    ipart                           ! particle index
+
+  eps=nxmax/3.e5
+
+  select case (wind_coord_type)
+
+    case ('ETA')
+
+      if (part(ipart)%zeta.le.real(uvheight(nz),kind=dp)) &
+        call set_zeta(ipart,uvheight(nz)+eps_eta)
+
+    case ('METER')
+
+      if (part(ipart)%z.ge.real(height(nz),kind=dp)) &
+        call set_z(ipart,height(nz)-100.*eps)
+
+    case default
+
+      if (part(ipart)%z.ge.real(height(nz),kind=dp)) &
+      call set_z(ipart,height(nz)-100.*eps)
+
+  end select
+
+end subroutine pushpartdown

src/adv_update_xy.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine update_xy(xchange,ychange,ipart)
+
+  implicit none
+
+  integer, intent(in) ::          &
+    ipart                           ! particle number
+  real, intent(in) ::             &
+    xchange,ychange                 ! change in position
+  real ::                         &
+    xlon,ylat,xpol,ypol,          & ! temporarily storing new particle positions
+    gridsize,cosfact                ! used to compute new positions of particles
+
+  eps=nxmax/3.e5
+
+  if (ngrid.ge.0) then
+
+    cosfact=dxconst/cos((real(part(ipart)%ylat)*dy+ylat0)*pi180)
+    call update_xlon(ipart,real(xchange*cosfact*real(ldirect),kind=dp))
+    call update_ylat(ipart,real(ychange*dyconst*real(ldirect),kind=dp))
+
+  else if (ngrid.eq.-1) then      ! around north pole
+
+    xlon=xlon0+real(part(ipart)%xlon)*dx !comment by MC: compute old part pos.
+    ylat=ylat0+real(part(ipart)%ylat)*dy
+    call cll2xy(northpolemap,ylat,xlon,xpol,ypol)
+      !convert old particle position in polar stereographic
+    gridsize=1000.*cgszll(northpolemap,ylat,xlon)
+      !calculate size in m of grid element in polar stereographic coordinate
+    xpol=xpol+xchange/gridsize*real(ldirect)
+      !position in grid unit polar stereographic
+    ypol=ypol+ychange/gridsize*real(ldirect)
+    call cxy2ll(northpolemap,xpol,ypol,ylat,xlon)
+      !convert to lat long coordinate
+    call set_xlon(ipart,real((xlon-xlon0)/dx,kind=dp))
+      !convert to grid units in lat long coordinate, comment by mc
+    call set_ylat(ipart,real((ylat-ylat0)/dy,kind=dp))
+
+  else if (ngrid.eq.-2) then    ! around south pole
+
+    xlon=xlon0+real(part(ipart)%xlon)*dx
+    ylat=ylat0+real(part(ipart)%ylat)*dy
+    call cll2xy(southpolemap,ylat,xlon,xpol,ypol)
+    gridsize=1000.*cgszll(southpolemap,ylat,xlon)
+    xpol=xpol+xchange/gridsize*real(ldirect)
+    ypol=ypol+ychange/gridsize*real(ldirect)
+    call cxy2ll(southpolemap,xpol,ypol,ylat,xlon)
+    call set_xlon(ipart,real((xlon-xlon0)/dx,kind=dp))
+    call set_ylat(ipart,real((ylat-ylat0)/dy,kind=dp))
+
+  endif
+
+  ! If global data are available, use cyclic boundary condition
+  !************************************************************
+
+  if (xglobal) then
+    if (part(ipart)%xlon .ge. real(nxmin1, kind=dp)) &
+      call update_xlon(ipart,-real(nxmin1, kind=dp))
+    if (part(ipart)%xlon .lt. 0.) call update_xlon(ipart,real(nxmin1, kind=dp))
+    if (part(ipart)%xlon .le. real(eps, kind=dp)) &
+      call set_xlon(ipart,real(eps, kind=dp))
+    if (abs( part(ipart)%xlon - real(nxmin1, kind=dp)) .le. eps) &
+      call set_xlon(ipart,real(nxmin1-eps,kind=dp))
+  endif
+
+  ! HSO/AL: Prevent particles from disappearing at the pole
+  !******************************************************************
+  if (sglobal .and. part(ipart)%ylat.lt.0. ) then
+    call set_xlon(ipart, &
+      mod( part(ipart)%xlon + real(nxmin1*0.5, kind=dp), real(nxmin1, kind=dp)))
+    call set_ylat(ipart,-part(ipart)%ylat)
+   ! In extremely rare cases, the ylat exceeds the bounds,
+   ! so we set it back into the domain here
+    if ( part(ipart)%ylat.gt.real(nymin1,kind=dp) ) &
+      call set_ylat(ipart, &
+        real(nymin1, kind=dp) - mod( part(ipart)%ylat, real(nymin1, kind=dp)))
+  else if (nglobal .and. part(ipart)%ylat .gt. real(nymin1, kind=dp) ) then
+    call set_xlon(ipart, &
+      mod( part(ipart)%xlon + real(nxmin1*0.5, kind=dp), real(nxmin1, kind=dp)))
+    call set_ylat(ipart,2.*real(nymin1,kind=dp) - part(ipart)%ylat)
+  endif
+
+  ! Check position: If trajectory outside model domain, terminate it
+  !*****************************************************************
+  ! Not necessary to check when using global domain, but some problems in the
+  ! meteo data could cause particles to go crazy.
+  ! if (gdomainfill) return
+
+  if (part(ipart)%xlon.lt.0. .or. part(ipart)%xlon.ge.real(nxmin1,kind=dp) &
+ .or. part(ipart)%ylat.lt.0. .or. part(ipart)%ylat.gt.real(nymin1,kind=dp)) then
+    part(ipart)%nstop=.true.
+    return
+  endif
+
+end subroutine update_xy

src/bin_concoutput_inv_nest.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine concoutput_inv_nest(itime,outnum)
+  !                        i     i
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the concentration grid and the receptor concentrations.      *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     24 May 1995                                                            *
+  !                                                                            *
+  !     13 April 1999, Major update: if output size is smaller, dump output    *
+  !                    in sparse matrix format; additional output of           *
+  !                    uncertainty                                             *
+  !                                                                            *
+  !     05 April 2000, Major update: output of age classes; output for backward*
+  !                    runs is time spent in grid cell times total mass of     *
+  !                    species.                                                *
+  !                                                                            *
+  !     17 February 2002, Appropriate dimensions for backward and forward runs *
+  !                       are now specified in file par_mod                    *
+  !                                                                            *
+  !     June 2006, write grid in sparse matrix with a single write command     *
+  !                in order to save disk space                                 *
+  !                                                                            *
+  !     2008 new sparse matrix format                                          *
+  !
+  !     January 2017,  Separate files by release but include all timesteps     *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! outnum          number of samples                                          *
+  ! ncells          number of cells with non-zero concentrations               *
+  ! sparse          .true. if in sparse matrix format, else .false.            *
+  ! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use unc_mod
+  use mean_mod
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
+  integer :: sp_count_i,sp_count_r
+  real :: sp_fact
+  real :: outnum,densityoutrecept(maxreceptor),xl,yl
+  ! RLT
+  real :: densitydryrecept(maxreceptor)
+  real :: factor_dryrecept(maxreceptor)
+
+  real(dep_prec) :: auxgrid(nclassunc)
+  real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
+  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+  real,parameter :: weightair=28.97
+  logical :: sp_zer
+  logical,save :: lnstart=.true.
+  logical,save,allocatable,dimension(:) :: lnstartrel
+  character :: adate*8,atime*6
+  character(len=3) :: anspec
+  logical :: lexist
+  character :: areldate*8,areltime*6
+
+  if(lnstart) then
+    allocate(lnstartrel(maxpointspec_act))
+    lnstartrel(:)=.true.
+  endif
+  print*, 'lnstartrel = ',lnstartrel
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+
+  print*, 'outnum:',outnum
+  print*, 'datetime:',adate//atime
+
+  ! For forward simulations, output fields have dimension MAXSPEC,
+  ! for backward simulations, output fields have dimension MAXPOINT.
+  ! Thus, make loops either about nspec, or about numpoint
+  !*****************************************************************
+
+
+    if (ldirect.eq.1) then
+       do ks=1,nspec
+         do kp=1,maxpointspec_act
+           tot_mu(ks,kp)=1
+         end do
+       end do
+   else
+      do ks=1,nspec
+             do kp=1,maxpointspec_act
+               tot_mu(ks,kp)=xmass(kp,ks)
+             end do
+      end do
+    endif
+
+
+  !*******************************************************************
+  ! Compute air density: sufficiently accurate to take it
+  ! from coarse grid at some time
+  ! Determine center altitude of output layer, and interpolate density
+  ! data to that altitude
+  !*******************************************************************
+
+  do kz=1,numzgrid
+    if (kz.eq.1) then
+      halfheight=outheight(1)/2.
+    else
+      halfheight=(outheight(kz)+outheight(kz-1))/2.
+    endif
+    do kzz=2,nz
+      if ((height(kzz-1).lt.halfheight).and. &
+           (height(kzz).gt.halfheight)) goto 46
+    end do
+46   kzz=max(min(kzz,nz),2)
+    dz1=halfheight-height(kzz-1)
+    dz2=height(kzz)-halfheight
+    dz=dz1+dz2
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        xl=outlon0n+real(ix)*dxoutn
+        yl=outlat0n+real(jy)*dyoutn
+        xl=(xl-xlon0)/dx
+        yl=(yl-ylat0)/dy
+        iix=max(min(nint(xl),nxmin1),0)
+        jjy=max(min(nint(yl),nymin1),0)
+        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
+             rho(iix,jjy,kzz-1,2)*dz2)/dz
+  ! RLT
+        densitydrygrid(ix,jy,kz)=(rho_dry(iix,jjy,kzz,2)*dz1+ &
+             rho_dry(iix,jjy,kzz-1,2)*dz2)/dz
+      end do
+    end do
+  end do
+
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    densityoutrecept(i)=rho(iix,jjy,1,2)
+  ! RLT
+    densitydryrecept(i)=rho_dry(iix,jjy,1,2)
+  end do
+
+  ! RLT
+  ! conversion factor for output relative to dry air
+  factor_drygrid=densityoutgrid/densitydrygrid
+  factor_dryrecept=densityoutrecept/densitydryrecept
+
+  ! Output is different for forward and backward simulations
+    do kz=1,numzgrid
+      do jy=0,numygridn-1
+        do ix=0,numxgridn-1
+          if (ldirect.eq.1) then
+            factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
+          else
+            factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
+          endif
+        end do
+      end do
+    end do
+
+  !*********************************************************************
+  ! Determine the standard deviation of the mean concentration or mixing
+  ! ratio (uncertainty of the output) and the dry and wet deposition
+  !*********************************************************************
+
+  do ks=1,nspec
+
+  write(anspec,'(i3.3)') ks
+
+    do kp=1,maxpointspec_act
+
+      print*, 'itime = ',itime
+      print*, 'lage(1) = ',lage(1)
+      print*, 'ireleasestart(kp) = ',ireleasestart(kp)
+      print*, 'ireleaseend(kp) = ',ireleaseend(kp)
+
+      ! check itime is within release and backward trajectory length
+      if (nageclass.eq.1) then
+        if ((itime.gt.ireleaseend(kp)).or.(itime.lt.(ireleasestart(kp)-lage(1)))) then
+          go to 10
+        endif
+      endif
+
+      ! calculate date of release
+      jul=bdate+real(ireleasestart(kp),kind=dp)/86400._dp    ! this is the current day
+      call caldate(jul,jjjjmmdd,ihmmss)
+      write(areldate,'(i8.8)') jjjjmmdd
+      write(areltime,'(i6.6)') ihmmss
+      print*, areldate//areltime
+
+      ! calculate date of field
+      jul=bdate+real(itime,kind=dp)/86400._dp
+      call caldate(jul,jjjjmmdd,ihmmss)
+      write(adate,'(i8.8)') jjjjmmdd
+      write(atime,'(i6.6)') ihmmss
+      print*, adate//atime
+
+      if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+        if (ldirect.eq.1) then
+          ! concentrations
+          inquire(file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
+                  areltime//'_'//anspec,exist=lexist)
+          if(lexist.and..not.lnstartrel(kp)) then
+            ! open and append to existing file
+            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
+                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
+          else
+            ! open new file
+            open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_'//areldate// &
+                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
+          endif
+        else
+          ! residence times
+          inquire(file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
+                  areltime//'_'//anspec,exist=lexist)
+          if(lexist.and..not.lnstartrel(kp)) then
+            ! open and append to existing file
+            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
+                 areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
+          else
+            ! open new file
+            open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_'//areldate// &
+                 areltime//'_'//anspec,form='unformatted',status='replace',action='write')
+          endif
+        endif
+        write(unitoutgrid) jjjjmmdd
+        write(unitoutgrid) ihmmss
+      endif
+
+      if ((iout.eq.2).or.(iout.eq.3)) then
+        ! mixing ratio
+        inquire(file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
+                areltime//'_'//anspec,exist=lexist)
+        if(lexist.and..not.lnstartrel(kp)) then
+          ! open and append to existing file
+          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
+               areltime//'_'//anspec,form='unformatted',status='old',action='write',access='append')
+        else
+          ! open new file
+          open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_'//areldate// &
+               areltime//'_'//anspec,form='unformatted',status='replace',action='write')
+        endif
+        write(unitoutgridppt) jjjjmmdd
+        write(unitoutgridppt) ihmmss
+      endif
+
+      lnstartrel(kp)=.false.
+
+      do nage=1,nageclass
+
+        do jy=0,numygridn-1
+          do ix=0,numxgridn-1
+
+  ! CONCENTRATION OR MIXING RATIO
+            do kz=1,numzgrid
+              do l=1,nclassunc
+                auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,grid(ix,jy,kz), &
+                     gridsigma(ix,jy,kz),nclassunc)
+  ! Multiply by number of classes to get total concentration
+              grid(ix,jy,kz)= &
+                   grid(ix,jy,kz)*nclassunc
+  ! Calculate standard deviation of the mean
+              gridsigma(ix,jy,kz)= &
+                   gridsigma(ix,jy,kz)* &
+                   sqrt(real(nclassunc))
+            end do
+          end do
+        end do
+
+
+  !*******************************************************************
+  ! Generate output: may be in concentration (ng/m3) or in mixing
+  ! ratio (ppt) or both
+  ! Output the position and the values alternated multiplied by
+  ! 1 or -1, first line is number of values, number of positions
+  ! For backward simulations, the unit is seconds, stored in grid_time
+  !*******************************************************************
+
+  ! Concentration output
+  !*********************
+
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+
+  ! Concentrations
+
+  ! surf_only write only 1st layer
+
+         sp_count_i=0
+         sp_count_r=0
+         sp_fact=-1.
+         sp_zer=.true.
+          do kz=1,1
+            do jy=0,numygridn-1
+              do ix=0,numxgridn-1
+                if (grid(ix,jy,kz).gt.smallnum) then
+                  if (sp_zer.eqv..true.) then ! first non zero value
+                    sp_count_i=sp_count_i+1
+                    sparse_dump_i(sp_count_i)= &
+                         ix+jy*numxgridn+kz*numxgridn*numygridn
+                    sp_zer=.false.
+                    sp_fact=sp_fact*(-1.)
+                   endif
+                   sp_count_r=sp_count_r+1
+                   sparse_dump_r(sp_count_r)= &
+                        sp_fact* &
+                        grid(ix,jy,kz)* &
+                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
+  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
+  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
+                   sparse_dump_u(sp_count_r)= &
+                        gridsigma(ix,jy,kz)* &
+                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
+              else ! concentration is zero
+                  sp_zer=.true.
+              endif
+              end do
+            end do
+          end do
+         write(unitoutgrid) sp_count_i
+         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
+         write(unitoutgrid) sp_count_r
+         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
+
+      endif !  concentration output
+
+  ! Mixing ratio output
+  !********************
+
+      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
+
+
+  ! Mixing ratios
+
+    ! surf_only write only 1st layer
+
+         sp_count_i=0
+         sp_count_r=0
+         sp_fact=-1.
+         sp_zer=.true.
+          do kz=1,1
+            do jy=0,numygridn-1
+              do ix=0,numxgridn-1
+                if (grid(ix,jy,kz).gt.smallnum) then
+                  if (sp_zer.eqv..true.) then ! first non zero value
+                    sp_count_i=sp_count_i+1
+                    sparse_dump_i(sp_count_i)= &
+                         ix+jy*numxgridn+kz*numxgridn*numygridn
+                    sp_zer=.false.
+                    sp_fact=sp_fact*(-1.)
+                 endif
+                 sp_count_r=sp_count_r+1
+                 sparse_dump_r(sp_count_r)= &
+                      sp_fact* &
+                      1.e12*grid(ix,jy,kz) &
+                      /volumen(ix,jy,kz)/outnum* &
+                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
+                 sparse_dump_u(sp_count_r)= &
+                      1.e12*gridsigma(ix,jy,kz)/volumen(ix,jy,kz)/ &
+                      outnum*weightair/weightmolar(ks)/ &
+                      densityoutgrid(ix,jy,kz)
+              else ! concentration is zero
+                  sp_zer=.true.
+              endif
+              end do
+            end do
+          end do
+          write(unitoutgridppt) sp_count_i
+          write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
+          write(unitoutgridppt) sp_count_r
+          write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
+
+        endif ! output for ppt
+
+      end do ! nageclass
+
+      close(unitoutgridppt)
+      close(unitoutgrid)
+
+      ! itime is outside range
+10    continue
+
+    end do ! maxpointspec_act
+
+  end do ! nspec
+
+
+  ! RLT Aug 2017
+  ! Write out conversion factor for dry air
+  inquire(file=path(2)(1:length(2))//'factor_drygrid_nest',exist=lexist)
+  if (lexist.and..not.lnstart) then
+    ! open and append
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='old',action='write',access='append')
+  else
+    ! create new
+    open(unitoutfactor,file=path(2)(1:length(2))//'factor_drygrid_nest',form='unformatted',&
+            status='replace',action='write')
+  endif
+  sp_count_i=0
+  sp_count_r=0
+  sp_fact=-1.
+  sp_zer=.true.
+  do kz=1,1
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        if (factor_drygrid(ix,jy,kz).gt.(1.+smallnum).or.factor_drygrid(ix,jy,kz).lt.(1.-smallnum)) then
+          if (sp_zer.eqv..true.) then ! first value not equal to one
+            sp_count_i=sp_count_i+1
+            sparse_dump_i(sp_count_i)= &
+                  ix+jy*numxgridn+kz*numxgridn*numygridn
+            sp_zer=.false.
+            sp_fact=sp_fact*(-1.)
+          endif
+          sp_count_r=sp_count_r+1
+          sparse_dump_r(sp_count_r)= &
+               sp_fact*factor_drygrid(ix,jy,kz)
+        else ! factor is one
+          sp_zer=.true.
+        endif
+      end do
+    end do
+  end do
+  write(unitoutfactor) sp_count_i
+  write(unitoutfactor) (sparse_dump_i(i),i=1,sp_count_i)
+  write(unitoutfactor) sp_count_r
+  write(unitoutfactor) (sparse_dump_r(i),i=1,sp_count_r)
+  close(unitoutfactor)
+
+  ! reset lnstart
+  if (lnstart) then
+    lnstart=.false.
+  endif
+
+  ! Reinitialization of grid
+  !*************************
+
+  do ks=1,nspec
+  do kp=1,maxpointspec_act
+    do i=1,numreceptor
+      creceptor(i,ks)=0.
+    end do
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        do l=1,nclassunc
+          do nage=1,nageclass
+            do kz=1,numzgrid
+              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
+            end do
+          end do
+        end do
+      end do
+    end do
+  end do
+  end do
+end subroutine concoutput_inv_nest

src/bin_initcond_output.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine initcond_output(itime)
+  !                                 i
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the initial condition sensitivity field.                     *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     24 May 1995                                                            *
+  !                                                                            *
+  !     13 April 1999, Major update: if output size is smaller, dump output    *
+  !                    in sparse matrix format; additional output of           *
+  !                    uncertainty                                             *
+  !                                                                            *
+  !     05 April 2000, Major update: output of age classes; output for backward*
+  !                    runs is time spent in grid cell times total mass of     *
+  !                    species.                                                *
+  !                                                                            *
+  !     17 February 2002, Appropriate dimensions for backward and forward runs *
+  !                       are now specified in file par_mod                    *
+  !                                                                            *
+  !     June 2006, write grid in sparse matrix with a single write command     *
+  !                in order to save disk space                                 *
+  !                                                                            *
+  !     2008 new sparse matrix format                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! ncells          number of cells with non-zero concentrations               *
+  ! sparse          .true. if in sparse matrix format, else .false.            *
+  !                                                                            *
+  !*****************************************************************************
+
+  use unc_mod
+
+  implicit none
+
+  integer :: itime,i,ix,jy,kz,ks,kp,sp_count_i,sp_count_r
+  real :: sp_fact,fact_recept
+  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
+  logical :: sp_zer
+  character(len=3) :: anspec
+
+
+  !*********************************************************************
+  ! Determine the standard deviation of the mean concentration or mixing
+  ! ratio (uncertainty of the output) and the dry and wet deposition
+  !*********************************************************************
+
+  do ks=1,nspec
+
+    write(anspec,'(i3.3)') ks
+    open(97,file=path(2)(1:length(2))//'grid_initial'// &
+         '_'//anspec,form='unformatted')
+    write(97) itime
+
+    do kp=1,maxpointspec_act
+
+      if (ind_rel.eq.1) then
+        fact_recept=rho_rel(kp)
+      else
+        fact_recept=1.
+      endif
+
+  !*******************************************************************
+  ! Generate output: may be in concentration (ng/m3) or in mixing
+  ! ratio (ppt) or both
+  ! Output the position and the values alternated multiplied by
+  ! 1 or -1, first line is number of values, number of positions
+  ! For backward simulations, the unit is seconds, stored in grid_time
+  !*******************************************************************
+
+  ! Write out dummy "wet and dry deposition" fields, to keep same format
+  ! as for concentration output
+      sp_count_i=0
+      sp_count_r=0
+      write(97) sp_count_i
+      write(97) (sparse_dump_i(i),i=1,sp_count_i)
+      write(97) sp_count_r
+      write(97) (sparse_dump_r(i),i=1,sp_count_r)
+      write(97) sp_count_i
+      write(97) (sparse_dump_i(i),i=1,sp_count_i)
+      write(97) sp_count_r
+      write(97) (sparse_dump_r(i),i=1,sp_count_r)
+
+
+  ! Write out sensitivity to initial conditions
+      sp_count_i=0
+      sp_count_r=0
+      sp_fact=-1.
+      sp_zer=.true.
+      do kz=1,numzgrid
+        do jy=0,numygrid-1
+          do ix=0,numxgrid-1
+            if (init_cond(ix,jy,kz,ks,kp).gt.smallnum) then
+              if (sp_zer.eqv..true.) then ! first non zero value
+                sp_count_i=sp_count_i+1
+                sparse_dump_i(sp_count_i)= &
+                     ix+jy*numxgrid+kz*numxgrid*numygrid
+                sp_zer=.false.
+                sp_fact=sp_fact*(-1.)
+              endif
+              sp_count_r=sp_count_r+1
+              sparse_dump_r(sp_count_r)=sp_fact* &
+                   init_cond(ix,jy,kz,ks,kp)/xmass(kp,ks)*fact_recept
+            else ! concentration is zero
+              sp_zer=.true.
+            endif
+          end do
+        end do
+      end do
+      write(97) sp_count_i
+      write(97) (sparse_dump_i(i),i=1,sp_count_i)
+      write(97) sp_count_r
+      write(97) (sparse_dump_r(i),i=1,sp_count_r)
+
+
+    end do
+
+    close(97)
+
+  end do
+end subroutine initcond_output

src/bin_openreceptors.f90

+! SPDX-FileCopyrightText: FLEXPART 1998-2023, see flexpart_license.txt
+! SPDX-License-Identifier: GPL-3.0-or-later
+
+subroutine openreceptors
+
+  !*****************************************************************************
+  !                                                                            *
+  !  This routine opens the receptor output files and writes out the receptor  *
+  !  names and the receptor locations. The receptor output files are not       *
+  !  closed, but kept open throughout the simulation. Concentrations are       *
+  !  continuously  dumped to these files.                                      *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     7 August 2002                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! numreceptor            actual number of receptor points specified          *
+  ! receptornames          names of the receptor points                        *
+  ! xreceptor,yreceptor    coordinates of the receptor points                  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: j
+
+  ! Open output file for receptor points and write out a short header
+  ! containing receptor names and locations
+  !******************************************************************
+
+  if (numreceptor.ge.1) then           ! do it only if receptors are specified
+
+  ! Concentration output
+  !*********************
+
+    if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+      open(unitoutrecept,file=path(2)(1:length(2))//'receptor_conc', &
+           form='unformatted',err=997)
+      write(unitoutrecept) (receptorname(j),j=1,numreceptor)
+      write(unitoutrecept) (xreceptor(j)*dx+xlon0, &
+           yreceptor(j)*dy+ylat0,j=1,numreceptor)
+    endif
+
+  ! Mixing ratio output
+  !********************
+
+    if ((iout.eq.2).or.(iout.eq.3)) then
+      open(unitoutreceptppt,file=path(2)(1:length(2))//'receptor_pptv', &
+           form='unformatted',err=998)
+      write(unitoutreceptppt) (receptorname(j),j=1,numreceptor)
+      write(unitoutreceptppt) (xreceptor(j)*dx+xlon0, &
+           yreceptor(j)*dy+ylat0,j=1,numreceptor)
+    endif
+  endif
+
+  return
+997   write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE           #### '
+  write(*,*) ' ####              receptor_conc               #### '
+  write(*,*) ' #### CANNOT BE OPENED.                        #### '
+  stop
+
+998   write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE           #### '
+  write(*,*) ' ####              receptor_pptv               #### '
+  write(*,*) ' #### CANNOT BE OPENED.                        #### '
+  stop
+end subroutine openreceptors