! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later

  !*****************************************************************************
  !                                                                            *
  !   L. Bakels 2021: This module contains subroutines related to the          *
  !                   initialisation of the particles                          *
  !                                                                            *
  !*****************************************************************************

module initialise_mod

  use com_mod
  use par_mod
  use date_mod
  use particle_mod
  use windfields_mod
  use random_mod
#ifdef ETA
  use coord_ecmwf_mod
#endif
#ifdef USE_NCF
  use netcdf_output_mod
#endif

  implicit none

  !**********************************************************
  ! Variables used for domain-filling trajectory calculations
  !**********************************************************

  integer ::    &
    nx_we(2),   & ! x indices of western and eastern boundary of domain-filling.
    ny_sn(2),   & ! y indices of southern and northern boundary of domain-filling.
    numcolumn     ! Maximum number of particles to be released within a single column.
  integer,allocatable,dimension(:,:) :: &
    numcolumn_we,   & ! Number of particles to be released within one column
                      ! at the western and eastern boundary surfaces.
    numcolumn_sn      ! Same as numcolumn_we, but for southern and northern domain boundary.
  real,allocatable,dimension(:,:,:) ::  &
    zcolumn_we,     & ! Altitudes where particles are to be released
                      ! at the western and eastern boundary surfaces.
    zcolumn_sn,     & ! Same as zcolumn_we, but for southern and northern domain boundary.
    acc_mass_we,    & ! Mass that has accumulated at the western and eastern boundary;
                      ! if it exceeds xmassperparticle, a particle is released and
                      ! acc_mass_we is reduced accordingly.
    acc_mass_sn       ! Same as acc_mass_we, but for southern and northern domain boundary
  real ::           &
    xmassperparticle  ! Air mass per particle in the domain-filling traj. option.

contains

subroutine alloc_domainfill
  implicit none
  allocate(numcolumn_we(2,0:nymax-1),numcolumn_sn(2,0:nxmax-1))
  allocate(zcolumn_we(2,0:nymax-1,maxcolumn),zcolumn_sn(2,0:nxmax-1,maxcolumn), &
    acc_mass_we(2,0:nymax-1,maxcolumn),acc_mass_sn(2,0:nxmax-1,maxcolumn))              
end subroutine alloc_domainfill

subroutine dealloc_domainfill
  if (mdomainfill.lt.1) return
  deallocate(numcolumn_we,numcolumn_sn,zcolumn_sn,zcolumn_we,acc_mass_sn, &
    acc_mass_we)
end subroutine dealloc_domainfill

subroutine releaseparticles(itime)
  !                              o
  !*****************************************************************************
  !                                                                            *
  !     This subroutine releases particles from the release locations.         *
  !                                                                            *
  !     It searches for a "vacant" storage space and assigns all particle      *
  !     information to that space. A space is vacant either when no particle   *
  !     is yet assigned to it, or when it's particle is expired and, thus,     *
  !     the storage space is made available to a new particle.                 *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     29 June 2002                                                           *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! itime [s]            current time                                          *
  ! ireleasestart, ireleaseend          start and end times of all releases    *
  ! npart(maxpoint)      number of particles to be released in total           *
  ! numrel               number of particles to be released during this time   *
  !                      step                                                  *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use xmass_mod
  use output_mod
  use interpol_mod

  implicit none

  !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)
  real :: xaux,yaux,zaux,rfraction
  real :: xlonav,timecorrect(maxspec)
  real :: average_timecorrect
  integer :: itime,numrel,i,j,k,ipart,minpart
  integer :: istart,iend,totpart,iterm_index
  integer :: nweeks,ndayofweek,nhour,jjjjmmdd,ihmmss,mm
  real(kind=dp) :: julmonday,jul,jullocal,juldiff

  integer :: idummy = -7
  !save idummy,xmasssave
  !data idummy/-7/,xmasssave/maxpoint*0./

  real :: eps
  eps=nxmax/3.e5


  ! Determine the actual date and time in Greenwich 
  ! (i.e., UTC + correction for daylight savings time)
  !***************************************************

  julmonday=juldate(19000101,0)          ! this is a Monday
  jul=bdate+real(itime,kind=dp)/86400._dp    ! this is the current day
  call caldate(jul,jjjjmmdd,ihmmss)
  mm=(jjjjmmdd-10000*(jjjjmmdd/10000))/100
  if ((mm.ge.4).and.(mm.le.9)) jul=jul+1._dp/24._dp   ! daylight savings time


  ! For every release point, check whether we are in the release time interval
  !***************************************************************************
  ! First allocate all particles that are going to be in the simulation
  ! If ipin==0,1, and ipout==0, then dead particles can be overwritten to save memory
  if (ipin.gt.1 .or. ipout.ne.0) then
    if (itime.eq.0) then
      totpart=0
      do i=1,numpoint
        totpart = totpart+npart(i)
      end do
      call alloc_particles(totpart)
    else if (itime.eq.itime_init) then !From restart point only allocate particles that are yet to be born
      totpart=0
      do i=1,numpoint
        totpart = totpart+npart(i)
      end do
      if (totpart.gt.count%allocated) call alloc_particles(totpart-count%allocated)
    end if 
  endif

  call get_totalpart_num(istart)
  if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
  minpart=1
  do i=1,numpoint
    if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
         (itime.le.ireleaseend(i))) then

  ! Determine the local day and time
  !*********************************

      xlonav=xlon0+(xpoint2(i)+xpoint1(i))*0.5*dx  ! longitude needed to determine local time
      if (xlonav.lt.-180.) xlonav=xlonav+360.
      if (xlonav.gt.180.) xlonav=xlonav-360.
      jullocal=jul+real(xlonav,kind=dp)/360._dp   ! correct approximately for time zone to obtain local time

      juldiff=jullocal-julmonday
      nweeks=int(juldiff/7._dp)
      juldiff=juldiff-real(nweeks,kind=dp)*7._dp
      ndayofweek=int(juldiff)+1              ! this is the current day of week, starting with Monday
      nhour=nint((juldiff-real(ndayofweek-1,kind=dp))*24._dp)    ! this is the current hour
      if (nhour.eq.0) then
        nhour=24
        ndayofweek=ndayofweek-1
        if (ndayofweek.eq.0) ndayofweek=7
      endif

  ! Calculate a species- and time-dependent correction factor, distinguishing between
  ! area (those with release starting at surface) and point (release starting above surface) sources
  ! Also, calculate an average time correction factor (species independent)
  !*****************************************************************************
      average_timecorrect=0.
      do k=1,nspec
        if(abs(xpoint2(i)-xpoint1(i)).lt.1.E-4.and.abs(ypoint2(i)-ypoint1(i)).lt.1.E-4) then
  !        if (zpoint1(i).gt.0.5) then      ! point source
          timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
        else                             ! area source
          timecorrect(k)=area_hour(k,nhour)*area_dow(k,ndayofweek)
        endif
        average_timecorrect=average_timecorrect+timecorrect(k)
      end do
      average_timecorrect=average_timecorrect/real(nspec)

  ! Determine number of particles to be released this time; at start and at end of release,
  ! only half the particles are released
  !*****************************************************************************

      if (ireleasestart(i).ne.ireleaseend(i)) then
        rfraction=abs(real(npart(i))*real(lsynctime)/ &
             real(ireleaseend(i)-ireleasestart(i)))
        if ((itime.eq.ireleasestart(i)).or. &
             (itime.eq.ireleaseend(i))) rfraction=rfraction*0.5

  ! Take the species-average time correction factor in order to scale the
  ! number of particles released this time
  !**********************************************************************
        rfraction=rfraction*average_timecorrect

        rfraction=rfraction+xmasssave(i)  ! number to be released at this time
        numrel=int(rfraction)
        xmasssave(i)=rfraction-real(numrel)
      else
        numrel=npart(i)
      endif

      xaux=xpoint2(i)-xpoint1(i)
      yaux=ypoint2(i)-ypoint1(i)
      zaux=zpoint2(i)-zpoint1(i)

      if (ipin.le.1 .and. ipout.eq.0) then
        call rewrite_iterm()
        totpart = numrel-count%iterm_max
        if (totpart.gt.0) call alloc_particles(totpart)
        call rewrite_iterm()
        iterm_index=1
      endif
      do j=1,numrel             ! loop over particles to be released this time
        call get_newpart_index(ipart,iterm_index)
        call spawn_particle(itime, ipart)

  ! Particle coordinates are determined by using a random position within the release volume
  !*****************************************************************************

  ! Determine horizontal particle position
  !***************************************
        call set_xlon(ipart,real(xpoint1(i)+ran1(idummy,0)*xaux,kind=dp))
        if (xglobal) then
          if (part(ipart)%xlon.gt.real(nxmin1,kind=dp)) &
            call set_xlon(ipart,-real(nxmin1,kind=dp))
          if (part(ipart)%xlon.lt.0.) &
            call set_xlon(ipart,real(nxmin1,kind=dp))
        endif
        call set_ylat(ipart,real(ypoint1(i)+ran1(idummy,0)*yaux,kind=dp))

  ! Assign mass to particle: Total mass divided by total number of particles.
  ! Time variation has partly been taken into account already by a species-average
  ! correction factor, by which the number of particles released this time has been
  ! scaled. Adjust the mass per particle by the species-dependent time correction factor
  ! divided by the species-average one
  ! for the scavenging calculation the mass needs to be multiplied with rho of the particle layer and
  ! divided by the sum of rho of all particles.
  !*****************************************************************************
        do k=1,nspec
          mass(ipart,k)=xmass(i,k)/real(npart(i)) &
                *timecorrect(k)/average_timecorrect
          mass_init(ipart,k)=mass(ipart,k)
        end do
  ! Assign certain properties to particle
  !**************************************
        part(ipart)%nclass=min(int(ran1(idummy,0)*real(nclassunc))+1, &
             nclassunc)
        numparticlecount=numparticlecount+1
        if (mquasilag.eq.0) then
          part(ipart)%npoint=i
        else
          part(ipart)%npoint=numparticlecount
        endif
        part(ipart)%idt=mintime               ! first time step

        ! Determine vertical particle position
        !*************************************
        call set_z(ipart,zpoint1(i)+ran1(idummy,0)*zaux)
        ! Interpolation of topography and density
        !****************************************

        ! Transform the verticle particle position from pressure or sea level to above ground
        ! if necessary
        !************************************************************************************
        call kindz_to_z(ipart)
