/[MITgcm]/MITgcm/verification/flt_example/extra/cvfloat.F
ViewVC logotype

Annotation of /MITgcm/verification/flt_example/extra/cvfloat.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Thu May 6 02:01:34 2004 UTC (20 years, 2 months ago) by edhill
Branch: MAIN
CVS Tags: checkpoint57t_post, checkpoint58l_post, checkpoint53f_post, checkpoint54a_pre, checkpoint55c_post, checkpoint53b_pre, checkpoint57o_post, checkpoint57m_post, checkpoint57i_post, checkpoint58e_post, checkpoint57r_post, checkpoint52n_post, checkpoint57f_post, checkpoint57s_post, checkpoint57j_post, checkpoint58b_post, checkpoint58m_post, checkpoint57b_post, checkpoint53c_post, checkpoint53d_post, checkpoint57f_pre, checkpoint57k_post, checkpoint55d_pre, checkpoint57d_post, checkpoint57g_post, checkpoint60, checkpoint61, checkpoint57c_pre, checkpoint58r_post, checkpoint55j_post, checkpoint56b_post, checkpoint57h_pre, checkpoint57y_post, checkpoint58g_post, checkpoint57x_post, checkpoint54a_post, checkpoint55h_post, checkpoint58n_post, checkpoint58x_post, checkpoint57g_pre, checkpoint54b_post, checkpoint58h_post, checkpoint57e_post, checkpoint58w_post, checkpoint54d_post, checkpoint56c_post, checkpoint54e_post, checkpoint58j_post, checkpoint55b_post, checkpoint57h_post, checkpoint57y_pre, checkpoint55, checkpoint53a_post, checkpoint55a_post, checkpoint57a_post, checkpoint57v_post, checkpoint59q, checkpoint59p, checkpoint55g_post, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint55f_post, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint59, checkpoint58, checkpoint57a_pre, checkpoint54, checkpoint57, checkpoint56, checkpoint53, checkpoint57h_done, checkpoint58f_post, checkpoint53g_post, checkpoint58d_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58o_post, checkpoint57z_post, checkpoint54f_post, checkpoint55e_post, eckpoint57e_pre, checkpoint58c_post, checkpoint58k_post, checkpoint57c_post, checkpoint58u_post, checkpoint58y_post, checkpoint53b_post, checkpoint58v_post, checkpoint53d_pre, checkpoint58s_post, checkpoint61f, checkpoint61g, checkpoint61d, checkpoint61e, checkpoint61b, checkpoint61c, checkpoint58p_post, checkpoint61a, checkpoint58t_post, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint58q_post, checkpoint55i_post, checkpoint54c_post, checkpoint57l_post, checkpoint56a_post, checkpoint55d_post
 o copying the contents of "aux" to "extra" so that MITgcm is easier
   to use within pathetic "OSes" that cannot deal with file or
   directory names such as "aux"

1 edhill 1.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

  ViewVC Help
Powered by ViewVC 1.1.22