1 |
program cvprofiles |
2 |
c |
3 |
c======================================================================= |
4 |
c converts binary float profiles to netCDF |
5 |
c |
6 |
c * must be compiled with a FORTRAN90 compiler and netcdf library |
7 |
c f90 cvprofiles.F /usr/local/lib/libnetcdf.a |
8 |
c f90 cvprofiles.F /net/ice/ecco/lib/libnetcdf.a (for the ECCO cluster) |
9 |
c * uses namelist data.profiles |
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 |
parameter (maxcoord=3000, maxdepth=100) |
19 |
character iotext*80, expnam*60, stamp*32 |
20 |
c |
21 |
c variables for filenames |
22 |
integer narg, npart |
23 |
character*6 cpart |
24 |
character*1 split |
25 |
integer factor(6) |
26 |
data factor / 1.e5, 1.e4, 1.e3, 1.e2, 1.e1, 1.e0 / |
27 |
character*(80) dataFName |
28 |
logical exst |
29 |
c |
30 |
c parameter(spval=-1.0e+23) |
31 |
parameter(spval=-999.) |
32 |
|
33 |
c number of variables per record |
34 |
c parameter(imax=10) |
35 |
integer narg |
36 |
logical preflag |
37 |
c |
38 |
c netCDF ids |
39 |
c |
40 |
integer iret, ncid, VARid |
41 |
integer partdim,Timedim,Dep_tdim, Dep_wdim, Dep_wm1dim |
42 |
integer partid, Timeid |
43 |
integer xpartid, ypartid, kpartid |
44 |
integer tempid, saltid, uvelid, vvelid, zetaid |
45 |
integer Dep_tid, Dep_wid, Dep_wm1id |
46 |
c |
47 |
c variable shapes, corner and edge lengths |
48 |
c |
49 |
integer dims(4), corner(4), edges(4) |
50 |
c |
51 |
character*1 inumber(4) |
52 |
c |
53 |
c attribute vectors |
54 |
c |
55 |
integer longval(1) |
56 |
real floatval(2) |
57 |
character*1 charval(1) |
58 |
character name*24, unit*16, grid*2 |
59 |
logical usingSphericalPolarGrid |
60 |
logical iedge |
61 |
c |
62 |
c data variables for NetCDF |
63 |
c |
64 |
real, dimension(:), allocatable :: Dep_t, Dep_w, Dep_wm1 |
65 |
real, dimension(:), allocatable :: pnum,time |
66 |
real, dimension(:,:), allocatable :: xpart,ypart,kpart,zeta |
67 |
real, dimension(:,:,:), allocatable :: uvel,vvel,temp,salt |
68 |
double precision, dimension(:), allocatable :: tmp |
69 |
c |
70 |
c these variables cannot be allocatable because they appear in the namelist |
71 |
c |
72 |
real delZ(maxdepth) |
73 |
c |
74 |
c namelist contains |
75 |
c |
76 |
data npart /10/ |
77 |
character*50 outname2 |
78 |
character*50 fName, outname |
79 |
data fName / 'float_profiles' / |
80 |
character*20 startDate |
81 |
data startDate / '01-JAN-1992:12:00:00' / |
82 |
data expnam /'Experimentname not set in data.profiles'/ |
83 |
data usingSphericalPolarGrid /.true./ |
84 |
namelist /dimensions/ expnam, startDate, usingSphericalPolarGrid |
85 |
namelist /floats/ fName |
86 |
namelist /coord/ Nr, delZ |
87 |
c |
88 |
c in most systems netcdf.inc should be side in /usr/local/include |
89 |
c include '/usr/local/include/netcdf.inc' |
90 |
c include '/users/guests2/ux451985/netcdf/include/netcdf.inc' |
91 |
include '/net/ice/ecco/include/netcdf.inc' |
92 |
|
93 |
ioun=11 |
94 |
open(ioun,file='data.profiles',status='old',form='formatted') |
95 |
read (unit=ioun, end=666, nml=dimensions) |
96 |
write (stdout,dimensions) |
97 |
close (ioun) |
98 |
666 continue |
99 |
open(ioun,file='data.profiles',status='old',form='formatted') |
100 |
read (unit=ioun, end=777, nml=coord) |
101 |
c write (stdout,coord) |
102 |
close (ioun) |
103 |
777 continue |
104 |
open(ioun,file='data.profiles',status='old',form='formatted') |
105 |
read (unit=ioun, end=999, nml=floats) |
106 |
write (stdout,floats) |
107 |
close (ioun) |
108 |
999 continue |
109 |
|
110 |
c |
111 |
c big data set: |
112 |
c if the data set contains a big number of particles and timesteps |
113 |
c it has to be read in chunks. This takes longer but fits better |
114 |
c into the memory. The argument preflag is used to indicate a big |
115 |
c data set. |
116 |
c |
117 |
preflag = .false. |
118 |
narg=iargc() |
119 |
if ( narg .gt. 0 ) preflag = .true. |
120 |
|
121 |
c |
122 |
c strip names |
123 |
c |
124 |
IL=ILNBLNK( fName ) |
125 |
|
126 |
c check existent files |
127 |
c |
128 |
iGmax=1 |
129 |
do m=1,100 |
130 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
131 |
& fName(1:IL),'.',iGmax,'.',1,'.data' |
132 |
inquire( file=dataFname, exist=exst ) |
133 |
if (exst) iGmax = iGmax + 1 |
134 |
enddo |
135 |
|
136 |
jGmax=1 |
137 |
do m=1,100 |
138 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
139 |
& fName(1:IL),'.',1,'.',jGmax,'.data' |
140 |
inquire( file=dataFname, exist=exst ) |
141 |
if (exst) jGmax = jGmax + 1 |
142 |
enddo |
143 |
|
144 |
iGmax = iGmax - 1 |
145 |
jGmax = jGmax - 1 |
146 |
print*, 'There are ',iGmax,' x ',jGmax,' files' |
147 |
|
148 |
c open first file and get dimensions (float number and time) |
149 |
c |
150 |
imax=(6+4*Nr) |
151 |
ilen=imax*8 |
152 |
|
153 |
allocate (tmp(imax)) |
154 |
c |
155 |
write(dataFname(1:80),'(2a,a)') |
156 |
& fName(1:IL),'.001.001.data' |
157 |
open(1,file=dataFname,status='old',form='unformatted' |
158 |
& ,access='direct',recl=ilen) |
159 |
c |
160 |
read(1,rec=1) tmp |
161 |
rcountstart = SNGL(tmp(2)) |
162 |
rcountdelta = SNGL(tmp(4)) |
163 |
icount = INT(tmp(5)) |
164 |
npart = INT(tmp(6)) |
165 |
close(1) |
166 |
print*, 'npart = ',npart |
167 |
print*, 'timesteps= ',icount |
168 |
if (preflag) then |
169 |
print*, 'big data set --> read in chunks' |
170 |
endif |
171 |
|
172 |
|
173 |
c----------------------------------------------------------------------- |
174 |
c allocate variables |
175 |
c----------------------------------------------------------------------- |
176 |
c |
177 |
allocate (pnum(npart)) |
178 |
allocate (time(icount)) |
179 |
allocate (Dep_t(Nr)) |
180 |
allocate (Dep_w(Nr+1)) |
181 |
allocate (Dep_wm1(Nr)) |
182 |
allocate (xpart(npart,icount)) |
183 |
allocate (ypart(npart,icount)) |
184 |
allocate (kpart(npart,icount)) |
185 |
allocate (temp(npart,Nr,icount)) |
186 |
if (.not. preflag) then |
187 |
allocate (uvel(npart,Nr,icount)) |
188 |
allocate (vvel(npart,Nr,icount)) |
189 |
allocate (salt(npart,Nr,icount)) |
190 |
endif |
191 |
allocate (zeta(npart,icount)) |
192 |
|
193 |
c initialize arrays |
194 |
c |
195 |
do m=1,npart |
196 |
do n=1,icount |
197 |
xpart(m,n) = spval |
198 |
ypart(m,n) = spval |
199 |
kpart(m,n) = spval |
200 |
do k=1,Nr |
201 |
if (.not. preflag) uvel(m,k,n) = spval |
202 |
if (.not. preflag) vvel(m,k,n) = spval |
203 |
temp(m,k,n) = spval |
204 |
if (.not. preflag) salt(m,k,n) = spval |
205 |
enddo |
206 |
zeta(m,n) = spval |
207 |
enddo |
208 |
enddo |
209 |
c |
210 |
c |
211 |
c test if depth axis is evenly spaced (in that case no edges have to |
212 |
c be set) |
213 |
c |
214 |
iedge=.false. |
215 |
do k=2,Nr |
216 |
if (delZ(k) .ne. delZ(k-1)) then |
217 |
iedge=.true. |
218 |
goto 20 |
219 |
endif |
220 |
enddo |
221 |
20 continue |
222 |
c |
223 |
Dep_w(1)=0. |
224 |
Dep_wm1(1)=0. |
225 |
do k=2,Nr+1 |
226 |
Dep_w(k)=Dep_w(k-1)+delZ(k-1) |
227 |
if (k .ne. Nr+1) Dep_wm1(k) = Dep_w(k) |
228 |
enddo |
229 |
c |
230 |
do k=1,Nr |
231 |
Dep_t(k)=(Dep_w(k)+Dep_w(k+1))*0.5 |
232 |
enddo |
233 |
c |
234 |
c generate axes |
235 |
c |
236 |
time(1)=rcountstart |
237 |
do m=2,icount |
238 |
time(m) = time(m-1)+rcountdelta |
239 |
enddo |
240 |
print*, 'time: ',time |
241 |
c |
242 |
do ip=1,npart |
243 |
pnum(ip) = FLOAT(ip) |
244 |
enddo |
245 |
c print*, 'pnum: ',pnum |
246 |
c |
247 |
c |
248 |
c----------------------------------------------------------------------- |
249 |
c open files and read input |
250 |
c----------------------------------------------------------------------- |
251 |
c |
252 |
itotalrecord = 0 |
253 |
|
254 |
do iG=1,iGmax |
255 |
do jG=1,jGmax |
256 |
c |
257 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
258 |
& fName(1:IL),'.',iG,'.',jG,'.data' |
259 |
open(1,file=dataFname,status='old',form='unformatted' |
260 |
& ,access='direct',recl=ilen) |
261 |
c |
262 |
read(1,rec=1) tmp |
263 |
imaxrecord = INT(tmp(1)) |
264 |
print*,'read ',dataFname,imaxrecord |
265 |
itotalrecord = itotalrecord + imaxrecord |
266 |
c goto 200 |
267 |
|
268 |
do irec=2,imaxrecord+1 |
269 |
|
270 |
read(1,rec=irec) tmp |
271 |
ip = INT(tmp(1)) |
272 |
if (ip .gt. npart) then |
273 |
print*,'ip out of order: ',ip,npart |
274 |
stop |
275 |
endif |
276 |
np = INT((tmp(2)-rcountstart)/rcountdelta+1) |
277 |
|
278 |
|
279 |
if (usingSphericalPolarGrid) then |
280 |
xpart(ip,np) = SNGL(tmp(3)) |
281 |
ypart(ip,np) = SNGL(tmp(4)) |
282 |
else |
283 |
xpart(ip,np) = SNGL(tmp(3))/1000. |
284 |
ypart(ip,np) = SNGL(tmp(4))/1000. |
285 |
endif |
286 |
kpart(ip,np) = SNGL(tmp(5)) |
287 |
zeta(ip,np) = SNGL(tmp(6)) |
288 |
do k=1,Nr |
289 |
if (.not. preflag) uvel(ip,k,np) = SNGL(tmp(6+k)) |
290 |
if (.not. preflag) vvel(ip,k,np) = SNGL(tmp(6+1*Nr+k)) |
291 |
temp(ip,k,np) = SNGL(tmp(6+2*Nr+k)) |
292 |
if (.not. preflag) salt(ip,k,np) = SNGL(tmp(6+3*Nr+k)) |
293 |
if (temp(ip,k,np) .eq. 0.) then |
294 |
if (.not. preflag) uvel(ip,k,np) = spval |
295 |
if (.not. preflag) vvel(ip,k,np) = spval |
296 |
temp(ip,k,np) = spval |
297 |
if (.not. preflag) salt(ip,k,np) = spval |
298 |
endif |
299 |
enddo |
300 |
c print*, 'rec= ',irec,' npart= ',ip,' timestep= ',np |
301 |
c & ,time(np),tmp |
302 |
c & ,xpart(ip,np),ypart(ip,np),kpart(ip,np),uvel(ip,np,1), |
303 |
c & vvel(ip,np,1),temp(ip,np,1),salt(ip,np,1),zeta(ip,np) |
304 |
100 continue |
305 |
enddo |
306 |
|
307 |
close(1) |
308 |
200 continue |
309 |
enddo |
310 |
enddo |
311 |
|
312 |
print*,icount,' x ',npart,' = ',icount*npart,' records expected,', |
313 |
& ' found ',itotalrecord,' float records' |
314 |
print*,'==> ',icount*npart-itotalrecord,' float records missing' |
315 |
c |
316 |
c----------------------------------------------------------------------- |
317 |
c define netCDF variables |
318 |
c----------------------------------------------------------------------- |
319 |
c |
320 |
write(stdout,*) ' Start Converting' |
321 |
c |
322 |
c enter define mode: NCCLOB=overwrite, NCNOCLOB=do not overwrite |
323 |
c |
324 |
IL=ILNBLNK( fname ) |
325 |
outname2=fname(1:IL)//'.cdf' |
326 |
write (stdout,*) |
327 |
& ' ==> Writing a profiles to file ',outname2(1:IL+4) |
328 |
ncid = nccre(outname2(1:IL+4), NCCLOB, iret) |
329 |
if (iret .ne. 0) write(stdout,*) 'Error: Open NetCDF file' |
330 |
c |
331 |
c define dimensions |
332 |
c |
333 |
partdim = ncddef(ncid, 'Particles', npart, iret) |
334 |
Dep_tdim = ncddef(ncid, 'Depth_t', Nr, iret) |
335 |
Dep_wm1dim= ncddef(ncid, 'Depth_wm1', Nr, iret) |
336 |
Dep_wdim = ncddef(ncid, 'Depth_w', Nr+1, iret) |
337 |
Timedim = ncddef(ncid, 'Time', NCUNLIM, iret) |
338 |
if (iret .ne. 0) write(stdout,*) 'Error: define dimensions' |
339 |
c |
340 |
c define variables |
341 |
c |
342 |
dims(1) = partdim |
343 |
partid = ncvdef (ncid,'Particles',NCFLOAT,1,dims,iret) |
344 |
dims(1) = Dep_tdim |
345 |
Dep_tid = ncvdef (ncid,'Depth_t', NCFLOAT,1,dims,iret) |
346 |
dims(1) = Dep_wdim |
347 |
Dep_wid = ncvdef (ncid,'Depth_w', NCFLOAT,1,dims,iret) |
348 |
dims(1) = Dep_wm1dim |
349 |
Dep_wm1id= ncvdef (ncid,'Depth_wm1', NCFLOAT,1,dims,iret) |
350 |
dims(1) = Timedim |
351 |
Timeid = ncvdef (ncid,'Time', NCFLOAT,1,dims,iret) |
352 |
if (iret .ne. 0) write(stdout,*) 'Error: define axis ids' |
353 |
c |
354 |
dims(1) = partdim |
355 |
dims(2) = Timedim |
356 |
xpartid = ncvdef (ncid,'xpart', NCFLOAT,2,dims,iret) |
357 |
ypartid = ncvdef (ncid,'ypart', NCFLOAT,2,dims,iret) |
358 |
kpartid = ncvdef (ncid,'kpart', NCFLOAT,2,dims,iret) |
359 |
zetaid = ncvdef (ncid,'zeta', NCFLOAT,2,dims,iret) |
360 |
c |
361 |
dims(1) = partdim |
362 |
dims(2) = Dep_tdim |
363 |
dims(3) = Timedim |
364 |
uvelid = ncvdef (ncid,'uvel', NCFLOAT,3,dims,iret) |
365 |
vvelid = ncvdef (ncid,'vvel', NCFLOAT,3,dims,iret) |
366 |
tempid = ncvdef (ncid,'temp', NCFLOAT,3,dims,iret) |
367 |
saltid = ncvdef (ncid,'salt', NCFLOAT,3,dims,iret) |
368 |
if (iret .ne. 0) write(stdout,*) 'Error: define variable ids' |
369 |
c |
370 |
c----------------------------------------------------------------------- |
371 |
c assign attributes |
372 |
c----------------------------------------------------------------------- |
373 |
c |
374 |
charval(1) = ' ' |
375 |
c |
376 |
name = 'Particles Number ' |
377 |
c unit = 'particle number ' |
378 |
call ncaptc(ncid, partid, 'long_name', NCCHAR, 24, name, iret) |
379 |
call ncaptc(ncid, partid, 'units', NCCHAR, 16, unit, iret) |
380 |
c |
381 |
name = 'Time' |
382 |
unit = 'seconds' |
383 |
call ncaptc(ncid, Timeid, 'long_name', NCCHAR, 24, name, iret) |
384 |
call ncaptc(ncid, Timeid, 'units', NCCHAR, 16, unit, iret) |
385 |
call ncaptc(ncid, Timeid,'time_origin',NCCHAR, 20,startDate, iret) |
386 |
if (iret .ne. 0) write(stdout,*) 'Error: assign axis attributes' |
387 |
c |
388 |
c |
389 |
name = 'Depth of T grid points ' |
390 |
unit = 'meters ' |
391 |
call ncaptc(ncid, Dep_tid, 'long_name', NCCHAR, 24, name, iret) |
392 |
call ncaptc(ncid, Dep_tid, 'units', NCCHAR, 16, unit, iret) |
393 |
call ncaptc(ncid, Dep_tid, 'positive', NCCHAR, 4, 'down',iret) |
394 |
c call ncaptc(ncid, Dep_tid, 'point_spacing',NCCHAR,6,'uneven',iret) |
395 |
if (iedge) |
396 |
& call ncaptc(ncid, Dep_tid, 'edges',NCCHAR, 7,'Depth_w',iret) |
397 |
c |
398 |
name = 'Depth at top of T box' |
399 |
unit = 'meters ' |
400 |
call ncaptc(ncid, Dep_wm1id, 'long_name', NCCHAR, 24, name, iret) |
401 |
call ncaptc(ncid, Dep_wm1id, 'units', NCCHAR, 16, unit, iret) |
402 |
call ncaptc(ncid, Dep_wm1id, 'positive', NCCHAR, 4, 'down',iret) |
403 |
c |
404 |
name = 'Depth at bottom of T box' |
405 |
unit = 'meters ' |
406 |
call ncaptc(ncid, Dep_wid, 'long_name', NCCHAR, 24, name, iret) |
407 |
call ncaptc(ncid, Dep_wid, 'units', NCCHAR, 16, unit, iret) |
408 |
call ncaptc(ncid, Dep_wid, 'positive', NCCHAR, 4, 'down',iret) |
409 |
floatval(1) = spval |
410 |
c |
411 |
if (usingSphericalPolarGrid) then |
412 |
name = 'LONGITUDE ' |
413 |
unit = 'degrees_W ' |
414 |
else |
415 |
name = 'X_t ' |
416 |
unit = 'kilometer ' |
417 |
endif |
418 |
call ncaptc(ncid, xpartid, 'long_name', NCCHAR, 24, name, iret) |
419 |
call ncaptc(ncid, xpartid, 'units', NCCHAR, 16, unit, iret) |
420 |
call ncapt (ncid, xpartid,'missing_value',NCFLOAT,1,floatval,iret) |
421 |
call ncapt (ncid, xpartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
422 |
c |
423 |
if (usingSphericalPolarGrid) then |
424 |
name = 'LATITUDE ' |
425 |
unit = 'degrees_N ' |
426 |
else |
427 |
name = 'Y_t ' |
428 |
unit = 'kilometer ' |
429 |
endif |
430 |
call ncaptc(ncid, ypartid, 'long_name', NCCHAR, 24, name, iret) |
431 |
call ncaptc(ncid, ypartid, 'units', NCCHAR, 16, unit, iret) |
432 |
call ncapt (ncid, ypartid,'missing_value',NCFLOAT,1,floatval,iret) |
433 |
call ncapt (ncid, ypartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
434 |
c |
435 |
name = 'LEVEL ' |
436 |
unit = 'LEVEL ' |
437 |
call ncaptc(ncid, kpartid, 'long_name', NCCHAR, 24, name, iret) |
438 |
call ncaptc(ncid, kpartid, 'units', NCCHAR, 16, unit, iret) |
439 |
call ncapt (ncid, kpartid,'missing_value',NCFLOAT,1,floatval,iret) |
440 |
call ncapt (ncid, kpartid,'_FillValue', NCFLOAT, 1,floatval, iret) |
441 |
c |
442 |
name = 'POTENTIAL TEMPERATURE ' |
443 |
unit = 'deg C ' |
444 |
call ncaptc(ncid, tempid, 'long_name', NCCHAR, 24, name, iret) |
445 |
call ncaptc(ncid, tempid, 'units', NCCHAR, 16, unit, iret) |
446 |
call ncapt (ncid, tempid, 'missing_value',NCFLOAT,1,floatval,iret) |
447 |
call ncapt (ncid, tempid, '_FillValue', NCFLOAT, 1,floatval, iret) |
448 |
c |
449 |
name = 'SALINITY ' |
450 |
unit = 'PSU ' |
451 |
call ncaptc(ncid, saltid, 'long_name', NCCHAR, 24, name, iret) |
452 |
call ncaptc(ncid, saltid, 'units', NCCHAR, 16, unit, iret) |
453 |
call ncapt (ncid, saltid, 'missing_value',NCFLOAT,1,floatval,iret) |
454 |
call ncapt (ncid, saltid, '_FillValue', NCFLOAT, 1,floatval, iret) |
455 |
c |
456 |
name = 'U VELOCITY ' |
457 |
unit = 'm/sec' |
458 |
call ncaptc(ncid, uvelid, 'long_name', NCCHAR, 24, name, iret) |
459 |
call ncaptc(ncid, uvelid, 'units', NCCHAR, 16, unit, iret) |
460 |
call ncapt (ncid, uvelid, 'missing_value',NCFLOAT,1,floatval,iret) |
461 |
call ncapt (ncid, uvelid, '_FillValue', NCFLOAT, 1,floatval, iret) |
462 |
c |
463 |
name = 'V VELOCITY ' |
464 |
unit = 'm/sec' |
465 |
call ncaptc(ncid, vvelid, 'long_name', NCCHAR, 24, name, iret) |
466 |
call ncaptc(ncid, vvelid, 'units', NCCHAR, 16, unit, iret) |
467 |
call ncapt (ncid, vvelid, 'missing_value',NCFLOAT,1,floatval,iret) |
468 |
call ncapt (ncid, vvelid, '_FillValue', NCFLOAT, 1,floatval, iret) |
469 |
c |
470 |
name = 'SEA SURFACE HEIGHT ' |
471 |
unit = 'm ' |
472 |
call ncaptc(ncid, zetaid, 'long_name', NCCHAR, 24, name, iret) |
473 |
call ncaptc(ncid, zetaid, 'units', NCCHAR, 16, unit, iret) |
474 |
call ncapt (ncid, zetaid,'missing_value',NCFLOAT, 1,floatval,iret) |
475 |
call ncapt (ncid, zetaid,'_FillValue', NCFLOAT, 1, floatval, iret) |
476 |
c |
477 |
if (iret .ne. 0) write(stdout,*) 'Error: define variable attrib.' |
478 |
c |
479 |
expname= ' ' |
480 |
stamp = ' ' |
481 |
call ncaptc(ncid, NCGLOBAL, 'title', NCCHAR, 60, expnam, iret) |
482 |
call ncaptc(ncid, NCGLOBAL, 'history', NCCHAR, 32, stamp, iret) |
483 |
c |
484 |
c----------------------------------------------------------------------- |
485 |
c leave define mode |
486 |
c----------------------------------------------------------------------- |
487 |
c |
488 |
call ncendf(ncid, iret) |
489 |
c |
490 |
c |
491 |
c----------------------------------------------------------------------- |
492 |
c put variables in netCDF file |
493 |
c----------------------------------------------------------------------- |
494 |
c |
495 |
c store Particles |
496 |
corner(1) = 1 |
497 |
edges(1) = npart |
498 |
call ncvpt(ncid, partid, corner, edges, pnum, iret) |
499 |
c |
500 |
c store Time |
501 |
corner(1) = 1 |
502 |
edges(1) = icount |
503 |
call ncvpt(ncid, Timeid, corner, edges, Time, iret) |
504 |
c |
505 |
c store Depth_t |
506 |
corner(1) = 1 |
507 |
edges(1) = Nr |
508 |
call ncvpt(ncid, Dep_tid, corner, edges, Dep_t, iret) |
509 |
c store Depth_w |
510 |
corner(1) = 1 |
511 |
edges(1) = Nr+1 |
512 |
call ncvpt(ncid, Dep_wid, corner, edges, Dep_w, iret) |
513 |
c store Depth_wm1 |
514 |
corner(1) = 1 |
515 |
edges(1) = Nr |
516 |
call ncvpt(ncid, Dep_wm1id, corner, edges, Dep_wm1, iret) |
517 |
c store 2D values |
518 |
corner(1) = 1 |
519 |
corner(2) = 1 |
520 |
edges(1) = npart |
521 |
edges(2) = icount |
522 |
VARid=xpartid |
523 |
call ncvpt(ncid, VARid, corner, edges, xpart, iret) |
524 |
VARid=ypartid |
525 |
call ncvpt(ncid, VARid, corner, edges, ypart, iret) |
526 |
VARid=kpartid |
527 |
call ncvpt(ncid, VARid, corner, edges, kpart, iret) |
528 |
VARid=zetaid |
529 |
call ncvpt(ncid, VARid, corner, edges, zeta, iret) |
530 |
c store values |
531 |
corner(1) = 1 |
532 |
corner(2) = 1 |
533 |
corner(3) = 1 |
534 |
edges(1) = npart |
535 |
edges(2) = Nr |
536 |
edges(3) = icount |
537 |
VARid=tempid |
538 |
call ncvpt(ncid, VARid, corner, edges, temp, iret) |
539 |
|
540 |
if (preflag) then |
541 |
c read in salt into temp array |
542 |
do iG=1,iGmax |
543 |
do jG=1,jGmax |
544 |
c |
545 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
546 |
& fName(1:IL),'.',iG,'.',jG,'.data' |
547 |
open(1,file=dataFname,status='old',form='unformatted' |
548 |
& ,access='direct',recl=ilen) |
549 |
c |
550 |
read(1,rec=1) tmp |
551 |
imaxrecord = INT(tmp(1)) |
552 |
print*,'read salt from ',dataFname |
553 |
do irec=2,imaxrecord+1 |
554 |
|
555 |
read(1,rec=irec) tmp |
556 |
ip = INT(tmp(1)) |
557 |
np = INT((tmp(2)-rcountstart)/rcountdelta+1) |
558 |
do k=1,Nr |
559 |
temp(ip,k,np) = SNGL(tmp(6+3*Nr+k)) |
560 |
if (SNGL(tmp(6+2*Nr+k)) .eq. 0.) then |
561 |
temp(ip,k,np) = spval |
562 |
endif |
563 |
enddo |
564 |
enddo |
565 |
close(1) |
566 |
enddo |
567 |
enddo |
568 |
VARid=saltid |
569 |
call ncvpt(ncid, VARid, corner, edges, temp, iret) |
570 |
|
571 |
c read in u into temp array |
572 |
do iG=1,iGmax |
573 |
do jG=1,jGmax |
574 |
c |
575 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
576 |
& fName(1:IL),'.',iG,'.',jG,'.data' |
577 |
open(1,file=dataFname,status='old',form='unformatted' |
578 |
& ,access='direct',recl=ilen) |
579 |
c |
580 |
read(1,rec=1) tmp |
581 |
imaxrecord = INT(tmp(1)) |
582 |
print*,'read uvel from ',dataFname |
583 |
do irec=2,imaxrecord+1 |
584 |
|
585 |
read(1,rec=irec) tmp |
586 |
ip = INT(tmp(1)) |
587 |
np = INT((tmp(2)-rcountstart)/rcountdelta+1) |
588 |
do k=1,Nr |
589 |
temp(ip,k,np) = SNGL(tmp(6+k)) |
590 |
if (SNGL(tmp(6+2*Nr+k)) .eq. 0.) then |
591 |
temp(ip,k,np) = spval |
592 |
endif |
593 |
enddo |
594 |
enddo |
595 |
close(1) |
596 |
enddo |
597 |
enddo |
598 |
VARid=uvelid |
599 |
call ncvpt(ncid, VARid, corner, edges, temp, iret) |
600 |
|
601 |
c read in v into temp array |
602 |
do iG=1,iGmax |
603 |
do jG=1,jGmax |
604 |
c |
605 |
write(dataFname(1:80),'(2a,i3.3,a,i3.3,a)') |
606 |
& fName(1:IL),'.',iG,'.',jG,'.data' |
607 |
open(1,file=dataFname,status='old',form='unformatted' |
608 |
& ,access='direct',recl=ilen) |
609 |
c |
610 |
read(1,rec=1) tmp |
611 |
imaxrecord = INT(tmp(1)) |
612 |
print*,'read vvel from ',dataFname |
613 |
do irec=2,imaxrecord+1 |
614 |
|
615 |
read(1,rec=irec) tmp |
616 |
ip = INT(tmp(1)) |
617 |
np = INT((tmp(2)-rcountstart)/rcountdelta+1) |
618 |
do k=1,Nr |
619 |
temp(ip,k,np) = SNGL(tmp(6+1*Nr+k)) |
620 |
if (SNGL(tmp(6+2*Nr+k)) .eq. 0.) then |
621 |
temp(ip,k,np) = spval |
622 |
endif |
623 |
enddo |
624 |
enddo |
625 |
close(1) |
626 |
enddo |
627 |
enddo |
628 |
VARid=vvelid |
629 |
call ncvpt(ncid, VARid, corner, edges, temp, iret) |
630 |
|
631 |
else |
632 |
|
633 |
VARid=saltid |
634 |
call ncvpt(ncid, VARid, corner, edges, salt, iret) |
635 |
VARid=uvelid |
636 |
call ncvpt(ncid, VARid, corner, edges, uvel, iret) |
637 |
VARid=vvelid |
638 |
call ncvpt(ncid, VARid, corner, edges, vvel, iret) |
639 |
|
640 |
endif |
641 |
c |
642 |
if (iret .ne. 0) write(stdout,*) 'Error: write variables' |
643 |
call ncclos (ncid, iret) |
644 |
c |
645 |
write(stdout,*) ' End ' |
646 |
|
647 |
end |
648 |
|
649 |
|
650 |
INTEGER FUNCTION ILNBLNK( string ) |
651 |
C /==========================================================\ |
652 |
C | FUNCTION ILNBLNK | |
653 |
C | o Find last non-blank in character string. | |
654 |
C \==========================================================/ |
655 |
IMPLICIT NONE |
656 |
CHARACTER*(*) string |
657 |
CEndOfInterface |
658 |
INTEGER L, LS |
659 |
C |
660 |
LS = LEN(string) |
661 |
ILNBLNK = LS |
662 |
DO 10 L = LS, 1, -1 |
663 |
IF ( string(L:L) .EQ. ' ' ) GOTO 10 |
664 |
ILNBLNK = L |
665 |
GOTO 11 |
666 |
10 CONTINUE |
667 |
11 CONTINUE |
668 |
C |
669 |
RETURN |
670 |
END |