#ifdef ETA
        call z_to_zeta(itime,part(ipart)%xlon,part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta)
        part(ipart)%etaupdate = .true. ! The z(meter) coordinate is up to date
        part(ipart)%meterupdate = .true.
#endif
        
        call init_mass_conversion(ipart,i)

        call get_totalpart_num(numpart)

      end do  ! numrel 
    endif ! releasepoint
  end do ! numpoint

  if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()

  call get_totalpart_num(iend)

  ! NetCDF only: write initial positions of new particles
! #ifdef USE_NCF
!   if ((iend-istart.gt.0).and.(ipout.ge.1)) then 
!     call wrt_part_initialpos(itime,istart,iend)
!     call output_particles(itime,.true.)
!   endif
! #endif
  return

! 996   continue
!   write(*,*) '#####################################################'
!   write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
!   write(*,*) '####                                             ####'
!   write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED  ####'
!   write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE  ####'
!   write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
!   write(*,*) '#### OR REDUCE NUMBER OF RELEASE POINTS.         ####'
!   write(*,*) '#####################################################'
!   stop

end subroutine releaseparticles

subroutine kindz_to_z(ipart)
  use point_mod
  use xmass_mod
  use output_mod
  use interpol_mod

  implicit none

  integer,intent(in) :: ipart
  integer :: kz
  real :: dp1,dp2,press,presspart,pressold,topo,r,t
  real,parameter :: eps2=1.e-6


  ! Determine the nest we are in
  !*****************************
  call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)

  ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
  !*****************************************************************************
  call find_grid_indices(real(part(ipart)%xlon),real(part(ipart)%ylat))
  call find_grid_distances(real(part(ipart)%xlon),real(part(ipart)%ylat))


  if (kindz(part(ipart)%npoint).eq.1) return ! Nothing needs to happen

  ! If starting height is in pressure coordinates, retrieve pressure profile and 
  ! convert zpart1 to meters
  !*****************************************************************************

  if (kindz(part(ipart)%npoint).eq.3) then
    presspart=real(part(ipart)%z)
    do kz=1,nz
      if (ngrid.gt.0) then
        r=p1*rhon(ix ,jy ,kz,2,ngrid) &
             +p2*rhon(ixp,jy ,kz,2,ngrid) &
             +p3*rhon(ix ,jyp,kz,2,ngrid) &
             +p4*rhon(ixp,jyp,kz,2,ngrid)
        t=p1*ttn(ix ,jy ,kz,2,ngrid) &
             +p2*ttn(ixp,jy ,kz,2,ngrid) &
             +p3*ttn(ix ,jyp,kz,2,ngrid) &
             +p4*ttn(ixp,jyp,kz,2,ngrid)
      else
        r=p1*rho(ix ,jy ,kz,2) &
             +p2*rho(ixp,jy ,kz,2) &
             +p3*rho(ix ,jyp,kz,2) &
             +p4*rho(ixp,jyp,kz,2)
        t=p1*tt(ix ,jy ,kz,2) &
             +p2*tt(ixp,jy ,kz,2) &
             +p3*tt(ix ,jyp,kz,2) &
             +p4*tt(ixp,jyp,kz,2)
      endif
      press=r*r_air*t/100.
      if (kz.eq.1) pressold=press

      if (press.lt.presspart) then
        if (kz.eq.1) then
          call set_z(ipart,height(1)*0.5)
        else
          dp1=pressold-presspart
          dp2=presspart-press
          call set_z(ipart,(height(kz-1)*dp2+height(kz)*dp1) &
               /(dp1+dp2))
        endif
        exit
      endif
      pressold=press
    end do

  ! If release positions are given in meters above sea level, subtract the
  ! topography from the starting height
  !***********************************************************************

  else if (kindz(part(ipart)%npoint).eq.2) then
    if (ngrid.gt.0) then
      topo=p1*oron(ix ,jy ,ngrid) &
           + p2*oron(ixp,jy ,ngrid) &
           + p3*oron(ix ,jyp,ngrid) &
           + p4*oron(ixp,jyp,ngrid)
    else
      topo=p1*oro(ix ,jy) &
           + p2*oro(ixp,jy) &
           + p3*oro(ix ,jyp) &
           + p4*oro(ixp,jyp)
    endif
    call update_z(ipart,-topo)
  endif
  if (part(ipart)%z.lt.eps2) call set_z(ipart,eps2)   ! Minimum starting height is eps2
  if (part(ipart)%z.gt.height(nz)-0.5) &
    call set_z(ipart,height(nz)-0.5) ! Maximum starting height is uppermost level - 0.5 meters

  if (ipin.eq.3 .or. ipin.eq.4) then
    if (part(ipart)%z.gt.zpoint2(part(ipart)%npoint)) &
      zpoint2(part(ipart)%npoint)=real(part(ipart)%z)
    if (part(ipart)%z.lt.zpoint1(part(ipart)%npoint)) &
      zpoint1(part(ipart)%npoint)=real(part(ipart)%z)
  endif    

end subroutine kindz_to_z

subroutine init_mass_conversion(ipart,ipoint)
  ! For special simulations, multiply particle concentration air density;
  ! Simply take the 2nd field in memory to do this (accurate enough)
  !***********************************************************************
  !AF IND_SOURCE switches between different units for concentrations at the source
  !Af    NOTE that in backward simulations the release of particles takes place at the
  !Af         receptor and the sampling at the source.
  !Af          1="mass"
  !Af          2="mass mixing ratio"
  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
  !Af          1="mass"
  !Af          2="mass mixing ratio"
  !            3 = wet deposition in outputfield
  !            4 = dry deposition in outputfield

  !Af switches for the releasefile:
  !Af IND_REL =  1 : xmass * rho
  !Af IND_REL =  0 : xmass * 1

  !Af ind_rel is defined in readcommand.f
  use point_mod
  use xmass_mod
  use interpol_mod

  implicit none

  integer,intent(in) :: ipart,ipoint
  integer :: n
  real :: rhoaux(2),rhoout
  real :: dz1,dz2


  if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then

    ! Interpolate the air density, horizontal values are computed in kindz_to_z
    !**************************************************************************
    call find_z_level_meters(real(part(ipart)%z))
    call find_vert_vars_lin(height,real(part(ipart)%z),indz,dz1,dz2,lbounds)

    if (ngrid.gt.0) then
      do n=1,2
        rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,2,ngrid) &
             +p2*rhon(ixp,jy ,indz+n-1,2,ngrid) &
             +p3*rhon(ix ,jyp,indz+n-1,2,ngrid) &
             +p4*rhon(ixp,jyp,indz+n-1,2,ngrid)
      end do
    else
      do n=1,2
        rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,2) &
             +p2*rho(ixp,jy ,indz+n-1,2) &
             +p3*rho(ix ,jyp,indz+n-1,2) &
             +p4*rho(ixp,jyp,indz+n-1,2)
      end do
    endif
    rhoout=dz2*rhoaux(1)+dz1*rhoaux(2)
    rho_rel(ipoint)=rhoout


    ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density
    !***********************************************************************
    mass(ipart,:)=mass(ipart,:)*rhoout
    mass_init(ipart,:)=mass(ipart,:)
  endif
end subroutine init_mass_conversion

subroutine readpartpositions

  !*****************************************************************************
  !                                                                            *
  !   This routine opens the particle dump file and reads all the particle     *
  !   positions from a previous run to initialize the current run.             *
  !                                                                            *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     24 March 2000                                                          *
  !                                                                            *
  !  Changes                                                                   *
  !     2022, L. Bakels: NetCDF option for reading particle information        *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  !*****************************************************************************
  implicit none

  integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,lix,ios
  integer :: id1,id2,it1,it2
  real :: xlonin,ylatin,topo,hmixi,pvi,qvi,rhoi,tri,tti
  character :: specin*7
  real(kind=dp) :: julin,julpartin

  integer :: idummy = -8

  numparticlecount=0

  ! Open header file of dumped particle data
  !*****************************************
#ifdef USE_NCF
    call readpartpositions_netcdf(ibtime,ibdate)
    call get_totalpart_num(numpart)
    numparticlecount=numpart
    return
#endif

  open(unitpartin,file=path(2)(1:length(2))//'header', &
       form='unformatted',err=998)

  read(unitpartin) ibdatein,ibtimein
  read(unitpartin)
  read(unitpartin)

  read(unitpartin)
  read(unitpartin)
  read(unitpartin) nspecin
  nspecin=nspecin/3
  if ((ldirect.eq.1).and.(nspec.ne.nspecin)) then
    write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
    write(*,*) ' #### THE NUMBER OF SPECIES TO BE READ IN DOES #### '
    write(*,*) ' #### NOT AGREE WITH CURRENT SETTINGS!         #### '
    stop
  end if

  do i=1,nspecin
    read(unitpartin)
    read(unitpartin)
    read(unitpartin) j,specin
    if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) then
      write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
      write(*,*) ' #### SPECIES NAMES TO BE READ IN DO NOT       #### '
      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
      stop
    end if
  end do

  read(unitpartin) numpointin
  if (numpointin.ne.numpoint) then
    write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
    write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
    write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
  end if 
  do i=1,numpointin
    read(unitpartin)
    read(unitpartin)
    read(unitpartin)
    read(unitpartin)
    do j=1,nspec
      read(unitpartin)
      read(unitpartin)
      read(unitpartin)
    end do
  end do
  read(unitpartin)
  read(unitpartin)

  do lix=0,numxgrid-1
    read(unitpartin)
  end do


  ! Open data file of dumped particle data
  !***************************************

  close(unitpartin)
  open(unitpartin,file=path(2)(1:length(2))//'partposit_end', &
       form='unformatted',err=998)
  

  do 
    read(unitpartin,iostat=ios) itimein
    if (ios.lt.0) exit
    i=0
    do
      i=i+1
      read(unitpartin) part(i)%npoint,xlonin,ylatin,part(i)%z,part(i)%tstart, &
           topo,pvi,qvi,rhoi,hmixi,tri,tti,(mass(i,j),j=1,nspec)
      ! For switching coordinates: this happens in timemanager.f90 after the first fields are read
      if (xlonin.eq.-9999.9) exit
      call set_xlon(i,real((xlonin-xlon0)/dx,kind=dp))
      call set_ylat(i,real((ylatin-ylat0)/dy,kind=dp))
      numparticlecount=max(numparticlecount,part(i)%npoint)
    end do
  end do

  numpart=i-1

  close(unitpartin)

  julin=juldate(ibdatein,ibtimein)+real(itimein,kind=dp)/86400._dp
  if (abs(julin-bdate).gt.1.e-5) then
    write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
    write(*,*) ' #### ENDING TIME OF PREVIOUS MODEL RUN DOES   #### '
    write(*,*) ' #### NOT AGREE WITH STARTING TIME OF THIS RUN.#### '
    call caldate(julin,id1,it1)
    call caldate(bdate,id2,it2)
    write(*,*) 'julin: ',julin,id1,it1
    write(*,*) 'bdate: ',bdate,id2,it2
    stop
  end if
  do i=1,numpart
    julpartin=juldate(ibdatein,ibtimein)+ &
         real(part(i)%tstart,kind=dp)/86400._dp
    part(i)%nclass=min(int(ran1(idummy,0)*real(nclassunc))+1, &
         nclassunc)
    part(i)%idt=mintime
    part(i)%tstart=nint((julpartin-bdate)*86400.)
  end do

  return

998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
  write(*,*) ' #### '//path(2)(1:length(2))//'partposit'//' #### '
  write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS    #### '
  write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
  write(*,*) ' #### THE PROGRAM AGAIN.                       #### '
  stop

end subroutine readpartpositions

subroutine init_particle(itime,ipart,ithread)
  !                        i    i   o  o  o
  !        o       o       o    i  i  i   o
  !*****************************************************************************
  !                                                                            *
  !  Calculation of trajectories utilizing a zero-acceleration scheme. The time*
  !  step is determined by the Courant-Friedrichs-Lewy (CFL) criterion. This   *
  !  means that the time step must be so small that the displacement within    *
  !  this time step is smaller than 1 grid distance. Additionally, a temporal  *
  !  CFL criterion is introduced: the time step must be smaller than the time  *
  !  interval of the wind fields used for interpolation.                       *
  !  For random walk simulations, these are the only time step criteria.       *
  !  For the other options, the time step is also limited by the Lagrangian    *
  !  time scale.                                                               *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     16 December 1997                                                       *
  !                                                                            *
  !  Literature:                                                               *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! h [m]              Mixing height                                           *
  ! lwindinterv [s]    time interval between two wind fields                   *
  ! itime [s]          current temporal position                               *
  ! ldt [s]            Suggested time step for next integration                *
  ! ladvance [s]       Total integration time period                           *
  ! rannumb(maxrand)   normally distributed random variables                   *
  ! usig,vsig,wsig     uncertainties of wind velocities due to interpolation   *
  ! xt,yt,zt           Next time step's spatial position of trajectory         *
  !                                                                            *
  !                                                                            *
  ! Constants:                                                                 *
  ! cfl                factor, by which the time step has to be smaller than   *
  !                    the spatial CFL-criterion                               *
  ! cflt               factor, by which the time step has to be smaller than   *
  !                    the temporal CFL-criterion                              *
  !                                                                            *
  !*****************************************************************************

  use turbulence_mod
  use random_mod, only: ran3
  use omp_lib
  use interpol_mod
  use cbl_mod

  implicit none

  integer,intent(in) ::  &
    ithread,             & ! OMP thread starting at 0
    itime,               &
    ipart
  integer :: nrand
  real :: dummy1,dummy2
  real :: xt,yt,zt,zteta

  ! ! Initialise scavenging for backward runs
  ! !****************************************
  ! if (DRYBKDEP.or.WETBKDEP) then ! if there is no scavenging in wetdepo it will be set to 0
  !   do k=1,nspec
  !     xscav_frac1(ipart,k)=-1.
  !   end do
  ! endif


  part(ipart)%icbt=1           ! initialize particle to no "reflection"

  nrand=int(ran3(iseed1(ithread),ithread)*real(maxrand-1))+1

  xt = real(part(ipart)%xlon)
  yt = real(part(ipart)%ylat)
  zt = real(part(ipart)%z)
