readpartpositions.f90 4.88 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
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                                                          *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  !                                                                            *
  !*****************************************************************************

  use par_mod
  use com_mod
21
  use random_mod
22 23 24 25 26

  implicit none

  integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,ix
  integer :: id1,id2,it1,it2
27
  real :: xlonin,ylatin,topo,hmixi,pvi,qvi,rhoi,tri,tti
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
  character :: specin*7
  real(kind=dp) :: julin,julpartin,juldate

  integer :: idummy = -8

  numparticlecount=0

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

  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)) goto 997

  do i=1,nspecin
    read(unitpartin)
    read(unitpartin)
    read(unitpartin) j,specin
    if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
  end do

  read(unitpartin) numpointin
59 60
  if (numpointin.ne.numpoint) goto 995
999 continue 
61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85
  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 ix=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)
86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
  

100 read(unitpartin,end=99) itimein
  i=0
200 i=i+1
  read(unitpartin) npoint(i),xlonin,ylatin,ztra1(i),itramem(i), &
       topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1(i,j),j=1,nspec)
  
  if (xlonin.eq.-9999.9) goto 100
  xtra1(i)=(xlonin-xlon0)/dx
  ytra1(i)=(ylatin-ylat0)/dy
  numparticlecount=max(numparticlecount,npoint(i))
  goto 200

99 numpart=i-1
101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128

  close(unitpartin)

  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

129 130 131 132 133 134
!995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
995   write(*,*) ' #### FLEXPART MODEL WARNING IN READPARTPOSITIONS#### '
  write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
  write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
!  stop
  goto 999 
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153

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