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

Last change on this file since 451 was 386, checked in by bellier, 16 years ago

Added CeCILL License information

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