readpartpositions_mpi.f90 6.97 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
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                                                                  *
  !     12/2014 eso: MPI version                                               *
  !                  Root process reads positions and distributes the data     *
  !                  :TODO: The above..                                        *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod
  use random_mod, only: ran1
26
  use mpi_mod !, only: mp_seed
27 28 29

  implicit none

30
  integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,ix,ip
31
  integer :: id1,id2,it1,it2
32
  integer :: addone,numparticlecount_all,numpart_all,lbnd,ubnd
33 34 35 36 37 38
  real :: xlonin,ylatin,topo,hmixi,pvi,qvi,rhoi,tri,tti
  character :: specin*7
  real(kind=dp) :: julin,julpartin,juldate

  integer :: idummy = -8

39 40 41 42 43
  ! These variables are allocated at the root process for all particles in file.
  real,dimension(maxpart) :: xtra1_all,ytra1_all,ztra1_all
  real,dimension(maxpart,maxspec) :: xmass1_all
  integer,dimension(maxpart) :: npoint_all,itramem_all

44 45 46 47
  ! Different seed for each process
  idummy=idummy+mp_seed

  numparticlecount=0
48 49
  numparticlecount_all=0
  numpart_all=0
50 51

  ! Open header file of dumped particle data
52 53
  ! Each MPI process sequentially access file (just in case)
  !*********************************************************
54

55 56 57 58
  do ip=0, mp_partgroup_np-1
    call mpif_mpi_barrier
    open(unitpartin,file=path(2)(1:length(2))//'header', &
         form='unformatted',err=998)
59

60
    read(unitpartin) ibdatein,ibtimein
61 62 63 64 65
    read(unitpartin)
    read(unitpartin)

    read(unitpartin)
    read(unitpartin)
66 67 68 69 70
    read(unitpartin) nspecin
    nspecin=nspecin/3
    if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997

    do i=1,nspecin
71 72
      read(unitpartin)
      read(unitpartin)
73 74
      read(unitpartin) j,specin
      if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
75
    end do
76 77 78 79 80 81 82 83 84 85 86 87
    
    read(unitpartin) numpointin
    if (numpointin.ne.numpoint) then ! goto 995

! eso 2016: moved this warning here to avoid out-of-block goto
!995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
      write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
      write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
!  stop
      goto 999 
    end if
88

89 90 91 92 93 94 95 96 97 98 99 100
999 continue 
    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
101
    read(unitpartin)
102 103 104 105 106 107
    read(unitpartin)

    do ix=0,numxgrid-1
      read(unitpartin)
    end do
    close(unitpartin)
108 109 110


  ! Open data file of dumped particle data
111
  ! All processes read the whole file
112 113
  !***************************************

114
    open(unitpartin,file=path(2)(1:length(2))//'partposit_end', &
115 116 117
         form='unformatted',err=998)
    
100 read(unitpartin,end=99) itimein
118
    i=0
119 120 121 122
200 i=i+1
    read(unitpartin) npoint_all(i),xlonin,ylatin,ztra1_all(i),itramem_all(i), &
         topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1_all(i,j),j=1,nspec)
    
123
    if (xlonin.eq.-9999.9) goto 100
124 125 126
    xtra1_all(i)=(xlonin-xlon0)/dx
    ytra1_all(i)=(ylatin-ylat0)/dy
    numparticlecount_all=max(numparticlecount_all,npoint(i))
127 128
    goto 200

129 130 131
99  numpart_all=i-1
    
    close(unitpartin)
132

133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172
  end do


  ! Distribute particles among processes
  !**************************************
  lbnd=1
  ubnd=0

  do ip=0, mp_partgroup_np-1
    
! Extra particle distributed in case remainder in the division
    if (ip.lt.mod(numpart_all,mp_partgroup_np)) then
      addone=1
    else
      addone=0
    end if

    if (ip==0) then 
      ubnd=ubnd + numpart_all/mp_partgroup_np + addone
    else 
      ubnd=lbnd + numpart_all/mp_partgroup_np + addone - 1
    end if
    
    if (ip==mp_pid) then
      
      numparticlecount=numparticlecount_all/mp_partgroup_np+addone
      numpart=numpart_all/mp_partgroup_np+addone

      xtra1(1:numpart) = xtra1_all(lbnd:ubnd)
      ytra1(1:numpart) = ytra1_all(lbnd:ubnd)
      ztra1(1:numpart) = ztra1_all(lbnd:ubnd)
      xmass1(1:numpart,:) = xmass1_all(lbnd:ubnd,:)
      npoint(1:numpart) = npoint_all(lbnd:ubnd)
      itramem(1:numpart) = itramem_all(lbnd:ubnd)

    end if

    lbnd=ubnd+1

  end do
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216

  julin=juldate(ibdatein,ibtimein)+real(itimein,kind=dp)/86400._dp
  if (abs(julin-bdate).gt.1.e-5) goto 994
  do i=1,numpart
    julpartin=juldate(ibdatein,ibtimein)+ &
         real(itramem(i),kind=dp)/86400._dp
    nclass(i)=min(int(ran1(idummy)*real(nclassunc))+1, &
         nclassunc)
    idt(i)=mintime
    itra1(i)=0
    itramem(i)=nint((julpartin-bdate)*86400.)
    itrasplit(i)=ldirect*itsplit
  end do

  return


994   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

996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
  write(*,*) ' #### SPECIES NAMES TO BE READ IN DO NOT       #### '
  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
  stop

997   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
  write(*,*) ' #### THE NUMBER OF SPECIES TO BE READ IN DOES #### '
  write(*,*) ' #### NOT AGREE WITH CURRENT SETTINGS!         #### '
  stop

998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
  write(*,*) ' #### '//path(2)(1:length(2))//'grid'//' #### '
  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