source: branches/jerome/src_flexwrf_v3.1/readpartpositions.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1
File size: 10.3 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23
24      subroutine readpartpositions
25!*******************************************************************************
26!                                                                              *
27!   Note:  This is the FLEXPART_WRF version of subroutine readpartpositions.   *
28!                                                                              *
29!   This routine opens the particle dump file and reads all the particle       *
30!   positions from a previous run to initialize the current run.               *
31!                                                                              *
32!                                                                              *
33!     Author: A. Stohl                                                         *
34!     24 March 2000                                                            *
35!                                                                              *
36!     Dec 2005, R. Easter                                                      *
37!             Changed names of "*lon0*" & "*lat0*" variables                   *
38!             Reads either binary or ascii output files from previous run.     *
39!             Particle positions may be in lat-lon or grid-meter units.        *
40!                                                                              *
41!*******************************************************************************
42!                                                                              *
43! Variables:                                                                   *
44!                                                                              *
45!*******************************************************************************
46
47!      include 'includepar'
48!     include 'includecom'
49  use par_mod
50  use com_mod
51
52  implicit none
53     
54      integer :: ibdatein,ibtimein,nspecin,itimein,numpointin,i,j,ix
55      integer :: iomode_xycoord_in, numpart_in
56      integer :: itmp,ntmp, numxgridin,numygridin
57      real :: xlonin,ylatin,ran1,topo,hmixi,pvi,qvi,rhoi,tri,tti
58      real :: xtmp
59      character :: specin*7
60      character :: ctmp*1
61  real(kind=dp) :: julin,julpartin,juldate
62  integer :: idummy = -8
63
64
65! Open and read header file of dumped particle data
66!*****************************************
67
68      if (iouttype .eq. 0) then
69
70      open(unitpartin,file=path(1)(1:length(1))//'header', &
71      form='unformatted',err=998)
72
73      read(unitpartin) ibdatein,ibtimein
74      read(unitpartin)
75      read(unitpartin)
76
77      read(unitpartin)
78      read(unitpartin)
79      read(unitpartin) nspecin
80      nspecin=nspecin/3
81      if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997
82
83      do i=1,nspecin
84        read(unitpartin)
85        read(unitpartin)
86        read(unitpartin) j,specin
87        if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
88      enddo
89 
90      read(unitpartin) numpointin
91!     if (numpointin.ne.numpoint) goto 995
92      do i=1,numpointin
93        read(unitpartin)
94        read(unitpartin)
95        read(unitpartin)
96        read(unitpartin)
97        do j=1,nspec
98          read(unitpartin)
99          read(unitpartin)
100          read(unitpartin)
101        enddo
102      enddo
103      read(unitpartin)
104      read(unitpartin)
105
106      do ix=0,numxgrid-1
107        read(unitpartin)
108      enddo
109
110      close(unitpartin)
111
112      else    ! (iouttype .eq. 1)
113
114      open(unitpartin,file=path(1)(1:length(1))//'header', &
115      form='formatted',err=998)
116
117! with formatted header file, need to read every variable
118! because some of the "writes" cover multiple lines
119! (and the number of lines is compiler dependent)
120      read(unitpartin,*) ibdatein,ibtimein
121      read(unitpartin,*) ctmp   ! version id
122      read(unitpartin,*) itmp,itmp,itmp,itmp   ! loutstep,loutaver,...
123      read(unitpartin,*) xtmp,xtmp,numxgridin,numygridin,xtmp,xtmp   ! outgrid x,y info
124
125      read(unitpartin,*) ntmp,(xtmp, i=1,ntmp)   ! numzgrid,(outheight(i),i=1,numzgrid)
126      read(unitpartin,*) itmp,itmp   ! jjjjmmdd,ihmmss
127      read(unitpartin,*) nspecin,itmp,itmp
128      nspecin=nspecin/3
129      if ((ldirect.eq.1).and.(nspec.ne.nspecin)) goto 997
130
131      do i=1,nspecin
132        read(unitpartin,*) itmp   ! 1
133        read(unitpartin,*) ctmp   ! "WD" name
134        read(unitpartin,*) itmp   ! 1
135        read(unitpartin,*) ctmp   ! "DD" name
136        read(unitpartin,*) j
137        read(unitpartin,*) specin
138        if ((ldirect.eq.1).and.(species(i)(1:7).ne.specin)) goto 996
139      enddo
140 
141      read(unitpartin,*) numpointin
142!     if (numpointin.ne.numpoint) goto 995
143      do i=1,numpointin
144        read(unitpartin,*) itmp,itmp,itmp   ! release start,end,kindz
145        read(unitpartin,*) xtmp,xtmp,xtmp,xtmp,xtmp,xtmp   ! release x,y,z info
146        read(unitpartin,*) itmp,itmp   ! npart(i), 1
147        read(unitpartin,*) ctmp   ! compoint(i)
148        do j=1,nspec
149          read(unitpartin,*) xtmp   ! xmass(i,j)
150          read(unitpartin,*) xtmp   ! xmass(i,j)
151          read(unitpartin,*) xtmp   ! xmass(i,j)
152        enddo
153      enddo
154      read(unitpartin,*) itmp,itmp,itmp   ! method,lsubgrid,lconvection
155      read(unitpartin,*) ntmp,(itmp, i=1,ntmp)   ! nageclass,(lage(i),i=1,nageclass)
156
157      do ix=0,numxgridin-1
158        read(unitpartin,*) (xtmp, j=0,numygridin-1)    ! oroout
159      enddo
160
161      close(unitpartin)
162
163      endif   ! (iouttype .eq. 0/1)
164
165
166! Open and read data file of dumped particle data
167!***************************************
168
169      if (iouttype .eq. 0) then
170        open(unitpartin,file=path(1)(1:length(1))//'partposit_end', &
171          form='unformatted',err=998)
172      else
173        open(unitpartin,file=path(1)(1:length(1))//'partposit_end', &
174          form='formatted',err=998)
175      endif
176
177100   continue
178      if (iouttype .eq. 0) then
179        read(unitpartin,end=99) itimein,numpart_in, &
180          iomode_xycoord_in
181      else
182        read(unitpartin,*,end=99) itimein,numpart_in, &
183          iomode_xycoord_in
184      endif
185
186! iomode_xycoord of previous & current runs must match
187      if (iomode_xycoord_in .ne. outgrid_option) then
188         write(*,'(/a/a/)') '*** readpartpositions fatal error', &
189           'outgrid_option from previous & current runs differ'
190         stop
191      end if
192
193      i=0
194200   i=i+1
195      if (iouttype .eq. 0) then
196        read(unitpartin) npoint(i),xlonin,ylatin,ztra1(i),itramem(i), &
197          topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1(i,j),j=1,nspec)
198      else
199        read(unitpartin,*) npoint(i),itramem(i),xlonin,ylatin,ztra1(i), &
200          topo,pvi,qvi,rhoi,hmixi,tri,tti,(xmass1(i,j),j=1,nspec)
201      endif
202         
203      if (xlonin.eq.-9999.9) goto 100
204
205      if (outgrid_option .eq. 1) then
206! convert from lat-lon to grid-index coordinates
207        call ll_to_xyindex_wrf( xlonin, ylatin, xtra1(i), ytra1(i) )
208      else
209! convert from grid-meter to grid-index coordinates
210        xtra1(i)=(xlonin-xmet0)/dx
211        ytra1(i)=(ylatin-ymet0)/dy
212      endif
213      goto 200
214
21599    numpart=i-1
216
217      close(unitpartin)
218
219
220! Set nclass, idt, itra1, itramem, itrasplit to be consistent
221! with current run
222!***************************************
223
224  julin=juldate(ibdatein,ibtimein)+real(itimein,kind=dp)/86400._dp
225
226      if (abs(julin-bdate).gt.1.e-5) goto 994
227      do i=1,numpart
228        julpartin=juldate(ibdatein,ibtimein)+ &
229         real(itramem(i),kind=dp)/86400._dp
230        nclass(i)=min(int(ran1(idummy)*real(nclassunc))+1, &
231        nclassunc)
232        idt(i)=mintime
233        itra1(i)=0
234        itramem(i)=nint((julpartin-bdate)*86400.)
235        itrasplit(i)=ldirect*itsplit
236      enddo
237      return
238
239
240994   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
241      write(*,*) ' #### ENDING TIME OF PREVIOUS MODEL RUN DOES   #### '
242      write(*,*) ' #### NOT AGREE WITH STARTING TIME OF THIS RUN.#### '
243      write(*,*) 'julin: ',julin
244      write(*,*) 'bdate: ',bdate
245      stop
246
247995   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
248      write(*,*) ' #### NUMBER OF RELEASE LOCATIONS DOES NOT     #### '
249      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
250      stop
251
252996   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
253      write(*,*) ' #### SPECIES NAMES TO BE READ IN DO NOT       #### '
254      write(*,*) ' #### AGREE WITH CURRENT SETTINGS!             #### '
255      stop
256
257997   write(*,*) ' #### FLEXPART MODEL ERROR IN READPARTPOSITIONS#### '
258      write(*,*) ' #### THE NUMBER OF SPECIES TO BE READ IN DOES #### '
259      write(*,*) ' #### NOT AGREE WITH CURRENT SETTINGS!         #### '
260      stop
261
262998   write(*,*) ' #### FLEXPART MODEL ERROR!   THE FILE         #### '
263      write(*,*) ' #### '//path(1)(1:length(1))//'grid'//' #### '
264      write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS    #### '
265      write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### '
266      write(*,*) ' #### THE PROGRAM AGAIN.                       #### '
267      stop
268
269end subroutine readpartpositions
270
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG