source: IOIPSL/trunk/tools/ncregular.f90 @ 15

Last change on this file since 15 was 15, checked in by bellier, 15 years ago

JB: add tools and examples
dd

File size: 9.4 KB
Line 
1PROGRAM ncregular
2!
3!$Id$
4!---------------------------------------------------------------------
5!- This code replaces a 2D surface grid by vectors.
6!- Obviously it only works if you have a regular grid.
7!-
8!- Jan Polcher (polcher@lmd.jussieu.fr)
9!- Jacques Bellier (jacques.bellier@cea.fr)
10!---------------------------------------------------------------------
11  USE netcdf
12!-
13  IMPLICIT NONE
14!-
15  INTEGER :: iread, if, in, iv, sz
16  INTEGER :: ier, nb_files, iret, ndims, nvars, nb_glat
17  INTEGER :: lon_dim_id, lat_dim_id
18  INTEGER :: lon_len, lat_len, lon_id, lat_id
19  INTEGER :: nav_lon_id, nav_lat_id
20  INTEGER :: alloc_stat_lon, alloc_stat_lat
21!-
22  INTEGER,ALLOCATABLE ::  file_id(:), tax_id(:)
23  CHARACTER(LEN=80),ALLOCATABLE :: names(:)
24  CHARACTER(LEN=80) :: dim_name
25  CHARACTER(LEN=80) :: varname
26  CHARACTER(LEN=20) :: xname, yname, lonname, latname
27  LOGICAL :: check, regular
28!-
29  REAL,ALLOCATABLE :: lon(:), lat(:), lon2(:), lat2(:)
30  REAL,ALLOCATABLE :: del_lon(:), del_lat(:)
31!-
32  INTEGER iargc, getarg
33  EXTERNAL iargc, getarg
34!---------------------------------------------------------------------
35  alloc_stat_lon = 0
36  alloc_stat_lat = 0
37!-
38  iread = iargc()
39!-
40  ALLOCATE (names(iread),stat=ier)
41  IF (ier /= 0) THEN
42    WRITE (*,*) ' Could not allocate names of size ', iread
43    STOP 'nctax'
44  ENDIF
45!-
46  CALL nct_getarg (iread, nb_files, names, check, &
47 &                 xname, yname, lonname, latname)
48!-
49! Allocate space
50!-
51  ALLOCATE (file_id(nb_files),stat=ier)
52  IF (ier /= 0) THEN
53    WRITE (*,*) ' Could not allocate file_id of size ', nb_files
54    STOP 'nctax'
55  ENDIF
56!-
57  ALLOCATE (tax_id(nb_files),stat=ier)
58  IF (ier /= 0) THEN
59    WRITE (*,*) ' Could not allocate tax_id of size ', nb_files
60    STOP 'nctax'
61  ENDIF
62!-
63  DO if=1,nb_files
64!---
65    IF (check) THEN
66      WRITE(*,*) 'ncregular : ', if,  names(if)
67    ENDIF
68!---
69    iret = NF90_OPEN (names(if),NF90_WRITE,file_id(if))
70    iret = NF90_INQUIRE (file_id(if),ndims,nvars,nb_glat,tax_id(if))
71!---
72!-- Get the IDs of the variables
73!---
74    lon_len = -9999
75    lat_len = -9999
76    DO in=1,ndims
77!-----
78      iret = NF90_INQUIRE_DIMENSION (file_id(if), in, dim_name, sz)
79!-----
80      IF (     (LEN_TRIM(dim_name) == 1) &
81 &        .AND.(INDEX(dim_name,TRIM(xname)) == 1) ) THEN
82        lon_dim_id = in
83        lon_len = sz
84      ENDIF
85!-----
86      IF (     (LEN_TRIM(dim_name) == 1) &
87 &        .AND.(INDEX(dim_name,TRIM(yname)) == 1) ) THEN
88        lat_dim_id = in
89        lat_len = sz
90      ENDIF
91!-----
92    ENDDO
93!---
94    IF ( (lon_len == -9999).OR.(lat_len == -9999) ) THEN
95      WRITE(*,*) 'ncregular : The specified dimensions were not'
96      WRITE(*,*) 'found in file : ',names(if)
97      iret = NF90_CLOSE (file_id(if))
98      STOP
99    ENDIF
100!---
101    IF (check) THEN
102      WRITE(*,*) 'ncregular : lon_dim_id, lon_len',lon_dim_id,lon_len
103      WRITE(*,*) 'ncregular : lat_dim_id, lat_len',lat_dim_id,lat_len
104    ENDIF
105!---
106!-- Look for the right variables
107!---
108    nav_lon_id = -9999
109    nav_lat_id = -9999
110    DO iv=1,nvars
111      iret = NF90_INQUIRE_VARIABLE (file_id(if),iv,name=varname)
112      IF (INDEX(varname,TRIM(lonname)) > 0) THEN
113        nav_lon_id = iv
114      ENDIF
115      IF (INDEX(varname,TRIM(latname)) > 0) THEN
116        nav_lat_id = iv
117      ENDIF
118    ENDDO
119!---
120    IF ( (nav_lon_id == -9999).OR.(nav_lat_id == -9999) ) THEN
121      WRITE(*,*) 'ncregular : The specified coordinate fields'
122      WRITE(*,*) 'were not found in file : ',names(if)
123      iret = NF90_CLOSE (file_id(if))
124      STOP
125    ENDIF
126!---
127    IF (check) THEN
128      WRITE(*,*) 'ncregular : nav_lon_id :', nav_lon_id
129      WRITE(*,*) 'ncregular : nav_lat_id :', nav_lat_id
130    ENDIF
131!---
132!-- Read variables from file and check if regular
133!---
134!-- Do we have the variable to read the
135!---
136    IF ( alloc_stat_lon < lon_len) THEN
137      IF ( alloc_stat_lon > 0) THEN
138        deallocate(lon)
139        deallocate(lon2)
140        deallocate(del_lon)
141      ENDIF
142      allocate(lon(lon_len))
143      allocate(lon2(lon_len))
144      allocate(del_lon(lon_len))
145      alloc_stat_lon = lon_len
146    ENDIF
147!---
148    IF ( alloc_stat_lat < lat_len) THEN
149      IF ( alloc_stat_lat > 0) THEN
150        deallocate(lat)
151        deallocate(lat2)
152        deallocate(del_lat)
153      ENDIF
154      allocate(lat(lat_len))
155      allocate(lat2(lat_len))
156      allocate(del_lat(lat_len))
157      alloc_stat_lat = lat_len
158    ENDIF
159!---
160!-- Read data
161!---
162    iret = NF90_GET_VAR (file_id(if),nav_lon_id,lon, &
163 & start=(/1,1/),count=(/lon_len,1/),stride=(/1,1/))
164    iret = NF90_GET_VAR (file_id(if),nav_lon_id,lon2, &
165 & start=(/1,int(lat_len/2)/),count=(/lon_len,1/),stride=(/1,1/))
166    del_lon = lon-lon2
167!-
168    iret = NF90_GET_VAR (file_id(if),nav_lat_id,lat, &
169 & start=(/1,1/),count=(/1,lat_len/),stride=(/lon_len,1/))
170    iret = NF90_GET_VAR (file_id(if),nav_lat_id,lat2, &
171 & start=(/int(lon_len/2),1/),count=(/1,lat_len/),stride=(/lon_len,1/))
172    del_lat = lat-lat2
173!-
174    regular = (    (MAXVAL(del_lon) < 0.001) &
175 &             .OR.(MAXVAL(del_lat) < 0.001) )
176!---
177!-- Create the new variables
178!---
179    IF (regular) THEN
180      IF (check) THEN
181        WRITE(*,*) 'Regular case'
182      ENDIF
183      iret = NF90_REDEF (file_id(if))
184      iret = NF90_RENAME_DIM (file_id(if), lon_dim_id, 'lon')
185      iret = NF90_RENAME_DIM (file_id(if), lat_dim_id, 'lat')
186      IF (check) THEN
187        WRITE(*,*) 'Dimensions renamed'
188      ENDIF
189      iret = NF90_DEF_VAR (file_id(if), 'lon', NF90_FLOAT, &
190 &                         lon_dim_id, lon_id)
191      iret = NF90_DEF_VAR (file_id(if), 'lat', NF90_FLOAT, &
192 &                         lat_dim_id, lat_id)
193      IF (check) THEN
194        WRITE(*,*) 'New variables defined'
195      ENDIF
196!-----
197!---- Copy attributes
198!-----
199      iret = NF90_COPY_ATT (file_id(if),nav_lon_id,'units', &
200 &                          file_id(if),lon_id)
201      iret = NF90_COPY_ATT (file_id(if),nav_lon_id,'title', &
202 &                          file_id(if),lon_id)
203      iret = NF90_COPY_ATT (file_id(if),nav_lon_id,'valid_max', &
204 &                          file_id(if),lon_id)
205      iret = NF90_COPY_ATT (file_id(if),nav_lon_id,'valid_min', &
206 &                          file_id(if),lon_id)
207!-----
208      iret = NF90_COPY_ATT (file_id(if),nav_lat_id,'units', &
209 &                          file_id(if),lat_id)
210      iret = NF90_COPY_ATT (file_id(if),nav_lat_id,'title', &
211 &                          file_id(if),lat_id)
212      iret = NF90_COPY_ATT (file_id(if),nav_lat_id,'valid_max', &
213 &                          file_id(if),lat_id)
214      iret = NF90_COPY_ATT (file_id(if),nav_lat_id,'valid_min', &
215 &                          file_id(if),lat_id)
216!-----
217!---- Go into write mode
218!-----
219      iret = NF90_ENDDEF (file_id(if))
220!-----
221!---- Write data
222!-----
223      iret = NF90_PUT_VAR (file_id(if),lon_id,lon(1:lon_len))
224      iret = NF90_PUT_VAR (file_id(if),lat_id,lat(1:lat_len))
225!-
226      iret = NF90_CLOSE (file_id(if))
227    ELSE
228      WRITE(*,*) 'ncregular : Your grid is not regular'
229      WRITE(*,*) names(if), 'remains unchanged'
230      iret = NF90_CLOSE (file_id(if))
231    ENDIF
232!-
233  ENDDO
234!--------------------
235END PROGRAM ncregular
236!-
237!===
238!-
239SUBROUTINE nct_getarg (argx, nb_files, names, check, &
240 &                     xname, yname, lonname, latname)
241!---------------------------------------------------------------------
242!- Read the arguments of nctax.
243!---------------------------------------------------------------------
244  INTEGER,INTENT(in) :: argx
245  INTEGER, INTENT(out) :: nb_files
246  CHARACTER(LEN=80),INTENT(out) :: names(argx)
247  CHARACTER(LEN=20) :: xname, yname, lonname, latname
248!-
249  CHARACTER(LEN=80) :: tmp, tmp_arg
250  LOGICAL :: check
251!---------------------------------------------------------------------
252  check = .FALSE.
253!-
254! Get the number of arguments
255!-
256  nb_files = 0
257!-
258  xname = 'x'
259  yname = 'y'
260  lonname = 'nav_lon'
261  latname = 'nav_lat'
262!-
263! Go through the arguments and analyse them one by one
264!-
265  IF (check) WRITE(*,*) 'Start going through the arguments'
266!-
267  IF (argx == 0) THEN
268    WRITE(*,*) 'To get usage : nctax -h '
269    STOP
270  ENDIF
271!-
272  iread = 1
273  DO WHILE (iread <= argx)
274    iret = getarg(iread,tmp)
275    IF (check) WRITE(*,*) ' iread, tmp :', iread, tmp
276    SELECTCASE(tmp)
277      CASE('-d')
278        WRITE(*,*) 'DEBUG MODE SELECTED'
279        check = .TRUE.
280        iread = iread+1
281      CASE('-h')
282        WRITE(*,*) 'Usage : nregular [options] file1 [file2 ...]'
283        WRITE(*,*) '    -d : Verbose mode'
284        WRITE(*,*) '    -h : This output'
285        STOP
286      CASE('-dim_lon')
287        iread = iread+1
288        iret  = getarg(iread,tmp_arg)
289        xname = TRIM(tmp_arg)
290        iread = iread+1
291      CASE('-dim_lat')
292        iread = iread+1
293        iret  = getarg(iread,tmp_arg)
294        yname = TRIM(tmp_arg)
295        iread = iread+1
296      CASE('-coo_lon')
297        iread = iread+1
298        iret  = getarg(iread,tmp_arg)
299        lonname = TRIM(tmp_arg)
300        iread = iread+1
301      CASE('-coo_lat')
302        iread = iread+1
303        iret  = getarg(iread,tmp_arg)
304        latname = TRIM(tmp_arg)
305        iread = iread+1
306      CASE DEFAULT
307        IF (check) WRITE(*,*) 'nct_getarg : CASE default'
308        IF (INDEX(tmp,'-') /= 1) THEN
309          nb_files = nb_files+1
310          names(nb_files) = tmp
311          iread = iread+1
312        ELSE
313          WRITE(*,*) "WARNING Unknown option ",tmp
314          WRITE(*,*) "For ore information : nctax -h"
315        ENDIF
316    END SELECT
317  ENDDO
318!-
319  IF (check) THEN
320    WRITE(*,*) ' nct_getarg : output >> '
321    WRITE(*,*) '>> nb_files : ', nb_files
322    WRITE(*,*) '>> names :', (names(ii), ii=1,nb_files)
323  ENDIF
324!------------------------
325END SUBROUTINE nct_getarg
Note: See TracBrowser for help on using the repository browser.