1 |
program cvfloat |
2 |
c |
3 |
c======================================================================= |
4 |
c converts binary float trajectories to netCDF |
5 |
c |
6 |
c * must be compiled with a FORTRAN90 compiler and netcdf library |
7 |
c f90 cvfloat.F /usr/local/lib/libnetcdf.a |
8 |
c f90 cvfloat.F /net/ice/ecco/lib/libnetcdf.a (for the ECCO cluster) |
9 |
c * uses namelist data.float |
10 |
c |
11 |
c Arne Biastoch, abiastoch@ucsd.edu, 11/16/2000 |
12 |
c |
13 |
c======================================================================= |
14 |
c |
15 |
integer stdin, stdout, stderr |
16 |
parameter (stdin = 5, stdout = 6, stderr = 6) |
17 |
c |
18 |
character iotext*80, expnam*60, stamp*32 |
19 |
c |
20 |
c variables for filenames |
21 |
integer narg, npart |
22 |
character*6 cpart |
23 |
character*1 split |
24 |
integer factor(6) |
25 |
data factor / 1.e5, 1.e4, 1.e3, 1.e2, 1.e1, 1.e0 / |
26 |
character*(80) dataFName |
27 |
logical exst |
28 |
c |
29 |
c parameter(spval=-1.0e+23) |
30 |
parameter(spval=-999.) |
31 |
c |
32 |
c number of variables per record |
33 |
parameter(imax=10) |
34 |
c |
35 |
integer narg |
36 |
logical preflag |
37 |
c |
38 |
c netCDF ids |
39 |
c |
40 |
integer iret, ncid, VARid |
41 |
integer partdim,Timedim |
42 |
integer partid, Timeid |
43 |
integer xpartid, ypartid, kpartid |
44 |
integer tempid, saltid, uvelid, vvelid, zetaid |
45 |
c |
46 |
c variable shapes, corner and edge lengths |
47 |
c |
48 |
integer dims(4), corner(4), edges(4) |
49 |
c |
50 |
character*1 inumber(4) |
51 |
c |
52 |
c attribute vectors |
53 |
c |
54 |
integer longval(1) |
55 |
real floatval(2) |
56 |
character*1 charval(1) |
57 |
character name*24, unit*16, grid*2 |
58 |
logical usingSphericalPolarGrid |
59 |
c |
60 |
c data variables for NetCDF |
61 |
c |
62 |
real, dimension(:), allocatable :: pnum,time |
63 |
real, dimension(:,:), allocatable :: xpart,ypart,kpart |
64 |
& ,uvel,vvel,temp,salt,zeta |
65 |
double precision, dimension(:), allocatable :: tmp |
66 |
c |
67 |
c namelist contains |
68 |
c |
69 |
data npart /10/ |
70 |
character*50 outname2 |
71 |
character*50 fName, outname |
72 |
data fName / 'float_trajectories' / |
73 |
character*20 startDate |
74 |
data startDate / '01-JAN-1992:12:00:00' / |
75 |
data expnam /'Experimentname not set in data.float'/ |
76 |
data usingSphericalPolarGrid /.true./ |
77 |
namelist /dimensions/ expnam, startDate, usingSphericalPolarGrid |
78 |
namelist /floats/ fName |
79 |
|
80 |
c |
81 |
c in most systems netcdf.inc should be side in /usr/local/include |
82 |
c include '/usr/local/include/netcdf.inc' |
83 |
c include '/users/guests2/ux451985/netcdf/include/netcdf.inc' |
84 |
include '/net/ice/ecco/include/netcdf.inc' |
85 |
c |
86 |
c call GETARG(1,cpart) |
87 |
c npart=0 |
88 |
c do i=1,6 |
89 |
c read(cpart(i:i),'(a1)') split |
90 |
c npart = npart + factor(i)*(ICHAR(split)-48) |
91 |
c enddo |
92 |
c print*, 'npart= ',npart |
93 |
c call GETARG(2,fName) |
94 |
|
95 |
ioun=11 |
96 |
open(ioun,file='data.float',status='old',form='formatted') |
97 |
read (unit=ioun, end=666, nml=dimensions) |
98 |
write (stdout,dimensions) |
99 |
close (ioun) |
100 |
666 continue |
101 |
open(ioun,file='data.float',status='old',form='formatted') |
102 |
read (unit=ioun, end=999, nml=floats) |
103 |
write (stdout,floats) |
104 |
close (ioun) |
105 |
999 continue |
106 |
|
107 |
c |
108 |
c preliminary use: |
109 |
c if floats should be viewed during a current model run the first |
110 |
c line of the file may not be updated correctly, i.e. there might |
111 |
c be more times than stated at the beginning. By giving a flag |
112 |
c only icount-1 timesteps are used |
113 |
c |
114 |
preflag = .false. |
115 |
narg=iargc() |
116 |
if ( narg .gt. 0 ) preflag = .true. |
117 |
c |
118 |
c strip names |
119 |
c |
120 |
IL=ILNBLNK( fName ) |
121 |
|
122 |
c check existent files |
123 |
c |
124 |
iGmax=1 |
125 |
do m=1,100 |
126 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
127 |
& fName(1:IL),'.',iGmax,'.',1,'.data' |
128 |
inquire( file=dataFname, exist=exst ) |
129 |
if (exst) iGmax = iGmax + 1 |
130 |
enddo |
131 |
|
132 |
jGmax=1 |
133 |
do m=1,100 |
134 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
135 |
& fName(1:IL),'.',1,'.',jGmax,'.data' |
136 |
inquire( file=dataFname, exist=exst ) |
137 |
if (exst) jGmax = jGmax + 1 |
138 |
enddo |
139 |
|
140 |
iGmax = iGmax - 1 |
141 |
jGmax = jGmax - 1 |
142 |
print*, 'There are ',iGmax,' x ',jGmax,' files' |
143 |
|
144 |
c open first file and get dimensions (float number and time) |
145 |
c |
146 |
ilen=10*8 |
147 |
allocate (tmp(imax)) |
148 |
c |
149 |
write(dataFname(1:80),'(2a,a)') |
150 |
& fName(1:IL),'.001.001.data' |
151 |
open(1,file=dataFname,status='old',form='unformatted' |
152 |
& ,access='direct',recl=ilen) |
153 |
c |
154 |
read(1,rec=1) tmp |
155 |
rcountstart = SNGL(tmp(2)) |
156 |
rcountdelta = SNGL(tmp(4)) |
157 |
icount = INT(tmp(5)) |
158 |
npart = INT(tmp(6)) |
159 |
close(1) |
160 |
print*, 'npart = ',npart |
161 |
print*, 'timesteps= ',icount |
162 |
if (preflag) then |
163 |
icount=icount-1 |
164 |
print*, 'preliminary --> use timesteps= ',icount |
165 |
endif |
166 |
|
167 |
c----------------------------------------------------------------------- |
168 |
c allocate variables |
169 |
c----------------------------------------------------------------------- |
170 |
c |
171 |
allocate (pnum(npart)) |
172 |
allocate (time(icount)) |
173 |
allocate (xpart(npart,icount)) |
174 |
allocate (ypart(npart,icount)) |
175 |
allocate (kpart(npart,icount)) |
176 |
allocate (uvel(npart,icount)) |
177 |
allocate (vvel(npart,icount)) |
178 |
allocate (temp(npart,icount)) |
179 |
allocate (salt(npart,icount)) |
180 |
allocate (zeta(npart,icount)) |
181 |
|
182 |
c initialize arrays |
183 |
c |
184 |
do m=1,npart |
185 |
do n=1,icount |
186 |
xpart(m,n) = spval |
187 |
ypart(m,n) = spval |
188 |
kpart(m,n) = spval |
189 |
uvel(m,n) = spval |
190 |
vvel(m,n) = spval |
191 |
temp(m,n) = spval |
192 |
salt(m,n) = spval |
193 |
zeta(m,n) = spval |
194 |
enddo |
195 |
enddo |
196 |
|
197 |
c generate axes |
198 |
c |
199 |
time(1)=rcountstart |
200 |
do m=2,icount |
201 |
time(m) = time(m-1)+rcountdelta |
202 |
enddo |
203 |
print*, 'time: ',time |
204 |
c |
205 |
do ip=1,npart |
206 |
pnum(ip) = FLOAT(ip) |
207 |
enddo |
208 |
c print*, 'pnum: ',pnum |
209 |
c |
210 |
c----------------------------------------------------------------------- |
211 |
c open files and read input |
212 |
c----------------------------------------------------------------------- |
213 |
c |
214 |
c |
215 |
itotalrecord = 0 |
216 |
|
217 |
do iG=1,iGmax |
218 |
do jG=1,jGmax |
219 |
c |
220 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
221 |
& fName(1:IL),'.',iG,'.',jG,'.data' |
222 |
open(1,file=dataFname,status='old',form='unformatted' |
223 |
& ,access='direct',recl=ilen) |
224 |
c |
225 |
read(1,rec=1) tmp |
226 |
imaxrecord = INT(tmp(1)) |
227 |
print*,'read ',dataFname |
228 |
itotalrecord = itotalrecord + imaxrecord |
229 |
|
230 |
do irec=2,imaxrecord+1 |
231 |
|
232 |
read(1,rec=irec) tmp |
233 |
ip = INT(tmp(1)) |
234 |
if (ip .gt. npart) then |
235 |
print*,'ip out of order: ',ip,npart |
236 |
stop |
237 |
endif |
238 |
np = INT((tmp(2)-rcountstart)/rcountdelta+1) |
239 |
|
240 |
c this is only for prelimiray results. Use only icount-1 timesteps |
241 |
if (preflag .and. (np .gt. icount .or. np .lt. 1)) |
242 |
& goto 100 |
243 |
|
244 |
if (usingSphericalPolarGrid) then |
245 |
xpart(ip,np) = SNGL(tmp(3)) |
246 |
ypart(ip,np) = SNGL(tmp(4)) |
247 |
else |
248 |
xpart(ip,np) = SNGL(tmp(3))/1000. |
249 |
ypart(ip,np) = SNGL(tmp(4))/1000. |
250 |
endif |
251 |
kpart(ip,np) = SNGL(tmp(5)) |
252 |
uvel(ip,np) = SNGL(tmp(6)) |
253 |
vvel(ip,np) = SNGL(tmp(7)) |
254 |
temp(ip,np) = SNGL(tmp(8)) |
255 |
salt(ip,np) = SNGL(tmp(9)) |
256 |
zeta(ip,np) = SNGL(tmp(10)) |
257 |
c if (ip .eq. 59) |
258 |
c & print*, 'rec= ',irec,' npart= ',ip,' timestep= ',np, |
259 |
c & time(np), |
260 |
c & xpart(ip,np),ypart(ip,np),kpart(ip,np) |
261 |
100 continue |
262 |
enddo |
263 |
|
264 |
close(1) |
265 |
enddo |
266 |
enddo |
267 |
|
268 |
print*,icount,' x ',npart,' = ',icount*npart,' records expected,', |
269 |
& ' found ',itotalrecord,' float records' |
270 |
print*,'==> ',icount*npart-itotalrecord,' float records missing' |
271 |
c |
272 |
c |
273 |
c----------------------------------------------------------------------- |
274 |
c define netCDF variables |
275 |
c----------------------------------------------------------------------- |
276 |
c |
277 |
write(stdout,*) ' Start Converting' |
278 |
c |
279 |
c enter define mode: NCCLOB=overwrite, NCNOCLOB=do not overwrite |
280 |
c |
281 |
IL=ILNBLNK( fname ) |
282 |
outname2=fname(1:IL)//'.cdf' |
283 |
write (stdout,*) |
284 |
& ' ==> Writing a trajectories to file ',outname2(1:IL+4) |
285 |
ncid = nccre(outname2(1:IL+4), NCCLOB, iret) |
286 |
if (iret .ne. 0) write(stdout,*) 'Error: Open NetCDF file' |
287 |
c |
288 |
c define dimensions |
289 |
c |
290 |
partdim = ncddef(ncid, 'Particles', npart, iret) |
291 |
Timedim = ncddef(ncid, 'Time', NCUNLIM, iret) |
292 |
if (iret .ne. 0) write(stdout,*) 'Error: define dimensions' |
293 |
c |
294 |
c define variables |
295 |
c |
296 |
dims(1) = partdim |
297 |
partid = ncvdef (ncid,'Particles',NCFLOAT,1,dims,iret) |
298 |
dims(1) = Timedim |
299 |
Timeid = ncvdef (ncid,'Time', NCFLOAT,1,dims,iret) |
300 |
if (iret .ne. 0) write(stdout,*) 'Error: define axis ids' |
301 |
c |
302 |
dims(1) = partdim |
303 |
dims(2) = Timedim |
304 |
xpartid = ncvdef (ncid,'xpart', NCFLOAT,2,dims,iret) |
305 |
ypartid = ncvdef (ncid,'ypart', NCFLOAT,2,dims,iret) |
306 |
kpartid = ncvdef (ncid,'kpart', NCFLOAT,2,dims,iret) |
307 |
uvelid = ncvdef (ncid,'uvel', NCFLOAT,2,dims,iret) |
308 |
vvelid = ncvdef (ncid,'vvel', NCFLOAT,2,dims,iret) |
309 |
tempid = ncvdef (ncid,'temp', NCFLOAT,2,dims,iret) |
310 |
saltid = ncvdef (ncid,'salt', NCFLOAT,2,dims,iret) |
311 |
zetaid = ncvdef (ncid,'zeta', NCFLOAT,2,dims,iret) |
312 |
if (iret .ne. 0) write(stdout,*) 'Error: define variable ids' |
313 |
c |
314 |
c----------------------------------------------------------------------- |
315 |
c assign attributes |
316 |
c----------------------------------------------------------------------- |
317 |
c |
318 |
charval(1) = ' ' |
319 |
c |
320 |
name = 'Particles Number ' |
321 |
c unit = 'particle number ' |
322 |
call ncaptc(ncid, partid, 'long_name', NCCHAR, 24, name, iret) |
323 |
call ncaptc(ncid, partid, 'units', NCCHAR, 16, unit, iret) |
324 |
c |
325 |
name = 'Time' |
326 |
unit = 'seconds' |
327 |
call ncaptc(ncid, Timeid, 'long_name', NCCHAR, 24, name, iret) |
328 |
call ncaptc(ncid, Timeid, 'units', NCCHAR, 16, unit, iret) |
329 |
call ncaptc(ncid, Timeid,'time_origin',NCCHAR, 20,startDate, iret) |
330 |
if (iret .ne. 0) write(stdout,*) 'Error: assign axis attributes' |
331 |
c |
332 |
floatval(1) = spval |
333 |
c |
334 |
if (usingSphericalPolarGrid) then |
335 |
name = 'LONGITUDE ' |
336 |
unit = 'degrees_W ' |
337 |
else |
338 |
name = 'X_t ' |
339 |
unit = 'kilometer ' |
340 |
endif |
341 |
call ncaptc(ncid, xpartid, 'long_name', NCCHAR, 24, name, iret) |
342 |
call ncaptc(ncid, xpartid, 'units', NCCHAR, 16, unit, iret) |
343 |
call ncapt (ncid, xpartid,'missing_value',NCFLOAT,1,floatval,iret) |
344 |
call ncapt (ncid, xpartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
345 |
c |
346 |
if (usingSphericalPolarGrid) then |
347 |
name = 'LATITUDE ' |
348 |
unit = 'degrees_N ' |
349 |
else |
350 |
name = 'Y_t ' |
351 |
unit = 'kilometer ' |
352 |
endif |
353 |
call ncaptc(ncid, ypartid, 'long_name', NCCHAR, 24, name, iret) |
354 |
call ncaptc(ncid, ypartid, 'units', NCCHAR, 16, unit, iret) |
355 |
call ncapt (ncid, ypartid,'missing_value',NCFLOAT,1,floatval,iret) |
356 |
call ncapt (ncid, ypartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
357 |
c |
358 |
name = 'LEVEL ' |
359 |
unit = 'LEVEL ' |
360 |
call ncaptc(ncid, kpartid, 'long_name', NCCHAR, 24, name, iret) |
361 |
call ncaptc(ncid, kpartid, 'units', NCCHAR, 16, unit, iret) |
362 |
call ncapt (ncid, kpartid,'missing_value',NCFLOAT,1,floatval,iret) |
363 |
call ncapt (ncid, kpartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
364 |
c |
365 |
name = 'POTENTIAL TEMPERATURE ' |
366 |
unit = 'deg C ' |
367 |
call ncaptc(ncid, tempid, 'long_name', NCCHAR, 24, name, iret) |
368 |
call ncaptc(ncid, tempid, 'units', NCCHAR, 16, unit, iret) |
369 |
call ncapt (ncid, tempid, 'missing_value',NCFLOAT,1,floatval,iret) |
370 |
call ncapt (ncid, tempid, '_FillValue', NCFLOAT, 1,floatval, iret) |
371 |
c |
372 |
name = 'SALINITY ' |
373 |
unit = 'PSU ' |
374 |
call ncaptc(ncid, saltid, 'long_name', NCCHAR, 24, name, iret) |
375 |
call ncaptc(ncid, saltid, 'units', NCCHAR, 16, unit, iret) |
376 |
call ncapt (ncid, saltid, 'missing_value',NCFLOAT,1,floatval,iret) |
377 |
call ncapt (ncid, saltid, '_FillValue', NCFLOAT, 1,floatval, iret) |
378 |
c |
379 |
name = 'U VELOCITY ' |
380 |
unit = 'm/sec' |
381 |
call ncaptc(ncid, uvelid, 'long_name', NCCHAR, 24, name, iret) |
382 |
call ncaptc(ncid, uvelid, 'units', NCCHAR, 16, unit, iret) |
383 |
call ncapt (ncid, uvelid, 'missing_value',NCFLOAT,1,floatval,iret) |
384 |
call ncapt (ncid, uvelid, '_FillValue', NCFLOAT, 1,floatval, iret) |
385 |
c |
386 |
name = 'V VELOCITY ' |
387 |
unit = 'm/sec' |
388 |
call ncaptc(ncid, vvelid, 'long_name', NCCHAR, 24, name, iret) |
389 |
call ncaptc(ncid, vvelid, 'units', NCCHAR, 16, unit, iret) |
390 |
call ncapt (ncid, vvelid, 'missing_value',NCFLOAT,1,floatval,iret) |
391 |
call ncapt (ncid, vvelid, '_FillValue', NCFLOAT, 1,floatval, iret) |
392 |
c |
393 |
name = 'SEA SURFACE HEIGHT ' |
394 |
unit = 'm ' |
395 |
call ncaptc(ncid, zetaid, 'long_name', NCCHAR, 24, name, iret) |
396 |
call ncaptc(ncid, zetaid, 'units', NCCHAR, 16, unit, iret) |
397 |
call ncapt (ncid, zetaid,'missing_value',NCFLOAT, 1,floatval,iret) |
398 |
call ncapt (ncid, zetaid,'_FillValue', NCFLOAT, 1, floatval, iret) |
399 |
c |
400 |
if (iret .ne. 0) write(stdout,*) 'Error: define variable attrib.' |
401 |
c |
402 |
expname= ' ' |
403 |
stamp = ' ' |
404 |
call ncaptc(ncid, NCGLOBAL, 'title', NCCHAR, 60, expnam, iret) |
405 |
call ncaptc(ncid, NCGLOBAL, 'history', NCCHAR, 32, stamp, iret) |
406 |
c |
407 |
c----------------------------------------------------------------------- |
408 |
c leave define mode |
409 |
c----------------------------------------------------------------------- |
410 |
c |
411 |
call ncendf(ncid, iret) |
412 |
c |
413 |
c |
414 |
c----------------------------------------------------------------------- |
415 |
c put variables in netCDF file |
416 |
c----------------------------------------------------------------------- |
417 |
c |
418 |
c store Particles |
419 |
corner(1) = 1 |
420 |
edges(1) = npart |
421 |
call ncvpt(ncid, partid, corner, edges, pnum, iret) |
422 |
c |
423 |
c store Time |
424 |
corner(1) = 1 |
425 |
edges(1) = icount |
426 |
call ncvpt(ncid, Timeid, corner, edges, Time, iret) |
427 |
c |
428 |
c store values |
429 |
corner(1) = 1 |
430 |
corner(2) = 1 |
431 |
edges(1) = npart |
432 |
edges(2) = icount |
433 |
VARid=xpartid |
434 |
call ncvpt(ncid, VARid, corner, edges, xpart, iret) |
435 |
VARid=ypartid |
436 |
call ncvpt(ncid, VARid, corner, edges, ypart, iret) |
437 |
VARid=kpartid |
438 |
call ncvpt(ncid, VARid, corner, edges, kpart, iret) |
439 |
VARid=tempid |
440 |
call ncvpt(ncid, VARid, corner, edges, temp, iret) |
441 |
VARid=saltid |
442 |
call ncvpt(ncid, VARid, corner, edges, salt, iret) |
443 |
VARid=uvelid |
444 |
call ncvpt(ncid, VARid, corner, edges, uvel, iret) |
445 |
VARid=vvelid |
446 |
call ncvpt(ncid, VARid, corner, edges, vvel, iret) |
447 |
VARid=zetaid |
448 |
call ncvpt(ncid, VARid, corner, edges, zeta, iret) |
449 |
c |
450 |
if (iret .ne. 0) write(stdout,*) 'Error: write variables' |
451 |
call ncclos (ncid, iret) |
452 |
c |
453 |
write(stdout,*) ' End ' |
454 |
|
455 |
end |
456 |
INTEGER FUNCTION ILNBLNK( string ) |
457 |
C /==========================================================\ |
458 |
C | FUNCTION ILNBLNK | |
459 |
C | o Find last non-blank in character string. | |
460 |
C \==========================================================/ |
461 |
IMPLICIT NONE |
462 |
CHARACTER*(*) string |
463 |
CEndOfInterface |
464 |
INTEGER L, LS |
465 |
C |
466 |
LS = LEN(string) |
467 |
ILNBLNK = LS |
468 |
DO 10 L = LS, 1, -1 |
469 |
IF ( string(L:L) .EQ. ' ' ) GOTO 10 |
470 |
ILNBLNK = L |
471 |
GOTO 11 |
472 |
10 CONTINUE |
473 |
11 CONTINUE |
474 |
C |
475 |
RETURN |
476 |
END |