| 19 |
michaesp |
1 |
PROGRAM lsl2rdf
|
|
|
2 |
|
|
|
3 |
c ***********************************************************************
|
|
|
4 |
c * Convert a trajectory file into a lat/lon/p RDF-netCDF file *
|
|
|
5 |
c * RDF = Reverse Domain Filling *
|
|
|
6 |
c * Michael Sprenger / Spring, summer 2010 *
|
|
|
7 |
c ***********************************************************************
|
|
|
8 |
|
|
|
9 |
use netcdf
|
|
|
10 |
|
|
|
11 |
implicit none
|
|
|
12 |
|
|
|
13 |
c ----------------------------------------------------------------------
|
|
|
14 |
c Declaration of parameters
|
|
|
15 |
c ----------------------------------------------------------------------
|
|
|
16 |
|
|
|
17 |
c Real comparison
|
|
|
18 |
real eps
|
|
|
19 |
parameter (eps=0.001)
|
|
|
20 |
|
|
|
21 |
c ----------------------------------------------------------------------
|
|
|
22 |
c Declaration of variables
|
|
|
23 |
c ----------------------------------------------------------------------
|
|
|
24 |
|
|
|
25 |
c Input and output format for trajectories (see iotra.f)
|
|
|
26 |
character*80 inpfile ! Input filename
|
|
|
27 |
character*80 outfile ! Output filename
|
|
|
28 |
|
|
|
29 |
c Trajectories
|
|
|
30 |
integer ntra ! Number of trajectories
|
|
|
31 |
integer ntim ! Number of times
|
|
|
32 |
integer ncol ! Number of columns
|
|
|
33 |
real,allocatable, dimension (:,:,:) :: tra ! Trajectories (ntra,ntim,ncol)
|
|
|
34 |
character*80 vars(100) ! Variable names
|
|
|
35 |
integer refdate(6) ! Reference date
|
|
|
36 |
|
|
|
37 |
c Output RDF-netCDF file
|
|
|
38 |
character*80 cdfname
|
|
|
39 |
character*80 varname
|
|
|
40 |
character*80 longname
|
|
|
41 |
character*80 unit
|
|
|
42 |
integer nx,ny,nz
|
|
|
43 |
real,allocatable, dimension (:) :: lon,lat,lev
|
|
|
44 |
real,allocatable, dimension (:,:,:) :: arr
|
|
|
45 |
real time
|
|
|
46 |
integer crefile,crevar,cretime
|
|
|
47 |
integer,allocatable, dimension (:,:,:) :: ind
|
|
|
48 |
integer indx,indy,indz
|
|
|
49 |
|
|
|
50 |
c Auxiliary variables
|
|
|
51 |
integer inpmode
|
|
|
52 |
integer stat
|
|
|
53 |
integer fid
|
|
|
54 |
integer i,j,k,ivar,itim
|
|
|
55 |
real tmp(10000)
|
|
|
56 |
integer ok
|
|
|
57 |
integer count
|
|
|
58 |
|
|
|
59 |
c ----------------------------------------------------------------------
|
|
|
60 |
c Read input data
|
|
|
61 |
c ----------------------------------------------------------------------
|
|
|
62 |
|
|
|
63 |
c Write start message
|
|
|
64 |
print*,'========================================================='
|
|
|
65 |
print*,' *** START OF PROGRAM LSL2RDF ***'
|
|
|
66 |
print*
|
|
|
67 |
|
|
|
68 |
c Read parameters
|
|
|
69 |
open(10,file='lsl2rdf.param')
|
|
|
70 |
read(10,*) inpfile
|
|
|
71 |
read(10,*) outfile
|
|
|
72 |
read(10,*) ntra,ntim,ncol
|
|
|
73 |
close(10)
|
|
|
74 |
|
|
|
75 |
c Determine the formats
|
|
|
76 |
call mode_tra(inpmode,inpfile)
|
|
|
77 |
if (inpmode.eq.-1) inpmode=1
|
|
|
78 |
|
|
|
79 |
c Allocate memory for input trajectory
|
|
|
80 |
allocate(tra(ntra,ntim,ncol),stat=stat)
|
|
|
81 |
if (stat.ne.0) print*,'*** error allocating array tra ***'
|
|
|
82 |
|
|
|
83 |
c Read inpufile
|
|
|
84 |
call ropen_tra(fid,inpfile,ntra,ntim,ncol,refdate,vars,inpmode)
|
|
|
85 |
call read_tra (fid,tra,ntra,ntim,ncol,inpmode)
|
|
|
86 |
call close_tra(fid,inpmode)
|
|
|
87 |
|
|
|
88 |
c ----------------------------------------------------------------------
|
|
|
89 |
c Determine output grid and get mapping {output array} <-> {trajectory}
|
|
|
90 |
c ----------------------------------------------------------------------
|
|
|
91 |
|
|
|
92 |
c Get a list of longitudes at time 0
|
|
|
93 |
nx = 0
|
|
|
94 |
do i=1,ntra
|
|
|
95 |
ok = 1
|
|
|
96 |
do j=1,nx
|
|
|
97 |
if ( abs(tmp(j)-tra(i,1,2)).lt.eps ) then
|
|
|
98 |
ok = 0
|
|
|
99 |
goto 100
|
|
|
100 |
endif
|
|
|
101 |
enddo
|
|
|
102 |
100 if ( ok.eq.1 ) then
|
|
|
103 |
nx = nx + 1
|
|
|
104 |
tmp(nx) = tra(i,1,2)
|
|
|
105 |
endif
|
|
|
106 |
enddo
|
|
|
107 |
|
|
|
108 |
if ( nx.gt.1 ) then
|
|
|
109 |
call sort(nx,tmp)
|
|
|
110 |
endif
|
|
|
111 |
|
|
|
112 |
allocate(lon(nx),stat=stat)
|
|
|
113 |
if (stat.ne.0) print*,'*** error allocating array lon ***'
|
|
|
114 |
do i=1,nx
|
|
|
115 |
lon(i) = tmp(i)
|
|
|
116 |
enddo
|
|
|
117 |
|
|
|
118 |
c Get a list of latitudes at time 0
|
|
|
119 |
ny = 0
|
|
|
120 |
do i=1,ntra
|
|
|
121 |
ok = 1
|
|
|
122 |
do j=1,ny
|
|
|
123 |
if ( abs(tmp(j)-tra(i,1,3)).lt.eps ) then
|
|
|
124 |
ok = 0
|
|
|
125 |
endif
|
|
|
126 |
enddo
|
|
|
127 |
101 if ( ok.eq.1 ) then
|
|
|
128 |
ny = ny + 1
|
|
|
129 |
tmp(ny) = tra(i,1,3)
|
|
|
130 |
endif
|
|
|
131 |
enddo
|
|
|
132 |
|
|
|
133 |
if ( ny.gt.1 ) then
|
|
|
134 |
call sort(ny,tmp)
|
|
|
135 |
endif
|
|
|
136 |
|
|
|
137 |
allocate(lat(ny),stat=stat)
|
|
|
138 |
if (stat.ne.0) print*,'*** error allocating array lat ***'
|
|
|
139 |
do i=1,ny
|
|
|
140 |
lat(i) = tmp(i)
|
|
|
141 |
enddo
|
|
|
142 |
|
|
|
143 |
c Get a list of pressure levels at time 0
|
|
|
144 |
nz = 0
|
|
|
145 |
do i=1,ntra
|
|
|
146 |
ok = 1
|
|
|
147 |
do j=1,nz
|
|
|
148 |
if ( abs(tmp(j)-tra(i,1,4)).lt.eps ) then
|
|
|
149 |
ok = 0
|
|
|
150 |
goto 102
|
|
|
151 |
endif
|
|
|
152 |
enddo
|
|
|
153 |
102 if ( ok.eq.1 ) then
|
|
|
154 |
nz = nz + 1
|
|
|
155 |
tmp(nz) = tra(i,1,4)
|
|
|
156 |
endif
|
|
|
157 |
enddo
|
|
|
158 |
|
|
|
159 |
if ( nz.gt.1 ) then
|
|
|
160 |
call sort(nz,tmp)
|
|
|
161 |
endif
|
|
|
162 |
|
|
|
163 |
allocate(lev(nz),stat=stat)
|
|
|
164 |
if (stat.ne.0) print*,'*** error allocating array lev ***'
|
|
|
165 |
do i=1,nz
|
|
|
166 |
lev(i) = tmp(i)
|
|
|
167 |
enddo
|
|
|
168 |
|
|
|
169 |
c Write output grid
|
|
|
170 |
print*,'---- OUTPUT GRID ----------------------------------------'
|
|
|
171 |
print*
|
|
|
172 |
write(*,'(a10,1000f5.0)') 'lon:',(lon(i),i=1,nx)
|
|
|
173 |
write(*,'(a10,1000f5.0)') 'lat:',(lat(i),i=1,ny)
|
|
|
174 |
write(*,'(a10,1000f5.0)') 'lev:',(lev(i),i=1,nz)
|
|
|
175 |
print*
|
|
|
176 |
|
|
|
177 |
c Decide whether several levels are really needed
|
|
|
178 |
if ( (nx*ny).eq.ntra ) then
|
|
|
179 |
if ( nz.ne.1 ) then
|
|
|
180 |
print*
|
|
|
181 |
print*,'All levels will be mapped to index level 1'
|
|
|
182 |
nz = 1
|
|
|
183 |
endif
|
|
|
184 |
endif
|
|
|
185 |
|
|
|
186 |
c Allocate memory for output RDF-netCDF
|
|
|
187 |
allocate(arr(nx,ny,nz),stat=stat)
|
|
|
188 |
if (stat.ne.0) print*,'*** error allocating array arr ***'
|
|
|
189 |
allocate(ind(nx,ny,nz),stat=stat)
|
|
|
190 |
if (stat.ne.0) print*,'*** error allocating array ind ***'
|
|
|
191 |
|
|
|
192 |
c Init the mapping array <-> trajectory
|
|
|
193 |
do i=1,nx
|
|
|
194 |
do j=1,ny
|
|
|
195 |
do k=1,nz
|
|
|
196 |
ind(i,j,k) = 0
|
|
|
197 |
enddo
|
|
|
198 |
enddo
|
|
|
199 |
enddo
|
|
|
200 |
|
|
|
201 |
c Now get for each grid point the right trajectory index
|
|
|
202 |
do i=1,ntra
|
|
|
203 |
|
|
|
204 |
indx = 0
|
|
|
205 |
do j=1,nx
|
|
|
206 |
if ( abs(lon(j)-tra(i,1,2)).lt.eps ) then
|
|
|
207 |
indx = j
|
|
|
208 |
goto 200
|
|
|
209 |
endif
|
|
|
210 |
enddo
|
|
|
211 |
200 continue
|
|
|
212 |
|
|
|
213 |
indy = 0
|
|
|
214 |
do j=1,ny
|
|
|
215 |
if ( abs(lat(j)-tra(i,1,3)).lt.eps ) then
|
|
|
216 |
indy = j
|
|
|
217 |
goto 201
|
|
|
218 |
endif
|
|
|
219 |
enddo
|
|
|
220 |
201 continue
|
|
|
221 |
|
|
|
222 |
if ( nz.gt.1 ) then
|
|
|
223 |
indz = 0
|
|
|
224 |
do j=1,nz
|
|
|
225 |
if ( abs(lev(j)-tra(i,1,4)).lt.eps ) then
|
|
|
226 |
indz = j
|
|
|
227 |
goto 202
|
|
|
228 |
endif
|
|
|
229 |
enddo
|
|
|
230 |
202 continue
|
|
|
231 |
else
|
|
|
232 |
indz = 1
|
|
|
233 |
endif
|
|
|
234 |
|
|
|
235 |
if ( (indx.ne.0).and.(indy.ne.0).and.(indz.ne.0) ) then
|
|
|
236 |
ind(indx,indy,indz) = i
|
|
|
237 |
endif
|
|
|
238 |
|
|
|
239 |
enddo
|
|
|
240 |
|
|
|
241 |
c Check whether all grid points are linked to a traejctory
|
|
|
242 |
count = 0
|
|
|
243 |
do i=1,nx
|
|
|
244 |
do j=1,ny
|
|
|
245 |
do k=1,nz
|
|
|
246 |
if ( ind(i,j,k).eq.0 ) count = count + 1
|
|
|
247 |
enddo
|
|
|
248 |
enddo
|
|
|
249 |
enddo
|
|
|
250 |
if ( count.ne.0 ) then
|
|
|
251 |
print*,' WARNING: not all grid points in RDF domain linked'
|
|
|
252 |
print*,' to a trajectory... proceed anyway.'
|
|
|
253 |
endif
|
|
|
254 |
|
|
|
255 |
|
|
|
256 |
c ----------------------------------------------------------------------
|
|
|
257 |
c Fill the output array
|
|
|
258 |
c ----------------------------------------------------------------------
|
|
|
259 |
|
|
|
260 |
print*,'---- MAPPING TO RDF GRID --------------------------------'
|
|
|
261 |
print*
|
|
|
262 |
|
|
|
263 |
c Loop over all variables
|
|
|
264 |
crefile = 1
|
|
|
265 |
do ivar=2,ncol
|
|
|
266 |
|
|
|
267 |
c Loop over all times
|
|
|
268 |
cretime = 1
|
|
|
269 |
crevar = 1
|
|
|
270 |
do itim=1,ntim
|
|
|
271 |
|
|
|
272 |
c Get the time from the first trajectory
|
|
|
273 |
time = tra(1,itim,1)
|
|
|
274 |
write(*,'(a20,a10,f8.2)')
|
|
|
275 |
> ' variable,time : ',trim( vars(ivar) ),time
|
|
|
276 |
|
|
|
277 |
c Do the remapping
|
|
|
278 |
do i=1,nx
|
|
|
279 |
do j=1,ny
|
|
|
280 |
do k=1,nz
|
|
|
281 |
if ( ind(i,j,k).ne.0 ) then
|
|
|
282 |
arr(i,j,k) = tra( ind(i,j,k), itim, ivar )
|
|
|
283 |
endif
|
|
|
284 |
enddo
|
|
|
285 |
enddo
|
|
|
286 |
enddo
|
|
|
287 |
|
|
|
288 |
c Save the array to RDF-netCDF
|
|
|
289 |
call writecdf2D_cf
|
|
|
290 |
> (outfile,vars(ivar),nx,lon,ny,lat,nz,lev,
|
|
|
291 |
> arr,time,crefile,crevar,cretime)
|
|
|
292 |
|
|
|
293 |
crefile = 0
|
|
|
294 |
crevar = 0
|
|
|
295 |
|
|
|
296 |
enddo
|
|
|
297 |
|
|
|
298 |
enddo
|
|
|
299 |
|
|
|
300 |
c Write end message
|
|
|
301 |
print*
|
|
|
302 |
print*,' *** END OF PROGRAM LSL2RDF ***'
|
|
|
303 |
print*,'========================================================='
|
|
|
304 |
|
|
|
305 |
end
|
|
|
306 |
|
|
|
307 |
c ********************************************************************
|
|
|
308 |
c * SUBROUTINES *
|
|
|
309 |
c ********************************************************************
|
|
|
310 |
|
|
|
311 |
c --------------------------------------------------------------------
|
|
|
312 |
c Subroutines to write the CF netcdf output file
|
|
|
313 |
c --------------------------------------------------------------------
|
|
|
314 |
|
|
|
315 |
subroutine writecdf2D_cf
|
|
|
316 |
> (cdfname,varname,nx,lon,ny,lat,nz,lev,
|
|
|
317 |
< arr,time,crefile,crevar,cretime)
|
|
|
318 |
|
|
|
319 |
USE netcdf
|
|
|
320 |
|
|
|
321 |
IMPLICIT NONE
|
|
|
322 |
|
|
|
323 |
c Declaration of input parameters
|
|
|
324 |
character*80 cdfname
|
|
|
325 |
character*80 varname,longname,unit
|
|
|
326 |
integer nx,ny,nz
|
|
|
327 |
real lon(nx),lat(ny),lev(nz)
|
|
|
328 |
real arr(nx,ny,nz)
|
|
|
329 |
real time
|
|
|
330 |
integer crefile,crevar,cretime
|
|
|
331 |
|
|
|
332 |
c Local variables
|
|
|
333 |
integer ierr
|
|
|
334 |
integer ncID
|
|
|
335 |
integer LonDimId, varLonID
|
|
|
336 |
integer LatDimID, varLatID
|
|
|
337 |
integer LevDimID, varLevID
|
|
|
338 |
integer TimeDimID, varTimeID
|
|
|
339 |
real timeindex
|
|
|
340 |
integer i
|
|
|
341 |
integer nvars,varids(100)
|
|
|
342 |
integer ndims,dimids(100)
|
|
|
343 |
real timelist(1000)
|
|
|
344 |
integer ntimes
|
|
|
345 |
integer ind
|
|
|
346 |
integer varID
|
|
|
347 |
real xmin,xmax,ymin,ymax,zmin,zmax
|
|
|
348 |
|
|
|
349 |
c Quick an dirty solution for fieldname conflict
|
|
|
350 |
if ( varname.eq.'time' ) varname = 'TIME'
|
|
|
351 |
|
|
|
352 |
c Initially set error to indicate no errors.
|
|
|
353 |
ierr = 0
|
|
|
354 |
|
|
|
355 |
c ---- Create the netCDF - skip if <crefile=0> ----------------------
|
|
|
356 |
if ( crefile.ne.1 ) goto 100
|
|
|
357 |
|
|
|
358 |
c Create the file
|
|
|
359 |
ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
|
|
|
360 |
|
|
|
361 |
c Define dimensions
|
|
|
362 |
ierr=nf90_def_dim(ncID,'longitude' ,nx , LonDimID )
|
|
|
363 |
ierr=nf90_def_dim(ncID,'latitude' ,ny , LatDimID )
|
|
|
364 |
ierr=nf90_def_dim(ncID,'level ' ,nz , LevDimID )
|
|
|
365 |
ierr=nf90_def_dim(ncID,'time' ,nf90_unlimited, TimeDimID)
|
|
|
366 |
|
|
|
367 |
c Define coordinate Variables
|
|
|
368 |
ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
|
|
|
369 |
> (/ LonDimID /),varLonID)
|
|
|
370 |
ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
|
|
|
371 |
ierr = nf90_put_att(ncID, varLonID, "units" ,"degree_east")
|
|
|
372 |
|
|
|
373 |
ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
|
|
|
374 |
> (/ LatDimID /),varLatID)
|
|
|
375 |
ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
|
|
|
376 |
ierr = nf90_put_att(ncID, varLatID, "units" ,"degree_north")
|
|
|
377 |
|
|
|
378 |
ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
|
|
|
379 |
> (/ LevDimID /),varLevID)
|
|
|
380 |
ierr = nf90_put_att(ncID, varLevID, "standard_name", "level")
|
|
|
381 |
ierr = nf90_put_att(ncID, varLevID, "units" ,"nil")
|
|
|
382 |
|
|
|
383 |
ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
|
|
|
384 |
> (/ TimeDimID /), varTimeID)
|
|
|
385 |
|
|
|
386 |
c Write general global attributes
|
|
|
387 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
|
|
|
388 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',
|
|
|
389 |
> 'RDF trajectory')
|
|
|
390 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source',
|
|
|
391 |
> 'Lagranto Trajectories')
|
|
|
392 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution',
|
|
|
393 |
> 'ETH Zurich, IACETH')
|
|
|
394 |
|
|
|
395 |
c Write grid extension
|
|
|
396 |
xmin = lon(1)
|
|
|
397 |
xmax = lon(1)
|
|
|
398 |
do i=2,nx
|
|
|
399 |
if ( lon(i).lt.xmin ) xmin = lon(i)
|
|
|
400 |
if ( lon(i).gt.xmax ) xmax = lon(i)
|
|
|
401 |
enddo
|
|
|
402 |
|
|
|
403 |
ymin = lat(1)
|
|
|
404 |
ymax = lat(1)
|
|
|
405 |
do i=2,ny
|
|
|
406 |
if ( lat(i).lt.ymin ) ymin = lat(i)
|
|
|
407 |
if ( lat(i).gt.ymax ) ymax = lat(i)
|
|
|
408 |
enddo
|
|
|
409 |
|
|
|
410 |
zmin = lev(1)
|
|
|
411 |
zmax = lev(1)
|
|
|
412 |
do i=2,nz
|
|
|
413 |
if ( lev(i).lt.zmin ) zmin = lev(i)
|
|
|
414 |
if ( lev(i).gt.zmax ) zmax = lev(i)
|
|
|
415 |
enddo
|
|
|
416 |
|
|
|
417 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
|
|
|
418 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
|
|
|
419 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nz',nz )
|
|
|
420 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
|
|
|
421 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmax',xmax )
|
|
|
422 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
|
|
|
423 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymax',ymax )
|
|
|
424 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'zmin',zmin )
|
|
|
425 |
ierr = nf90_put_att(ncID, NF90_GLOBAL, 'zmax',zmax )
|
|
|
426 |
|
|
|
427 |
c Check whether the definition was successful
|
|
|
428 |
ierr = nf90_enddef(ncID)
|
|
|
429 |
if (ierr.gt.0) then
|
|
|
430 |
print*, 'An error occurred while attempting to ',
|
|
|
431 |
> 'finish definition mode.'
|
|
|
432 |
stop
|
|
|
433 |
endif
|
|
|
434 |
|
|
|
435 |
c Write coordinate data
|
|
|
436 |
ierr = nf90_put_var(ncID,varLonID ,lon )
|
|
|
437 |
ierr = nf90_put_var(ncID,varLatID ,lat )
|
|
|
438 |
ierr = nf90_put_var(ncID,varLevID ,lev )
|
|
|
439 |
|
|
|
440 |
c Close netCDF file
|
|
|
441 |
ierr = nf90_close(ncID)
|
|
|
442 |
|
|
|
443 |
100 continue
|
|
|
444 |
|
|
|
445 |
c ---- Define a new variable - skip if <crevar=0> -----------------------
|
|
|
446 |
|
|
|
447 |
if ( crevar.ne.1 ) goto 110
|
|
|
448 |
|
|
|
449 |
c Open the file for read(write access
|
|
|
450 |
ierr = nf90_open (trim(cdfname), NF90_WRITE , ncID)
|
|
|
451 |
|
|
|
452 |
c Get the IDs for dimensions
|
|
|
453 |
ierr = nf90_inq_dimid(ncID,'longitude' , LonDimID )
|
|
|
454 |
ierr = nf90_inq_dimid(ncID,'latitude' , LatDimID )
|
|
|
455 |
ierr = nf90_inq_dimid(ncID,'level' , LevDimID )
|
|
|
456 |
ierr = nf90_inq_dimid(ncID,'time' , TimeDimID)
|
|
|
457 |
|
|
|
458 |
c Enter define mode
|
|
|
459 |
ierr = nf90_redef(ncID)
|
|
|
460 |
|
|
|
461 |
c Write definition and add attributes
|
|
|
462 |
ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
|
|
|
463 |
> (/ LonDimID, LatDimID, LevDimID, TimeDimID /),varID)
|
|
|
464 |
ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99 )
|
|
|
465 |
|
|
|
466 |
c Check whether definition was successful
|
|
|
467 |
ierr = nf90_enddef(ncID)
|
|
|
468 |
if (ierr.gt.0) then
|
|
|
469 |
print*, 'An error occurred while attempting to ',
|
|
|
470 |
> 'finish definition mode.'
|
|
|
471 |
stop
|
|
|
472 |
endif
|
|
|
473 |
|
|
|
474 |
c Close netCDF file
|
|
|
475 |
ierr = nf90_close(ncID)
|
|
|
476 |
|
|
|
477 |
110 continue
|
|
|
478 |
|
|
|
479 |
c ---- Create a new time (unlimited dimension) - skip if <cretime=0> ------
|
|
|
480 |
|
|
|
481 |
if ( cretime.ne.1 ) goto 120
|
|
|
482 |
|
|
|
483 |
c Open the file for read/write access
|
|
|
484 |
ierr = nf90_open (trim(cdfname), NF90_WRITE, ncID)
|
|
|
485 |
|
|
|
486 |
c Get the list of times on the netCDF file
|
|
|
487 |
ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
|
|
|
488 |
if ( ierr.ne.0 ) then
|
|
|
489 |
print*,'Time dimension is not defined on ',trim(cdfname),
|
|
|
490 |
> ' .... Stop'
|
|
|
491 |
stop
|
|
|
492 |
endif
|
|
|
493 |
ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
|
|
|
494 |
ierr = nf90_inq_varid(ncID,'time', varTimeID)
|
|
|
495 |
if ( ierr.ne.0 ) then
|
|
|
496 |
print*,'Variable time is not defined on ',trim(cdfname),
|
|
|
497 |
> ' ... Stop'
|
|
|
498 |
stop
|
|
|
499 |
endif
|
|
|
500 |
ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
501 |
|
|
|
502 |
c Decide whether a new time must be written
|
|
|
503 |
ind = 0
|
|
|
504 |
do i=1,ntimes
|
|
|
505 |
if ( time.eq.timelist(i) ) ind = i
|
|
|
506 |
enddo
|
|
|
507 |
|
|
|
508 |
c Extend the time list if required
|
|
|
509 |
if ( ind.eq.0 ) then
|
|
|
510 |
ntimes = ntimes + 1
|
|
|
511 |
timelist(ntimes) = time
|
|
|
512 |
ierr = nf90_put_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
513 |
endif
|
|
|
514 |
|
|
|
515 |
c Close netCDF file
|
|
|
516 |
ierr = nf90_close(ncID)
|
|
|
517 |
|
|
|
518 |
120 continue
|
|
|
519 |
|
|
|
520 |
c ---- Write data --------------------------------------------------
|
|
|
521 |
|
|
|
522 |
c Open the file for read/write access
|
|
|
523 |
ierr = nf90_open (trim(cdfname), NF90_WRITE , ncID)
|
|
|
524 |
|
|
|
525 |
c Get the varID
|
|
|
526 |
ierr = nf90_inq_varid(ncID,varname, varID )
|
|
|
527 |
if (ierr.ne.0) then
|
|
|
528 |
print*,'Variable ',trim(varname),' is not defined on ',
|
|
|
529 |
> trim(cdfname)
|
|
|
530 |
stop
|
|
|
531 |
endif
|
|
|
532 |
|
|
|
533 |
c Get the time index
|
|
|
534 |
ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
|
|
|
535 |
if ( ierr.ne.0 ) then
|
|
|
536 |
print*,'Time dimension is not defined on ',trim(cdfname),
|
|
|
537 |
> ' .... Stop'
|
|
|
538 |
stop
|
|
|
539 |
endif
|
|
|
540 |
ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
|
|
|
541 |
ierr = nf90_inq_varid(ncID,'time', varTimeID)
|
|
|
542 |
if ( ierr.ne.0 ) then
|
|
|
543 |
print*,'Variable time is not defined on ',trim(cdfname),
|
|
|
544 |
> ' ... Stop'
|
|
|
545 |
stop
|
|
|
546 |
endif
|
|
|
547 |
ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
|
|
|
548 |
ind = 0
|
|
|
549 |
do i=1,ntimes
|
|
|
550 |
if ( time.eq.timelist(i) ) ind = i
|
|
|
551 |
enddo
|
|
|
552 |
if (ind.eq.0) then
|
|
|
553 |
print*,'Time',time,' is not defined on the netCDF file',
|
|
|
554 |
> trim(cdfname),' ... Stop'
|
|
|
555 |
stop
|
|
|
556 |
endif
|
|
|
557 |
|
|
|
558 |
c Write data block
|
|
|
559 |
ierr = nf90_put_var(ncID,varID,arr,
|
|
|
560 |
> start = (/ 1, 1, 1, ind /),
|
|
|
561 |
> count = (/ nx, ny, nz, 1 /) )
|
|
|
562 |
|
|
|
563 |
c Check whether writing was successful
|
|
|
564 |
ierr = nf90_close(ncID)
|
|
|
565 |
if (ierr.ne.0) then
|
|
|
566 |
write(*,*) trim(nf90_strerror(ierr))
|
|
|
567 |
write(*,*) 'An error occurred while attempting to ',
|
|
|
568 |
> 'close the netcdf file.'
|
|
|
569 |
write(*,*) 'in clscdf_CF'
|
|
|
570 |
endif
|
|
|
571 |
|
|
|
572 |
end
|
|
|
573 |
|
|
|
574 |
c --------------------------------------------------------------------
|
|
|
575 |
c Sort an array (from Muerical Recipes)
|
|
|
576 |
c --------------------------------------------------------------------
|
|
|
577 |
|
|
|
578 |
SUBROUTINE SORT(N,RA)
|
|
|
579 |
DIMENSION RA(N)
|
|
|
580 |
L=N/2+1
|
|
|
581 |
IR=N
|
|
|
582 |
10 CONTINUE
|
|
|
583 |
IF(L.GT.1)THEN
|
|
|
584 |
L=L-1
|
|
|
585 |
RRA=RA(L)
|
|
|
586 |
ELSE
|
|
|
587 |
RRA=RA(IR)
|
|
|
588 |
RA(IR)=RA(1)
|
|
|
589 |
IR=IR-1
|
|
|
590 |
IF(IR.EQ.1)THEN
|
|
|
591 |
RA(1)=RRA
|
|
|
592 |
RETURN
|
|
|
593 |
ENDIF
|
|
|
594 |
ENDIF
|
|
|
595 |
I=L
|
|
|
596 |
J=L+L
|
|
|
597 |
20 IF(J.LE.IR)THEN
|
|
|
598 |
IF(J.LT.IR)THEN
|
|
|
599 |
IF(RA(J).LT.RA(J+1))J=J+1
|
|
|
600 |
ENDIF
|
|
|
601 |
IF(RRA.LT.RA(J))THEN
|
|
|
602 |
RA(I)=RA(J)
|
|
|
603 |
I=J
|
|
|
604 |
J=J+J
|
|
|
605 |
ELSE
|
|
|
606 |
J=IR+1
|
|
|
607 |
ENDIF
|
|
|
608 |
GO TO 20
|
|
|
609 |
ENDIF
|
|
|
610 |
RA(I)=RRA
|
|
|
611 |
GO TO 10
|
|
|
612 |
END
|
|
|
613 |
|
|
|
614 |
|
|
|
615 |
|