#ifdef ETA
  zteta = real(part(ipart)%zeta)
#else
  zteta = 0.
#endif

  !******************************
  ! 2. Interpolate necessary data
  !******************************

  ! Where in the grid? Stereographic (ngrid<0) or nested (ngrid>0)
  !***************************************************************
  call find_ngrid(xt,yt)
  ! Compute maximum mixing height around particle position
  !*******************************************************
  call find_grid_indices(xt,yt)

  ! If part_ic.nc is used, convert particle mass to mixing mass in forward,
  ! or internal mass conversion for wet and dry backward depostion.
  !************************************************************************
  if ((iout.eq.4).or.(iout.eq.5)) call init_mass_conversion(ipart,part(ipart)%npoint)
  
  h=max(hmix(ix ,jy,1,memind(1)), &
       hmix(ixp,jy ,1,memind(1)), &
       hmix(ix ,jyp,1,memind(1)), &
       hmix(ixp,jyp,1,memind(1)), &
       hmix(ix ,jy ,1,memind(2)), &
       hmix(ixp,jy ,1,memind(2)), &
       hmix(ix ,jyp,1,memind(2)), &
       hmix(ixp,jyp,1,memind(2)))

  zeta=zt/h


  !*************************************************************
  ! If particle is in the PBL, interpolate once and then make a
  ! time loop until end of interval is reached
  !*************************************************************

  if (zeta.le.1.) then

    call interpol_pbl(itime,xt,yt,zt,zteta,ithread+1)

  ! 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(zt,dummy1,dummy2,ithread+1)

  ! Compute the turbulent disturbances

  ! Determine the sigmas and the timescales
  !****************************************

    if (turbswitch) then
      call hanna(zt)
    else
      call hanna1(zt)
    endif


  ! Determine the new diffusivity velocities
  !*****************************************

    if (nrand+2.gt.maxrand) nrand=1
    part(ipart)%turbvel%u=rannumb(nrand)*sigu
    part(ipart)%turbvel%v=rannumb(nrand+1)*sigv
    part(ipart)%turbvel%w=rannumb(nrand+2)
    if (.not.turbswitch) then     ! modified by mc
      part(ipart)%turbvel%w=part(ipart)%turbvel%w*sigw
    else if (cblflag.eq.1) then   ! modified by mc
      if(-h/ol.gt.5) then
  !if (ol.lt.0.) then
  !if (ol.gt.0.) then !by mc : only for test correct is lt.0
        call init_cbl_vel( &
          iseed1(ithread),zt,wst,h,sigw,part(ipart)%turbvel%w,ol,ithread)
      else
        part(ipart)%turbvel%w=part(ipart)%turbvel%w*sigw
      end if
    end if


  ! 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), &
          600.) &
        *ctl)
    else
      part(ipart)%idt = int( &
        min(tlw, &
          h/max(2.*abs(part(ipart)%turbvel%w),1.e-5), &
          600.) &
        *ctl)
    endif
    part(ipart)%idt=max(part(ipart)%idt,mintime)

    ! call interpol_average()
    ! usig=(usigprof(indzp)+usigprof(indz))/2.
    ! vsig=(vsigprof(indzp)+vsigprof(indz))/2.
    ! wsig=(wsigprof(indzp)+wsigprof(indz))/2.

    ! wsigeta=(wsigprofeta(indzpeta)+wsigprofeta(indzeta))/2.

  else



  !**********************************************************
  ! 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
  !*********************

    call interpol_wind(itime,xt,yt,zt,zteta)


  ! Compute everything for above the PBL

  ! Assume constant turbulent perturbations
  !****************************************

    part(ipart)%idt=abs(lsynctime)

    if (nrand+1.gt.maxrand) nrand=1
    part(ipart)%turbvel%u=rannumb(nrand)*0.3
    part(ipart)%turbvel%v=rannumb(nrand+1)*0.3
    nrand=nrand+2
    part(ipart)%turbvel%w=0.
    sigw=0.

  endif

  !****************************************************************
  ! Add mesoscale random disturbances
  ! This is done only once for the whole lsynctime interval to save
  ! computation time
  !****************************************************************


  ! It is assumed that the average interpolation error is 1/2 sigma
  ! of the surrounding points, autocorrelation time constant is
  ! 1/2 of time interval between wind fields
  !****************************************************************
  if (lmesoscale_turb) then
    call interpol_mesoscale(xt,yt,zt,zteta)
    if (nrand+2.gt.maxrand) nrand=1
    part(ipart)%mesovel%u=rannumb(nrand)*usig
    part(ipart)%mesovel%v=rannumb(nrand+1)*vsig
#ifdef ETA
    part(ipart)%mesovel%w=rannumb(nrand+2)*wsigeta
#else
    part(ipart)%mesovel%w=rannumb(nrand+2)*wsig
#endif
  endif
end subroutine init_particle

subroutine init_domainfill
  !
  !*****************************************************************************
  !                                                                            *
  ! Initializes particles equally distributed over the first release location  *
  ! specified in file RELEASES. This box is assumed to be the domain for doing *
  ! domain-filling trajectory calculations.                                    *
  ! All particles carry the same amount of mass which alltogether comprises the*
  ! mass of air within the box.                                                *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     15 October 2002                                                        *
  !                                                                            *
  !  Changes                                                                   *
  !     2022, L. Bakels: OpenMP parallelisation                                *
  !     2023, L. Bakels: smooth vertical particle distribution instead of      *
  !                      distributing particles on fixed vertical layers       *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! numparticlecount    consecutively counts the number of particles released  *
  ! nx_we(2)       grid indices for western and eastern boundary of domain-    *
  !                filling trajectory calculations                             *
  ! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
  !                filling trajectory calculations                             *
  !                                                                            *
  !*****************************************************************************

  use point_mod
  use particle_mod

  implicit none

  integer :: j,kz,lix,ljy,ncolumn,numparttot,stat,iterminate
  real :: ylat,ylatp,ylatm,hzone
  real :: cosfactm,cosfactp,deltacol,dz1,dz2,dz,pnew,pnew_temp,fractus
  real,parameter :: pih=pi/180.
  real :: colmasstotal,zposition
  real,allocatable,dimension(:) :: pp

  integer :: ixm,ixp,jym,jyp,indzm,indzh,indzp,i,jj,ii
  ! integer :: alive_tmp,allocated_tmp,spawned_tmp,terminated_tmp
  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2)
  integer :: idummy = -11

  real :: height_tmp

  logical :: deall

  real,allocatable,dimension(:)   :: gridarea !
  real,allocatable,dimension(:,:) :: colmass !

  ! Determine the release region (only full grid cells), over which particles
  ! shall be initialized
  ! Use 2 fields for west/east and south/north boundary
  !**************************************************************************
  call alloc_domainfill

  nx_we(1)=max(int(xpoint1(1)),0)
  nx_we(2)=min((ceiling(xpoint2(1))),nxmin1)
  ny_sn(1)=max(int(ypoint1(1)),0)
  ny_sn(2)=min((ceiling(ypoint2(1))),nymin1)

  ! For global simulations (both global wind data and global domain-filling),
  ! set a switch, such that no boundary conditions are used
  !**************************************************************************
  if (xglobal.and.sglobal.and.nglobal) then
    if ((nx_we(1).eq.0).and.(nx_we(2).eq.nxmin1).and. &
         (ny_sn(1).eq.0).and.(ny_sn(2).eq.nymin1)) then
      gdomainfill=.true.
    else
      gdomainfill=.false.
    endif
  endif
  write(*,*) 'Global domain: ', gdomainfill

  ! Exit here if resuming a run from particle dump
  !***********************************************
  if (gdomainfill.and.ipin.ne.0) return

  ! Allocate grid and column mass
  !*******************************
  allocate(gridarea(0:nymax-1),colmass(0:nxmax-1,0:nymax-1),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate gridarea or colmass'

  ! Do not release particles twice (i.e., not at both in the leftmost and rightmost
  ! grid cell) for a global domain
  !*****************************************************************************
  if (xglobal) nx_we(2)=min(nx_we(2),nx-2)


  ! Calculate area of grid cell with formula M=2*pi*R*h*dx/360,
  ! see Netz, Formeln der Mathematik, 5. Auflage (1983), p.90
  !************************************************************
  ! First for the south pole
  
  if (sglobal) then
    ylat=ylat0
    ylatp=ylat+0.5*dy
    ylatm=ylat
    cosfactm=0.
    cosfactp=cos(ylatp*pih)*r_earth
    hzone=sqrt(r_earth**2-cosfactm**2)- &
         sqrt(r_earth**2-cosfactp**2)
    gridarea(0)=2.*pi*r_earth*hzone*dx/360.
  endif

  ! Do the same for the north pole

  if (nglobal) then
    ylat=ylat0+real(nymin1)*dy
    ylatp=ylat
    ylatm=ylat-0.5*dy
    cosfactp=0.
    cosfactm=cos(ylatm*pih)*r_earth
    hzone=sqrt(r_earth**2-cosfactp**2)- &
         sqrt(r_earth**2-cosfactm**2)
    gridarea(nymin1)=2.*pi*r_earth*hzone*dx/360.
  endif



  ! Allocate memory for storing the particles
  !******************************************
  call alloc_particles(int(npart(1)*1.1)) ! A bit more to avoid single part alloc

  ! Initialise total particle number
  numparttot=0
  ! Initialise max column number
  numcolumn=0

  ! Initialise the sum over the total mass of the atmosphere
  colmasstotal=0.

!$OMP PARALLEL PRIVATE(ljy,ylat,ylatp,ylatm,hzone,cosfactp,cosfactm,pp,lix) &
!$OMP REDUCTION(+:colmasstotal)
  
  allocate( pp(nzmax),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate pp inside of OMP loop'

!$OMP DO
  do ljy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    ylat=ylat0+real(ljy)*dy
    ylatp=ylat+0.5*dy
    ylatm=ylat-0.5*dy
    if ((ylatm.lt.0).and.(ylatp.gt.0.)) then
      hzone=1./dyconst
    else
      cosfactp=cos(ylatp*pih)*r_earth
      cosfactm=cos(ylatm*pih)*r_earth
      if (cosfactp.lt.cosfactm) then
        hzone=sqrt(r_earth**2-cosfactp**2)- &
             sqrt(r_earth**2-cosfactm**2)
      else
        hzone=sqrt(r_earth**2-cosfactm**2)- &
             sqrt(r_earth**2-cosfactp**2)
      endif
    endif
    gridarea(ljy)=2.*pi*r_earth*hzone*dx/360.
  end do
!$OMP END DO
!$OMP BARRIER

  ! Calculate total mass of each grid column and of the whole atmosphere
  !*********************************************************************
!$OMP DO 
  do ljy=ny_sn(1),ny_sn(2)          ! loop about latitudes
    do lix=nx_we(1),nx_we(2)      ! loop about longitudes
      pp(1)=prs(lix,ljy,1,1) !rho(lix,ljy,1,1)*r_air*tt(lix,ljy,1,1)
      pp(nz)=prs(lix,ljy,nz,1) !rho(lix,ljy,nz,1)*r_air*tt(lix,ljy,nz,1)
      colmass(lix,ljy)=(pp(1)-pp(nz))/ga*gridarea(ljy)
      colmasstotal=colmasstotal+colmass(lix,ljy)
    end do
  end do
!$OMP END DO
  deallocate(pp)
!$OMP END PARALLEL

  write(*,*) 'Atm. mass: ',colmasstotal

  allocate( pp(nzmax),stat=stat)
  if (stat.ne.0) write(*,*)'ERROR: could not allocate pp'

  if (ipin.eq.0) numpart=0

  ! Determine the particle positions
  !*********************************
  iterminate=0
  do ljy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    ylat=ylat0+real(ljy)*dy
    do lix=nx_we(1),nx_we(2)      ! loop about longitudes
      ncolumn=nint(0.999*real(npart(1))*colmass(lix,ljy)/colmasstotal)
      if (ncolumn.eq.0) cycle
      if (ncolumn.gt.numcolumn) numcolumn=ncolumn

  ! Calculate pressure at the altitudes of model surfaces, using the air density
  ! information, which is stored as a 3-d field
  !*****************************************************************************

      pp(:)=prs(lix,ljy,:,1)!rho(lix,ljy,kz,1)*r_air*tt(lix,ljy,kz,1)


      deltacol=(pp(1)-pp(nz))/real(ncolumn)
      pnew=pp(1)+deltacol*0.5
      jj=0
      do j=1,ncolumn ! looping over the number of particles within the column

        ! For columns with many particles (i.e. around the equator), distribute
        ! the particles equally (1 on a random position within the deltacol range), 
        ! for columns with few particles (i.e. around the poles), 
        ! distribute the particles randomly
        !***********************************************************************

        if ((ncolumn.gt.20).and.(ncolumn-j.gt.20)) then
          pnew_temp=pnew-ran1(idummy,0)*deltacol
          pnew=pnew-deltacol
        else if ((ncolumn.gt.20).and.(ncolumn-j.le.20)) then
          ! When only few particles are left, distribute them randomly above pnew
          pnew_temp=pnew-ran1(idummy,0)*(pnew-pp(nz))
        else
          pnew_temp=pp(1)-ran1(idummy,0)*(pp(1)-pp(nz))
        endif

        do kz=1,nz-1
          if ((pp(kz).ge.pnew_temp).and.(pp(kz+1).lt.pnew_temp)) then
            dz1=log(pnew_temp)-log(pp(kz))
            dz=1./log(pp(kz+1)/pp(kz))

            ! Assign particle position
            !*************************
            ! Do the following steps only if particles are not read in 
            ! from previous model run
            !**********************************************************
            if (ipin.eq.0) then
              ! First spawn the particle into existence
              !****************************************
              jj=jj+1
              !THIS WILL CAUSE PROBLEMS WITH OMP! because of dynamical allocation
              call spawn_particle(0,numpart+jj)
              call set_xlon(numpart+jj,real(real(lix)-0.5+ran1(idummy,0),kind=dp))
              if (lix.eq.0) call set_xlon(numpart+jj,real(ran1(idummy,0),kind=dp))
              if (lix.eq.nxmin1) &
                call set_xlon(numpart+jj,real(real(nxmin1)-ran1(idummy,0),kind=dp))
              call set_ylat(numpart+jj,real(real(ljy)-0.5+ran1(idummy,0),kind=dp))
              ! Logarithmic distribution of particles along pressure levels:
              ! hx=h1+(h2-h1)/log(p2/p1)*log(px/p1)
              height_tmp=height(kz)+(height(kz+1)-height(kz))*dz*dz1
              call set_z(numpart+jj,height_tmp)
              if (real(part(numpart+jj)%z).gt.height(nz)-0.5) &
                call set_z(numpart+jj, height(nz)-0.5)
#ifdef ETA
              call update_z_to_zeta(0, numpart+jj)
#endif
              ! Interpolate PV to the particle position
              !****************************************
              ixm=int(part(numpart+jj)%xlon)
              jym=int(part(numpart+jj)%ylat)
              ixp=ixm+1
              jyp=jym+1
              ddx=real(part(numpart+jj)%xlon)-real(ixm)
              ddy=real(part(numpart+jj)%ylat)-real(jym)
              rddx=1.-ddx
              rddy=1.-ddy
              p1=rddx*rddy
              p2=ddx*rddy
              p3=rddx*ddy
              p4=ddx*ddy

              !***********************************************
              indzm=nz-1
              indzp=nz
              do i=2,nz
                if (real(height(i),kind=dp).gt.part(numpart+jj)%z) then
                  indzm=i-1
                  indzp=i
                  exit
                endif
              end do
              dz1=real(part(numpart+jj)%z)-height(indzm)
              dz2=height(indzp)-real(part(numpart+jj)%z)
              dz=1./(dz1+dz2)
              do ii=1,2
                indzh=indzm+ii-1
                y1(ii)=p1*pv(ixm,jym,indzh,1) &
                     + p2*pv(ixp,jym,indzh,1) &
                     + p3*pv(ixm,jyp,indzh,1) &
                     + p4*pv(ixp,jyp,indzh,1)
              end do
              pvpart=(dz2*y1(1)+dz1*y1(2))*dz
              if (ylat.lt.0.) pvpart=-1.*pvpart


              ! For domain-filling option 2 (stratospheric O3), 
              ! do the rest only in the stratosphere
              !************************************************

              if (((part(numpart+jj)%z.gt.3000.).and. &
                   (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then

                ! Assign certain properties to the particle
                !******************************************
                part(numpart+jj)%nclass=min( &
                  int(ran1(idummy,0)*real(nclassunc))+1,nclassunc)
                numparticlecount=numparticlecount+1
                part(numpart+jj)%npoint=numparticlecount
                part(numpart+jj)%idt=mintime
                mass(numpart+jj,1)=colmass(lix,ljy)/real(ncolumn)
                if (mdomainfill.eq.2) mass(numpart+jj,1)= &
                     mass(numpart+jj,1)*pvpart*48./29.*ozonescale/10.**9
                mass_init(numpart+jj,1)=mass(numpart+jj,1)
              else
                call terminate_particle(numpart+jj, 0)
                jj=jj-1
                iterminate=iterminate+1
              endif
            endif
          endif
        end do
      end do
      numparttot=numparttot+ncolumn
      if (ipin.eq.0) numpart=numpart+jj
    end do
  end do


  ! alive_tmp=count%alive
  ! spawned_tmp=count%spawned
  ! terminated_tmp=count%terminated

! !$OMP PARALLEL PRIVATE(j) REDUCTION(+:alive_tmp,spawned_tmp,allocated_tmp,terminated_tmp)
 
  ! Make sure that all particles are within domain
  !***********************************************
! !$OMP DO
  do j=1,numpart
    if ((part(j)%xlon.lt.0.).or.(part(j)%xlon.ge.real(nxmin1,kind=dp)).or. &
         (part(j)%ylat.lt.0.).or.(part(j)%ylat.ge.real(nymin1,kind=dp))) then
      call terminate_particle(j,0) ! Cannot be within an OMP region
      iterminate=iterminate+1
      ! alive_tmp=alive_tmp-1
      ! terminated_tmp=terminated_tmp+1
    endif
  end do

  if (iterminate.gt.0) call rewrite_ialive()
! !$OMP END DO
! !$OMP END PARALLEL

  ! count%alive=alive_tmp
  ! count%spawned=spawned_tmp
  ! count%terminated=terminated_tmp
  ! Check whether numpart is really smaller than maxpart
  !*****************************************************

  ! ! ESO :TODO: this warning need to be moved further up, else out-of-bounds error earlier
  !   if (numpart.gt.maxpart) then
  !     write(*,*) 'numpart too large: change source in init_atm_mass.f'
  !     write(*,*) 'numpart: ',numpart,' maxpart: ',maxpart
  !   endif


  xmassperparticle=colmasstotal/real(numparttot)


  ! For boundary conditions, we need fewer particle release heights per column,
  ! because otherwise it takes too long until enough mass has accumulated to
  ! release a particle at the boundary (would take dx/u seconds), leading to
  ! relatively large position errors of the order of one grid distance.
  ! It's better to release fewer particles per column, but to do so more often.
  ! Thus, use on the order of nz starting heights per column.
  ! We thus repeat the above to determine fewer starting heights, that are
  ! used furtheron in subroutine boundcond_domainfill.f.
  !****************************************************************************

  fractus=real(numcolumn)/real(nz)
  write(*,*) 'Total number of particles at model start: ',numpart
  write(*,*) 'Maximum number of particles per column: ',numcolumn
  write(*,*) 'If ',fractus,' <1, better use more particles'
  fractus=sqrt(max(fractus,1.))*0.5

  do ljy=ny_sn(1),ny_sn(2)      ! loop about latitudes
    do lix=nx_we(1),nx_we(2)      ! loop about longitudes
      ncolumn=nint(0.999/fractus*real(npart(1))*colmass(lix,ljy) &
           /colmasstotal)
      if (ncolumn.gt.maxcolumn) stop 'maxcolumn too small'
      if (ncolumn.eq.0) cycle


  ! Memorize how many particles per column shall be used for all boundaries
  ! This is further used in subroutine boundcond_domainfill.f
  ! Use 2 fields for west/east and south/north boundary
  !************************************************************************

      if (lix.eq.nx_we(1)) numcolumn_we(1,ljy)=ncolumn
      if (lix.eq.nx_we(2)) numcolumn_we(2,ljy)=ncolumn
      if (ljy.eq.ny_sn(1)) numcolumn_sn(1,lix)=ncolumn
      if (ljy.eq.ny_sn(2)) numcolumn_sn(2,lix)=ncolumn

  ! Calculate pressure at the altitudes of model surfaces, using the air density
  ! information, which is stored as a 3-d field
  !*****************************************************************************

      do kz=1,nz
        pp(kz)=prs(lix,ljy,kz,1) !rho(lix,ljy,kz,1)*r_air*tt(lix,ljy,kz,1)
      end do

  ! Determine the reference starting altitudes
  !*******************************************

      deltacol=(pp(1)-pp(nz))/real(ncolumn)
      pnew=pp(1)+deltacol*0.5
      do j=1,ncolumn
        pnew=pnew-deltacol
        do kz=1,nz-1
          if ((pp(kz).ge.pnew).and.(pp(kz+1).lt.pnew)) then
            dz1=pp(kz)-pnew
            dz2=pnew-pp(kz+1)
            dz=1./(dz1+dz2)
            zposition=(height(kz)*dz2+height(kz+1)*dz1)*dz
            if (zposition.gt.height(nz)-0.5) zposition=height(nz)-0.5

  ! Memorize vertical positions where particles are introduced
  ! This is further used in subroutine boundcond_domainfill.f
  !***********************************************************

            if (lix.eq.nx_we(1)) zcolumn_we(1,ljy,j)=zposition
            if (lix.eq.nx_we(2)) zcolumn_we(2,ljy,j)=zposition
            if (ljy.eq.ny_sn(1)) zcolumn_sn(1,lix,j)=zposition
            if (ljy.eq.ny_sn(2)) zcolumn_sn(2,lix,j)=zposition

  ! Initialize mass that has accumulated at boundary to zero
  !*********************************************************

            acc_mass_we(1,ljy,j)=0.
            acc_mass_we(2,ljy,j)=0.
            acc_mass_sn(1,ljy,j)=0.
            acc_mass_sn(2,ljy,j)=0.
          endif
        end do
      end do
    end do
  end do

  ! If there were more particles allocated than used,
  ! Deallocate unused memory and update numpart
  !**************************************************
  deall=.false.
  do i=numpart, 1, -1
    if (.not. part(i)%alive) then
      deall=.true.
      numpart = numpart - 1
    else
      exit
    endif
  end do

  if (deall) call dealloc_particle(numpart) !Deallocates everything above numpart (F2008)


  ! If particles shall be read in to continue an existing run,
  ! then the accumulated masses at the domain boundaries must be read in, too.
  ! This overrides any previous calculations.
  !***************************************************************************

  if ((ipin.eq.1).and.(.not.gdomainfill)) then
    open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
         form='unformatted')
    read(unitboundcond) numcolumn_we,numcolumn_sn, &
         zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
    close(unitboundcond)
  endif

  deallocate(gridarea,colmass)
end subroutine init_domainfill

subroutine boundcond_domainfill(itime,loutend)
  !                                  i      i
  !*****************************************************************************
  !                                                                            *
  ! Particles are created by this subroutine continuously throughout the       *
  ! simulation at the boundaries of the domain-filling box.                    *
  ! All particles carry the same amount of mass which alltogether comprises the*
  ! mass of air within the box, which remains (more or less) constant.         *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     16 October 2002                                                        *
  !                                                                            *
  !  Changes                                                                   *
  !     2022, L. Bakels: OpenMP parallelisation                                *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  ! nx_we(2)       grid indices for western and eastern boundary of domain-    *
  !                filling trajectory calculations                             *
  ! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
  !                filling trajectory calculations                             *
  !                                                                            *
  !*****************************************************************************

  use point_mod
#ifdef _OPENMP
  use omp_lib
#endif
  implicit none

  real :: dz,dz1,dz2,dt1,dt2,dtt,ylat,cosfact,accmasst
  integer :: itime,ii,indz,indzp,i,loutend,numparticlecount_tmp
  integer :: j,k,ix,jy,m,indzh,indexh,ipart,mmass,ithread
  integer :: numactiveparticles,iterminate

  real :: windl(2),rhol(2)
  real :: windhl(2),rhohl(2)
  real :: windx,rhox
  real :: deltaz,boundarea,fluxofmass

  integer :: ixm,ixp,jym,jyp,indzm,mm,iterm_index
  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)

  integer :: idummy = -11


  ! If domain-filling is global, no boundary conditions are needed
  !***************************************************************

  if (gdomainfill) return

  ! Determine auxiliary variables for time interpolation
  !*****************************************************

  dt1=real(itime-memtime(1))
  dt2=real(memtime(2)-itime)
  dtt=1./(dt1+dt2)

  numactiveparticles=0
  numparticlecount_tmp=numparticlecount
  accmasst=0.
  ! Terminate trajectories that have left the domain, if domain-filling
  ! trajectory calculation domain is not global
  !********************************************************************
  if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()

  iterminate=0
  do i=1,count%allocated
    if (.not. part(i)%alive) cycle

    if ((part(i)%ylat.gt.real(ny_sn(2))).or.(part(i)%ylat.lt.real(ny_sn(1)))) then
      call terminate_particle(i,itime)
      iterminate=iterminate+1
    endif
    if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
         ((part(i)%xlon.lt.real(nx_we(1))).or. &
         (part(i)%xlon.gt.real(nx_we(2))))) then
      call terminate_particle(i,itime)
      iterminate=iterminate+1
    endif
    if (part(i)%alive) numactiveparticles = numactiveparticles+1
  end do
  if (iterminate.gt.0) call rewrite_ialive()
  !***************************************
  ! Western and eastern boundary condition
  !***************************************

  ! Loop from south to north
  !*************************
! OMP doesn't work here yet because particles are being spawned inside of the loop
! !$OMP PARALLEL PRIVATE(i,jy,k,j,ii,deltaz,boundarea,indz,indzp,indexh,windl,rhol, &
! !$OMP windhl,rhohl,windx,rhox,fluxofmass,mmass,ixm,jym,ixp,jyp,ddx,ddy,rddx, &
! !$OMP rddy,p1,p2,p3,p4,indzm,mm,indzh,pvpart,ylat,ix,cosfact,ipart) &
! !$OMP REDUCTION(+:numactiveparticles,numparticlecount_tmp,accmasst)

! #ifdef _OPENMP
!   ithread = OMP_GET_THREAD_NUM()
! #else
!   ithread = 0
! #endif
  ithread=0
  iterm_index=1
! !$OMP DO
  do jy=ny_sn(1),ny_sn(2)

  ! Loop over western (index 1) and eastern (index 2) boundary
  !***********************************************************

    do k=1,2

  ! Loop over all release locations in a column
  !********************************************

      do j=1,numcolumn_we(k,jy)

  ! Determine, for each release location, the area of the corresponding boundary
  !*****************************************************************************

        if (j.eq.1) then
          deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))*0.5
        else if (j.eq.numcolumn_we(k,jy)) then
  ! In order to avoid taking a very high column for very many particles,
  ! use the deltaz from one particle below instead
          deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))*0.5
        else
          deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))*0.5
        endif
        if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then
          boundarea=deltaz*111198.5/2.*dy
        else
          boundarea=deltaz*111198.5*dy
        endif


  ! Interpolate the wind velocity and density to the release location
  !******************************************************************

  ! Determine the model level below the release position
  !*****************************************************
        indz=nz-1
        indzp=nz
        do i=2,nz
          if (height(i).gt.zcolumn_we(k,jy,j)) then
            indz=i-1
            indzp=i
            exit
          endif
        end do

  ! Vertical distance to the level below and above current position
  !****************************************************************

        dz1=zcolumn_we(k,jy,j)-height(indz)
        dz2=height(indzp)-zcolumn_we(k,jy,j)
        dz=1./(dz1+dz2)

  ! Vertical and temporal interpolation
  !************************************

        do m=1,2
          indexh=memind(m)
          do ii=1,2
            indzh=indz+ii-1
            windl(ii)=uu(nx_we(k),jy,indzh,indexh)
            rhol(ii)=rho(nx_we(k),jy,indzh,indexh)
          end do

          windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
          rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
        end do

        windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
        rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt

  ! Calculate mass flux
  !********************

        fluxofmass=windx*rhox*boundarea*real(lsynctime)


  ! If the mass flux is directed into the domain, add it to previous mass fluxes;
  ! if it is out of the domain, set accumulated mass flux to zero
  !******************************************************************************

        if (k.eq.1) then
          if (fluxofmass.ge.0.) then
            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass
          else
            acc_mass_we(k,jy,j)=0.
          endif
        else
          if (fluxofmass.le.0.) then
            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass)
          else
            acc_mass_we(k,jy,j)=0.
          endif
        endif
        accmasst=accmasst+acc_mass_we(k,jy,j)

  ! If the accumulated mass exceeds half the mass that each particle shall carry,
  ! one (or more) particle(s) is (are) released and the accumulated mass is
  ! reduced by the mass of this (these) particle(s)
  !******************************************************************************

        if (acc_mass_we(k,jy,j).ge.xmassperparticle*0.5) then
          mmass=int((acc_mass_we(k,jy,j)+xmassperparticle*0.5)/ &
               xmassperparticle)
          acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
               real(mmass)*xmassperparticle
        else
          mmass=0
        endif

        do m=1,mmass
          !THIS WILL CAUSE PROBLEMS WITH OMP! because of dynamical allocation
          call get_newpart_index(ipart,iterm_index)
          call spawn_particle(itime, ipart)

  ! Assign particle positions
  !**************************

          call set_xlon(ipart,real(nx_we(k),kind=dp))
          if (jy.eq.ny_sn(1)) then
            call set_ylat(ipart,real(real(jy)+0.5*ran1(idummy,ithread),kind=dp))
          else if (jy.eq.ny_sn(2)) then
            call set_ylat(ipart,real(real(jy)-0.5*ran1(idummy,ithread),kind=dp))
          else
            call set_ylat(ipart,real(real(jy)+(ran1(idummy,ithread)-.5),kind=dp))
          endif
          if (j.eq.1) then
            call set_z(ipart,zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- &
                 zcolumn_we(k,jy,1))/4.)
          else if (j.eq.numcolumn_we(k,jy)) then
            call set_z(ipart,(2.*zcolumn_we(k,jy,j)+ &
                 zcolumn_we(k,jy,j-1)+height(nz))/4.)
          else
            call set_z(ipart,zcolumn_we(k,jy,j-1)+ran1(idummy,ithread)* &
                 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1)))
          endif
#ifdef ETA
          call update_z_to_zeta(itime, ipart)
#endif
  ! Interpolate PV to the particle position
  !****************************************
          ixm=int(part(ipart)%xlon)
          jym=int(part(ipart)%ylat)
          ixp=ixm+1
          jyp=jym+1
          ddx=real(part(ipart)%xlon)-real(ixm)
          ddy=real(part(ipart)%ylat)-real(jym)
          rddx=1.-ddx
          rddy=1.-ddy
          p1=rddx*rddy
          p2=ddx*rddy
          p3=rddx*ddy
          p4=ddx*ddy
          indzm=nz-1
          indzp=nz
          do i=2,nz
            if (real(height(i),kind=dp).gt.part(ipart)%z) then
              indzm=i-1
              indzp=i
              exit
            endif
          end do
          dz1=real(part(ipart)%z)-height(indzm)
          dz2=height(indzp)-real(part(ipart)%z)
          dz=1./(dz1+dz2)
          do mm=1,2
            indexh=memind(mm)
            do ii=1,2
              indzh=indzm+ii-1
              y1(ii)=p1*pv(ixm,jym,indzh,indexh) &
                   +p2*pv(ixp,jym,indzh,indexh) &
                   +p3*pv(ixm,jyp,indzh,indexh) &
                   +p4*pv(ixp,jyp,indzh,indexh)
            end do
            yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
          end do
          pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
          ylat=ylat0+real(part(ipart)%ylat)*dy
          if (ylat.lt.0.) pvpart=-1.*pvpart


  ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
  !*****************************************************************************

          if (((part(ipart)%z.gt.3000.).and. &
               (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
            part(ipart)%nclass=min(int(ran1(idummy,ithread)* &
                 real(nclassunc))+1,nclassunc)
            numactiveparticles=numactiveparticles+1
            numparticlecount_tmp=numparticlecount_tmp+1
            part(ipart)%npoint=numparticlecount_tmp
            part(ipart)%idt=mintime
            part(ipart)%tstart=itime
            mass(ipart,1)=xmassperparticle
            if (mdomainfill.eq.2) mass(ipart,1)= &
                 mass(ipart,1)*pvpart*48./29.*ozonescale/10.**9
            mass_init(ipart,1)=mass(ipart,1)
          else
            stop 'boundcond_domainfill error: look into original to understand what should happen here'
          endif
        end do ! particles
      end do ! release locations in column
    end do ! western and eastern boundary
  end do ! south to north
! !$OMP END DO

  !*****************************************
  ! Southern and northern boundary condition
  !*****************************************

  ! Loop from west to east
  !***********************
! !$OMP DO
  do ix=nx_we(1),nx_we(2)

  ! Loop over southern (index 1) and northern (index 2) boundary
  !*************************************************************

    do k=1,2
      ylat=ylat0+real(ny_sn(k))*dy
      cosfact=cos(ylat*pi180)

  ! Loop over all release locations in a column
  !********************************************

      do j=1,numcolumn_sn(k,ix)

  ! Determine, for each release location, the area of the corresponding boundary
  !*****************************************************************************

        if (j.eq.1) then
          deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))*0.5
        else if (j.eq.numcolumn_sn(k,ix)) then
  !        deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+
  !    +        zcolumn_sn(k,ix,j))/2.
  ! In order to avoid taking a very high column for very many particles,
  ! use the deltaz from one particle below instead
          deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))*0.5
        else
          deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))*0.5
        endif
        if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then
          boundarea=deltaz*111198.5*0.5*cosfact*dx
        else
          boundarea=deltaz*111198.5*cosfact*dx
        endif


  ! Interpolate the wind velocity and density to the release location
  !******************************************************************

  ! Determine the model level below the release position
  !*****************************************************
        indz=nz-1
        indzp=nz
        do i=2,nz
          if (height(i).gt.zcolumn_sn(k,ix,j)) then
            indz=i-1
            indzp=i
            exit
          endif
        end do

  ! Vertical distance to the level below and above current position
  !****************************************************************

        dz1=zcolumn_sn(k,ix,j)-height(indz)
        dz2=height(indzp)-zcolumn_sn(k,ix,j)
        dz=1./(dz1+dz2)

  ! Vertical and temporal interpolation
  !************************************

        do m=1,2
          indexh=memind(m)
          do ii=1,2
            indzh=indz+ii-1
            windl(ii)=vv(ix,ny_sn(k),indzh,indexh)
            rhol(ii)=rho(ix,ny_sn(k),indzh,indexh)
          end do

          windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
          rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
        end do

        windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
        rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt

  ! Calculate mass flux
  !********************

        fluxofmass=windx*rhox*boundarea*real(lsynctime)

  ! If the mass flux is directed into the domain, add it to previous mass fluxes;
  ! if it is out of the domain, set accumulated mass flux to zero
  !******************************************************************************

        if (k.eq.1) then
          if (fluxofmass.ge.0.) then
            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass
          else
            acc_mass_sn(k,ix,j)=0.
          endif
        else
          if (fluxofmass.le.0.) then
            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass)
          else
            acc_mass_sn(k,ix,j)=0.
          endif
        endif
        accmasst=accmasst+acc_mass_sn(k,ix,j)

  ! If the accumulated mass exceeds half the mass that each particle shall carry,
  ! one (or more) particle(s) is (are) released and the accumulated mass is
  ! reduced by the mass of this (these) particle(s)
  !******************************************************************************

        if (acc_mass_sn(k,ix,j).ge.xmassperparticle*0.5) then
          mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle*0.5)/ &
               xmassperparticle)
          acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
               real(mmass)*xmassperparticle
        else
          mmass=0
        endif

        do m=1,mmass
          call get_newpart_index(ipart,iterm_index)
          call spawn_particle(itime, ipart)
  
  ! Assign particle positions
  !**************************
          call set_ylat(ipart,real(ny_sn(k),kind=dp))
          if (ix.eq.nx_we(1)) then
            call set_xlon(ipart,real(real(ix)+0.5*ran1(idummy,ithread),kind=dp))
          else if (ix.eq.nx_we(2)) then
            call set_xlon(ipart,real(real(ix)-0.5*ran1(idummy,ithread),kind=dp))
          else
            call set_xlon(ipart,real(real(ix)+(ran1(idummy,ithread)-.5),kind=dp))
          endif
          if (j.eq.1) then
            call set_z(ipart,zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- &
                 zcolumn_sn(k,ix,1))/4.)
          else if (j.eq.numcolumn_sn(k,ix)) then
            call set_z(ipart,(2.*zcolumn_sn(k,ix,j)+ &
                 zcolumn_sn(k,ix,j-1)+height(nz))/4.)
          else
            call set_z(ipart,zcolumn_sn(k,ix,j-1)+ran1(idummy,ithread)* &
                 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1)))
          endif
