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

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

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


Revision 1.1 - (hide annotations) (download)
Thu Apr 2 17:49:59 2009 UTC (15 years, 2 months ago) by dfer
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint62, checkpoint63, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint61n, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61q, checkpoint61z, checkpoint61x, checkpoint61y, HEAD
Christopher Wolfe's new version of cvfloat.F,
update Makefile and data.float accordingly

1 dfer 1.1 program cvfloat
2     !
3     !=======================================================================
4     !
5     ! converts binary float trajectories to netCDF
6     !
7     !=======================================================
8     !
9     ! * uses namelist data.float
10     !
11     ! Arne Biastoch, abiastoch@ucsd.edu, 11/16/2000
12     ! Updated 03/16/2009 Christopher Wolfe, clwolfe@ucsd.edu
13     !
14     !=======================================================================
15     !
16     use netcdf
17     implicit none
18     integer, parameter :: STDIN = 5, STDOUT = 6, STDERR = 6
19    
20     character(80) :: iotext
21     character(32) :: stamp
22     integer :: ioun,iargc,m,n,ilen,irec,ip,np
23     integer :: icount,itotalrecord,imaxrecord
24    
25     ! tile counters
26     integer :: iGmax,jGmax,iG,jG
27    
28     ! variables for filenames
29     integer :: narg, npart
30     character(80) :: dataFName
31     logical :: exst
32    
33     ! number of variables per record
34     integer, parameter :: IMAX = 13
35     !
36     ! integer narg
37     logical preflag
38     !
39     ! netCDF ids
40     !
41     integer :: ncid
42     integer :: partdim,Timedim
43     integer :: partid, Timeid
44     integer :: xpartid, ypartid, zpartid
45     integer :: ipartid, jpartid, kpartid
46     integer tempid, saltid, uvelid, vvelid, presid
47     !
48     ! attribute vectors
49     !
50     character(24) :: name*24
51     character(16) :: unit*16
52     !
53     ! data variables for NetCDF
54     !
55     double precision :: rcountstart, rcountdelta
56     double precision, dimension(:), allocatable :: pnum,time
57     double precision, dimension(:,:), allocatable :: xpart,ypart,zpart, &
58     ipart,jpart,kpart, &
59     pres,uvel,vvel,temp,salt
60     double precision, dimension(:), allocatable :: tmp
61     !
62     ! namelist contains
63     !
64     character(50) :: outname2, fName="float_trajectories", outname
65     character(20) :: startDate="01-JAN-1992:12:00:00"
66     character(60) :: expnam = "Experiment name not set"
67     logical :: usingSphericalPolarGrid=.true.
68    
69     namelist /dimensions/ expnam, startDate, usingSphericalPolarGrid
70     namelist /floats/ fName
71    
72     ioun=11
73     open(ioun,file="data.float",status="old",form="formatted")
74     read (unit=ioun, end=666, nml=dimensions)
75     write (STDOUT,dimensions)
76     close (ioun)
77     666 continue
78     open(ioun,file="data.float",status="old",form="formatted")
79     read (unit=ioun, end=999, nml=floats)
80     write (STDOUT,floats)
81     close (ioun)
82     999 continue
83    
84     !
85     ! PRELIMINARY USE:
86     ! IF FLOATS SHOULD BE VIEWED DURING A CURRENT MODEL RUN THE FIRST
87     ! LINE OF THE FILE MAY NOT BE UPDATED CORRECTLY, I.E. THERE MIGHT
88     ! BE MORE TIMES THAN STATED AT THE BEGINNING. BY GIVING A FLAG
89     ! ONLY ICOUNT-1 TIMESTEPS ARE USED
90     !
91     preflag = .false.
92     narg=iargc()
93     if ( narg > 0 ) preflag = .true.
94     !
95     ! check existent files
96     !
97     iGmax=1
98     do m=1,100
99     write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",iGmax,".",1,".data"
100     inquire( file=dataFname, exist=exst )
101     if (exst) iGmax = iGmax + 1
102     enddo
103    
104     jGmax=1
105     do m=1,100
106     write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",1,".",jGmax,".data"
107     inquire( file=dataFname, exist=exst )
108     if (exst) jGmax = jGmax + 1
109     enddo
110    
111     iGmax = iGmax - 1
112     jGmax = jGmax - 1
113     print*, "There are ",iGmax," x ",jGmax," files"
114    
115     ! open first file and get dimensions (float number and time)
116     ilen=IMAX*8
117     allocate (tmp(IMAX))
118    
119     write(dataFname(1:80),"(2a,a)") trim(fName),".001.001.data"
120     open(1,file=dataFname,status="old",form="unformatted",access="direct",recl=ilen)
121    
122     read(1,rec=1) tmp
123     ! print*,"tmp:", tmp
124    
125     rcountstart = tmp(2)
126     rcountdelta = tmp(4)
127     icount = INT(tmp(5))
128     npart = INT(tmp(6))
129     close(1)
130    
131     print*, "npart = ",npart
132     print*, "timesteps= ",icount
133     ! print*,"rcountstart=",rcountstart,"rcountdelta",rcountdelta
134     if (preflag) then
135     icount=icount-1
136     print*, "preliminary --> use timesteps= ",icount
137     endif
138    
139     !-----------------------------------------------------------------------
140     ! allocate variables
141     !-----------------------------------------------------------------------
142     !
143     allocate (pnum(npart))
144     allocate (time(icount))
145     allocate (xpart(npart,icount))
146     allocate (ypart(npart,icount))
147     allocate (zpart(npart,icount))
148     allocate (ipart(npart,icount))
149     allocate (jpart(npart,icount))
150     allocate (kpart(npart,icount))
151     allocate (uvel(npart,icount))
152     allocate (vvel(npart,icount))
153     allocate (temp(npart,icount))
154     allocate (salt(npart,icount))
155     allocate (pres(npart,icount))
156    
157     ! initialize arrays
158     !
159     do m=1,npart
160     do n=1,icount
161     xpart(m,n) = NF90_FILL_DOUBLE
162     ypart(m,n) = NF90_FILL_DOUBLE
163     zpart(m,n) = NF90_FILL_DOUBLE
164     ipart(m,n) = NF90_FILL_DOUBLE
165     jpart(m,n) = NF90_FILL_DOUBLE
166     kpart(m,n) = NF90_FILL_DOUBLE
167     uvel(m,n) = NF90_FILL_DOUBLE
168     vvel(m,n) = NF90_FILL_DOUBLE
169     temp(m,n) = NF90_FILL_DOUBLE
170     salt(m,n) = NF90_FILL_DOUBLE
171     pres(m,n) = NF90_FILL_DOUBLE
172     enddo
173     enddo
174    
175     ! generate axes
176     !
177     time(1)=rcountstart
178     do m=2,icount
179     time(m) = time(m-1)+rcountdelta
180     enddo
181     print*, "time: ",time(1:icount)
182    
183     do ip=1,npart
184     pnum(ip) = FLOAT(ip)
185     enddo
186     !
187     !-----------------------------------------------------------------------
188     ! open files and read input
189     !-----------------------------------------------------------------------
190     !
191    
192     itotalrecord = 0
193    
194     do iG=1,iGmax
195     do jG=1,jGmax
196    
197     write(dataFname(1:80),"(2a,i3.3,a,i3.3,a)") trim(fName),".",iG,".",jG,".data"
198     open(1,file=dataFname,status="old",form="unformatted",access="direct",recl=ilen)
199    
200     read(1,rec=1) tmp
201     imaxrecord = INT(tmp(1))
202     print "(1x,2a)","read ",dataFname
203     itotalrecord = itotalrecord + imaxrecord
204    
205     do irec=2,imaxrecord+1
206    
207     read(1,rec=irec) tmp
208     ip = INT(tmp(1))
209     if (ip > npart) then
210     print*,"ip out of order: ",ip,npart
211     stop
212     endif
213    
214     ! Note: If the float initial conditions are also written to the data files the time interval
215     ! between record 1 and 2 may be somewhat less than rcountdelta. The +0.9999 deals with
216     ! this case without upsetting the rest of the indexing.
217     np = INT((tmp(2)-rcountstart)/rcountdelta + 1 + 0.9999)
218    
219     ! this is only for prelimiray results. Use only icount-1 timesteps
220     if (preflag .and. (np > icount .or. np < 1)) cycle
221    
222     if (usingSphericalPolarGrid) then
223     xpart(ip,np) = tmp(3)
224     ypart(ip,np) = tmp(4)
225     else
226     xpart(ip,np) = tmp(3)
227     ypart(ip,np) = tmp(4)
228     endif
229     zpart(ip,np) = tmp(5)
230     ipart(ip,np) = tmp(6)
231     jpart(ip,np) = tmp(7)
232     kpart(ip,np) = tmp(8)
233     pres(ip,np) = tmp(9)
234     uvel(ip,np) = tmp(10)
235     vvel(ip,np) = tmp(11)
236     temp(ip,np) = tmp(12)
237     salt(ip,np) = tmp(13)
238     enddo
239    
240     close(1)
241     enddo
242     enddo
243    
244     print*,icount," x ",npart," = ",icount*npart," records expected,",&
245     " found ",itotalrecord," float records"
246     print*,"==> ",icount*npart-itotalrecord," float records missing"
247     !
248     !
249     !-----------------------------------------------------------------------
250     ! define netCDF variables
251     !-----------------------------------------------------------------------
252     !
253     write(STDOUT,*) " Start Converting"
254     !
255     ! enter define mode
256     !
257     outname2=trim(fname)//".cdf"
258     write (STDOUT,*)" ==> Writing a trajectories to file ",trim(outname2)
259     call nc_check(nf90_create(trim(outname2),NF90_CLOBBER,ncid),__LINE__)
260     !
261     ! define dimensions
262     !
263     call nc_check(nf90_def_dim(ncid, "Particles", npart, partdim),__LINE__)
264     call nc_check(nf90_def_dim(ncid, "Time", NF90_UNLIMITED, Timedim),__LINE__)
265     !
266     ! define variables
267     !
268     call nc_check(nf90_def_var(ncid, "Particles",NF90_DOUBLE, (/partdim/), partid),__LINE__)
269     call nc_check(nf90_def_var(ncid, "Time", NF90_DOUBLE, (/Timedim/), Timeid),__LINE__)
270    
271     call nc_check(nf90_def_var(ncid, "x", NF90_DOUBLE, (/partdim,Timedim/), xpartid),__LINE__)
272     call nc_check(nf90_def_var(ncid, "y", NF90_DOUBLE, (/partdim,Timedim/), ypartid),__LINE__)
273     call nc_check(nf90_def_var(ncid, "z", NF90_DOUBLE, (/partdim,Timedim/), zpartid),__LINE__)
274     call nc_check(nf90_def_var(ncid, "i", NF90_DOUBLE, (/partdim,Timedim/), ipartid),__LINE__)
275     call nc_check(nf90_def_var(ncid, "j", NF90_DOUBLE, (/partdim,Timedim/), jpartid),__LINE__)
276     call nc_check(nf90_def_var(ncid, "k", NF90_DOUBLE, (/partdim,Timedim/), kpartid),__LINE__)
277     call nc_check(nf90_def_var(ncid, "pressure", NF90_DOUBLE, (/partdim,Timedim/), presid),__LINE__)
278     call nc_check(nf90_def_var(ncid, "u", NF90_DOUBLE, (/partdim,Timedim/), uvelid),__LINE__)
279     call nc_check(nf90_def_var(ncid, "v", NF90_DOUBLE, (/partdim,Timedim/), vvelid),__LINE__)
280     call nc_check(nf90_def_var(ncid, "T", NF90_DOUBLE, (/partdim,Timedim/), tempid),__LINE__)
281     call nc_check(nf90_def_var(ncid, "S", NF90_DOUBLE, (/partdim,Timedim/), saltid),__LINE__)
282     !
283     !-----------------------------------------------------------------------
284     ! assign attributes
285     !-----------------------------------------------------------------------
286     !
287     name = "Particle Number"
288     unit = "particle number"
289     call nc_check(nf90_put_att(ncid,partid,"long_name",name),__LINE__)
290     call nc_check(nf90_put_att(ncid,partid,"units", unit),__LINE__)
291    
292     name = "Time"
293     unit = "seconds"
294     call nc_check(nf90_put_att(ncid,Timeid,"long_name", name),__LINE__)
295     call nc_check(nf90_put_att(ncid,Timeid,"units", unit),__LINE__)
296     call nc_check(nf90_put_att(ncid,Timeid,"time_origin",startDate),__LINE__)
297    
298     if (usingSphericalPolarGrid) then
299     name = "LONGITUDE "
300     unit = "degrees_W "
301     else
302     name = "X_t "
303     unit = "meter "
304     endif
305     call nc_check(nf90_put_att(ncid,xpartid,"long_name", name),__LINE__)
306     call nc_check(nf90_put_att(ncid,xpartid,"units", unit),__LINE__)
307     call nc_check(nf90_put_att(ncid,xpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
308     call nc_check(nf90_put_att(ncid,xpartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
309    
310     if (usingSphericalPolarGrid) then
311     name = "LATITUDE "
312     unit = "degrees_N "
313     else
314     name = "Y_t "
315     unit = "meter "
316     endif
317     call nc_check(nf90_put_att(ncid,ypartid,"long_name", name),__LINE__)
318     call nc_check(nf90_put_att(ncid,ypartid,"units", unit),__LINE__)
319     call nc_check(nf90_put_att(ncid,ypartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
320     call nc_check(nf90_put_att(ncid,ypartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
321    
322     name = "DEPTH "
323     unit = "meter "
324     call nc_check(nf90_put_att(ncid,zpartid,"long_name", name),__LINE__)
325     call nc_check(nf90_put_att(ncid,zpartid,"units", unit),__LINE__)
326     call nc_check(nf90_put_att(ncid,zpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
327     call nc_check(nf90_put_att(ncid,zpartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
328    
329     name = "iINDEX "
330     unit = "index "
331     call nc_check(nf90_put_att(ncid,ipartid,"long_name", name),__LINE__)
332     call nc_check(nf90_put_att(ncid,ipartid,"units", unit),__LINE__)
333     call nc_check(nf90_put_att(ncid,ipartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
334     call nc_check(nf90_put_att(ncid,ipartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
335    
336     name = "jINDEX "
337     unit = "index "
338     call nc_check(nf90_put_att(ncid,jpartid,"long_name", name),__LINE__)
339     call nc_check(nf90_put_att(ncid,jpartid,"units", unit),__LINE__)
340     call nc_check(nf90_put_att(ncid,jpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
341     call nc_check(nf90_put_att(ncid,jpartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
342    
343     name = "LEVEL "
344     unit = "level "
345     call nc_check(nf90_put_att(ncid,kpartid,"long_name", name),__LINE__)
346     call nc_check(nf90_put_att(ncid,kpartid,"units", unit),__LINE__)
347     call nc_check(nf90_put_att(ncid,kpartid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
348     call nc_check(nf90_put_att(ncid,kpartid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
349    
350     name = "POTENTIAL TEMPERATURE "
351     unit = "deg C "
352     call nc_check(nf90_put_att(ncid,tempid,"long_name", name),__LINE__)
353     call nc_check(nf90_put_att(ncid,tempid,"units", unit),__LINE__)
354     call nc_check(nf90_put_att(ncid,tempid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
355     call nc_check(nf90_put_att(ncid,tempid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
356    
357     name = "SALINITY "
358     unit = "PSU "
359     call nc_check(nf90_put_att(ncid,saltid,"long_name", name),__LINE__)
360     call nc_check(nf90_put_att(ncid,saltid,"units", unit),__LINE__)
361     call nc_check(nf90_put_att(ncid,saltid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
362     call nc_check(nf90_put_att(ncid,saltid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
363    
364     name = "U VELOCITY "
365     unit = "m/sec"
366     call nc_check(nf90_put_att(ncid,uvelid,"long_name", name),__LINE__)
367     call nc_check(nf90_put_att(ncid,uvelid,"units", unit),__LINE__)
368     call nc_check(nf90_put_att(ncid,uvelid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
369     call nc_check(nf90_put_att(ncid,uvelid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
370    
371     name = "V VELOCITY "
372     unit = "m/sec"
373     call nc_check(nf90_put_att(ncid,vvelid,"long_name", name),__LINE__)
374     call nc_check(nf90_put_att(ncid,vvelid,"units", unit),__LINE__)
375     call nc_check(nf90_put_att(ncid,vvelid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
376     call nc_check(nf90_put_att(ncid,vvelid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
377    
378     name = "PRESSURE "
379     unit = "N/m^2 "
380     call nc_check(nf90_put_att(ncid,presid,"long_name", name),__LINE__)
381     call nc_check(nf90_put_att(ncid,presid,"units", unit),__LINE__)
382     call nc_check(nf90_put_att(ncid,presid,"missing_value",NF90_FILL_DOUBLE),__LINE__)
383     call nc_check(nf90_put_att(ncid,presid,"_FillValue", NF90_FILL_DOUBLE),__LINE__)
384    
385    
386     ! expnam= " "
387     ! stamp = " "
388     call nc_check(nf90_put_att(ncid,NF90_GLOBAL,"title", trim(expnam)),__LINE__)
389     ! call nc_check(nf90_put_att(ncid,NF90_GLOBAL,"history",stamp),__LINE__)
390     !
391     !-----------------------------------------------------------------------
392     ! leave define mode
393     !-----------------------------------------------------------------------
394     !
395     call nc_check(nf90_enddef(ncid),__LINE__)
396     !
397     !
398     !-----------------------------------------------------------------------
399     ! put variables in netCDF file
400     !-----------------------------------------------------------------------
401     !
402     ! store Particles
403     call nc_check(nf90_put_var(ncid,partid,pnum),__LINE__)
404     !
405     ! store Time
406     call nc_check(nf90_put_var(ncid,Timeid,Time),__LINE__)
407     !
408     ! store values
409     call nc_check(nf90_put_var(ncid,xpartid,xpart),__LINE__)
410     call nc_check(nf90_put_var(ncid,ypartid,ypart),__LINE__)
411     call nc_check(nf90_put_var(ncid,zpartid,zpart),__LINE__)
412     call nc_check(nf90_put_var(ncid,ipartid,ipart),__LINE__)
413     call nc_check(nf90_put_var(ncid,jpartid,jpart),__LINE__)
414     call nc_check(nf90_put_var(ncid,kpartid,kpart),__LINE__)
415     call nc_check(nf90_put_var(ncid,tempid, temp),__LINE__)
416     call nc_check(nf90_put_var(ncid,saltid, salt),__LINE__)
417     call nc_check(nf90_put_var(ncid,uvelid, uvel),__LINE__)
418     call nc_check(nf90_put_var(ncid,vvelid, vvel),__LINE__)
419     call nc_check(nf90_put_var(ncid,presid, pres),__LINE__)
420    
421     call nc_check(nf90_close(ncid),__LINE__)
422    
423     write(STDOUT,*) " End "
424    
425     !-----------------------------------------------------------------------
426     ! internal subroutines
427     !-----------------------------------------------------------------------
428     contains
429     subroutine nc_check(status,lineno)
430     integer, intent ( in) :: status
431     integer, intent ( in) :: lineno
432    
433     if(status /= nf90_noerr) then
434     print *, "Error at line number ",lineno,":"
435     print *, trim(nf90_strerror(status))
436     stop 2
437     end if
438     end subroutine nc_check
439     end program cvfloat

  ViewVC Help
Powered by ViewVC 1.1.22