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
26
27
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
59
60
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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
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
129
130
131
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
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
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
|
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
! *
! This file is part of FLEXPART. *
! *
! FLEXPART is free software: you can redistribute it and/or modify *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or *
! (at your option) any later version. *
! *
! FLEXPART is distributed in the hope that it will be useful, *
! but WITHOUT ANY WARRANTY; without even the implied warranty of *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
! GNU General Public License for more details. *
! *
! You should have received a copy of the GNU General Public License *
! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
!**********************************************************************
subroutine readavailable
!*****************************************************************************
! *
! This routine reads the dates and times for which windfields are *
! available. *
! *
! Authors: A. Stohl *
! *
! 6 February 1994 *
! 8 February 1999, Use of nested fields, A. Stohl *
! *
!*****************************************************************************
! *
! Variables: *
! bdate beginning date as Julian date *
! beg beginning date for windfields *
! end ending date for windfields *
! fname filename of wind field, help variable *
! ideltas [s] duration of modelling period *
! idiff time difference between 2 wind fields *
! idiffnorm normal time difference between 2 wind fields *
! idiffmax [s] maximum allowable time between 2 wind fields *
! jul julian date, help variable *
! numbwf actual number of wind fields *
! wfname(maxwf) file names of needed wind fields *
! wfspec(maxwf) file specifications of wind fields (e.g., if on disc) *
! wftime(maxwf) [s]times of wind fields relative to beginning time *
! wfname1,wfspec1,wftime1 = same as above, but only local (help variables) *
! *
! Constants: *
! maxwf maximum number of wind fields *
! unitavailab unit connected to file AVAILABLE *
! *
!*****************************************************************************
use par_mod
use com_mod
implicit none
integer :: i,idiff,ldat,ltim,wftime1(maxwf),numbwfn(maxnests),k
integer :: wftime1n(maxnests,maxwf),wftimen(maxnests,maxwf)
real(kind=dp) :: juldate,jul,beg,end
character(len=255) :: fname,spec,wfname1(maxwf),wfspec1(maxwf)
character(len=255) :: wfname1n(maxnests,maxwf)
character(len=40) :: wfspec1n(maxnests,maxwf)
! Windfields are only used, if they are within the modelling period.
! However, 1 additional day at the beginning and at the end is used for
! interpolation. -> Compute beginning and ending date for the windfields.
!************************************************************************
if (ideltas.gt.0) then ! forward trajectories
beg=bdate-1._dp
end=bdate+real(ideltas,kind=dp)/86400._dp+real(idiffmax,kind=dp)/ &
86400._dp
else ! backward trajectories
beg=bdate+real(ideltas,kind=dp)/86400._dp-real(idiffmax,kind=dp)/ &
86400._dp
end=bdate+1._dp
endif
! Open the wind field availability file and read available wind fields
! within the modelling period.
!*********************************************************************
open(unitavailab,file=path(4)(1:length(4)),status='old', &
err=999)
do i=1,3
read(unitavailab,*)
end do
numbwf=0
100 read(unitavailab,'(i8,1x,i6,2(6x,a255))',end=99) &
ldat,ltim,fname,spec
jul=juldate(ldat,ltim)
if ((jul.ge.beg).and.(jul.le.end)) then
numbwf=numbwf+1
if (numbwf.gt.maxwf) then ! check exceedance of dimension
write(*,*) 'Number of wind fields needed is too great.'
write(*,*) 'Reduce modelling period (file "COMMAND") or'
write(*,*) 'reduce number of wind fields (file "AVAILABLE").'
stop
endif
wfname1(numbwf)=fname(1:index(fname,' '))
wfspec1(numbwf)=spec
wftime1(numbwf)=nint((jul-bdate)*86400._dp)
endif
goto 100 ! next wind field
99 continue
close(unitavailab)
! Open the wind field availability file and read available wind fields
! within the modelling period (nested grids)
!*********************************************************************
do k=1,numbnests
open(unitavailab,file=path(numpath+2*(k-1)+2) &
(1:length(numpath+2*(k-1)+2)),status='old',err=998)
do i=1,3
read(unitavailab,*)
end do
numbwfn(k)=0
700 read(unitavailab,'(i8,1x,i6,2(6x,a255))',end=699) ldat, &
ltim,fname,spec
jul=juldate(ldat,ltim)
if ((jul.ge.beg).and.(jul.le.end)) then
numbwfn(k)=numbwfn(k)+1
if (numbwfn(k).gt.maxwf) then ! check exceedance of dimension
write(*,*) 'Number of nested wind fields is too great.'
write(*,*) 'Reduce modelling period (file "COMMAND") or'
write(*,*) 'reduce number of wind fields (file "AVAILABLE").'
stop
endif
wfname1n(k,numbwfn(k))=fname
wfspec1n(k,numbwfn(k))=spec
wftime1n(k,numbwfn(k))=nint((jul-bdate)*86400._dp)
endif
goto 700 ! next wind field
699 continue
close(unitavailab)
end do
! Check wind field times of file AVAILABLE (expected to be in temporal order)
!****************************************************************************
if (numbwf.eq.0) then
write(*,*) ' #### FLEXPART MODEL ERROR! NO WIND FIELDS #### '
write(*,*) ' #### AVAILABLE FOR SELECTED TIME PERIOD. #### '
stop
endif
do i=2,numbwf
if (wftime1(i).le.wftime1(i-1)) then
write(*,*) 'FLEXPART ERROR: FILE AVAILABLE IS CORRUPT.'
write(*,*) 'THE WIND FIELDS ARE NOT IN TEMPORAL ORDER.'
write(*,*) 'PLEASE CHECK FIELD ',wfname1(i)
stop
endif
end do
! Check wind field times of file AVAILABLE for the nested fields
! (expected to be in temporal order)
!***************************************************************
do k=1,numbnests
if (numbwfn(k).eq.0) then
write(*,*) '#### FLEXPART MODEL ERROR! NO WIND FIELDS ####'
write(*,*) '#### AVAILABLE FOR SELECTED TIME PERIOD. ####'
stop
endif
do i=2,numbwfn(k)
if (wftime1n(k,i).le.wftime1n(k,i-1)) then
write(*,*) 'FLEXPART ERROR: FILE AVAILABLE IS CORRUPT. '
write(*,*) 'THE NESTED WIND FIELDS ARE NOT IN TEMPORAL ORDER.'
write(*,*) 'PLEASE CHECK FIELD ',wfname1n(k,i)
write(*,*) 'AT NESTING LEVEL ',k
stop
endif
end do
end do
! For backward trajectories, reverse the order of the windfields
!***************************************************************
if (ideltas.ge.0) then
do i=1,numbwf
wfname(i)=wfname1(i)
wfspec(i)=wfspec1(i)
wftime(i)=wftime1(i)
end do
do k=1,numbnests
do i=1,numbwfn(k)
wfnamen(k,i)=wfname1n(k,i)
wfspecn(k,i)=wfspec1n(k,i)
wftimen(k,i)=wftime1n(k,i)
end do
end do
else
do i=1,numbwf
wfname(numbwf-i+1)=wfname1(i)
wfspec(numbwf-i+1)=wfspec1(i)
wftime(numbwf-i+1)=wftime1(i)
end do
do k=1,numbnests
do i=1,numbwfn(k)
wfnamen(k,numbwfn(k)-i+1)=wfname1n(k,i)
wfspecn(k,numbwfn(k)-i+1)=wfspec1n(k,i)
wftimen(k,numbwfn(k)-i+1)=wftime1n(k,i)
end do
end do
endif
! Check the time difference between the wind fields. If it is big,
! write a warning message. If it is too big, terminate the trajectory.
!*********************************************************************
do i=2,numbwf
idiff=abs(wftime(i)-wftime(i-1))
if (idiff.gt.idiffmax) then
write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
write(*,*) 'WIND FIELDS IS TOO BIG FOR TRANSPORT CALCULATION.&
&'
write(*,*) 'THEREFORE, TRAJECTORIES HAVE TO BE SKIPPED.'
else if (idiff.gt.idiffnorm) then
write(*,*) 'FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO'
write(*,*) 'WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION'
write(*,*) 'OF SIMULATION QUALITY.'
endif
end do
do k=1,numbnests
if (numbwfn(k).ne.numbwf) then
write(*,*) 'FLEXPART ERROR: THE AVAILABLE FILES FOR THE'
write(*,*) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH'
write(*,*) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. '
write(*,*) 'ERROR AT NEST LEVEL: ',k
stop
endif
do i=1,numbwf
if (wftimen(k,i).ne.wftime(i)) then
write(*,*) 'FLEXPART ERROR: THE AVAILABLE FILES FOR THE'
write(*,*) 'NESTED WIND FIELDS ARE NOT CONSISTENT WITH'
write(*,*) 'THE AVAILABLE FILE OF THE MOTHER DOMAIN. '
write(*,*) 'ERROR AT NEST LEVEL: ',k
stop
endif
end do
end do
! Reset the times of the wind fields that are kept in memory to no time
!**********************************************************************
do i=1,2
memind(i)=i
memtime(i)=999999999
end do
return
998 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
write(*,'(a)') ' '//path(numpath+2*(k-1)+2) &
(1:length(numpath+2*(k-1)+2))
write(*,*) ' #### CANNOT BE OPENED #### '
stop
999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE #### '
write(*,'(a)') ' '//path(4)(1:length(4))
write(*,*) ' #### CANNOT BE OPENED #### '
stop
end subroutine readavailable
|