#ifdef ETA
          call update_z_to_zeta(itime, ipart)
#endif
  ! Interpolate PV to the particle position
  !****************************************
          ixm=int(part(ipart)%xlon)
          jym=int(part(ipart)%ylat)
          ixp=ixm+1
          jyp=jym+1
          ddx=real(part(ipart)%xlon)-real(ixm)
          ddy=real(part(ipart)%ylat)-real(jym)
          rddx=1.-ddx
          rddy=1.-ddy
          p1=rddx*rddy
          p2=ddx*rddy
          p3=rddx*ddy
          p4=ddx*ddy
          indzm=nz-1
          indzp=nz
          do i=2,nz
            if (real(height(i),kind=dp).gt.part(ipart)%z) then
              indzm=i-1
              indzp=i
              exit
            endif
          end do
          dz1=real(part(ipart)%z)-height(indzm)
          dz2=height(indzp)-real(part(ipart)%z)
          dz=1./(dz1+dz2)
          do mm=1,2
            indexh=memind(mm)
            do ii=1,2
              indzh=indzm+ii-1
              y1(ii)=p1*pv(ixm,jym,indzh,indexh) &
                   +p2*pv(ixp,jym,indzh,indexh) &
                   +p3*pv(ixm,jyp,indzh,indexh) &
                   +p4*pv(ixp,jyp,indzh,indexh)
            end do
            yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
          end do
          pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
          if (ylat.lt.0.) pvpart=-1.*pvpart


  ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
  !*****************************************************************************

          if (((part(ipart)%z.gt.3000.).and. &
               (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
            part(ipart)%nclass=min(int(ran1(idummy,ithread)* &
                 real(nclassunc))+1,nclassunc)
            numactiveparticles=numactiveparticles+1
            numparticlecount_tmp=numparticlecount_tmp+1
            part(ipart)%npoint=numparticlecount_tmp
            part(ipart)%idt=mintime
            mass(ipart,1)=xmassperparticle
            if (mdomainfill.eq.2) mass(ipart,1)= &
                 mass(ipart,1)*pvpart*48./29.*ozonescale/10.**9
            mass_init(ipart,1)=mass(ipart,1)
          else
            stop 'boundcond_domainfill error: look into original to understand what should happen here'
          endif
        end do ! particles
      end do ! releases per column
    end do ! east west
  end do ! north south
! !$OMP END DO
! !$OMP END PARALLEL
  if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
  numparticlecount = numparticlecount_tmp
  ! If particles shall be dumped, then accumulated masses at the domain boundaries
  ! must be dumped, too, to be used for later runs
  !*****************************************************************************

  if ((ipout.gt.0).and.(itime.eq.loutend)) then
    open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
         form='unformatted')
    write(unitboundcond) numcolumn_we,numcolumn_sn, &
         zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
    close(unitboundcond)
  endif
end subroutine boundcond_domainfill

end module initialise_mod
