source: BOL/IPCC_AR4/ts2IPCC.F90 @ 5008

Last change on this file since 5008 was 605, checked in by Laurent Fairhead, 20 years ago

Un probleme sur la boucle de lecture du fichier de correspondance
variable IPSL - variable IPCC empechait l'execution sur PC
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.3 KB
Line 
1!
2! $Header$
3!
4PROGRAM ts2IPCC
5
6!
7! Filtre permettant de transformer les fichiers de serie temporelle a une
8! variable de l'IPSL en fichiers acceptables par le PCMDI/IPCC.
9!
10! Utilisation de la bibliothèque CMOR du PCMDI
11!
12! L. Fairhead 2004/08
13!
14! Ce programme est appelé avec un argument, le nom du fichier à traiter
15! Il nécessite aussi un fichier config.def contenant diverses informations
16! (voir plus bas) et l'accès à un tableau faisant la correspondance entre
17! les noms de variables du modèle et les noms imposés par l'IPCC
18! Pour l'instant on ne traite que les fichiers contenant la serie
19! temporelle d'une seule variable
20
21  use cmor_users_functions
22  use netcdf
23
24  implicit none
25
26#include "netcdf.inc"
27
28  CHARACTER (len=256)        :: orig_file   ! nom du fichier à traiter
29  character (len=512)        :: line_read
30  CHARACTER (len=128)        :: inpath, contact, hist_gen,repert,instit
31  CHARACTER (len=128)        :: hist_var,expt_id,source,comment,refs
32  CHARACTER (len=20)         :: action
33  character (len=20)         :: first_part
34  character (len=1004)       :: second_part
35  CHARACTER (len=20), DIMENSION(100) :: ipsl_name, ipsl_units, ipsl_pos
36  CHARACTER (len=20), DIMENSION(100) :: ipcc_name, ipcc_table
37  CHARACTER (len=80)         :: varname, units, namedim
38
39  INTEGER                    :: orig_file_id, nvars, ndims
40  INTEGER                    :: verbos, exit_ctl, realis, indice,index_table
41  INTEGER                    :: iargc, iostat, ierr
42  INTEGER                    :: i,idim
43  INTEGER, ALLOCATABLE, DIMENSION(:)       :: dimids,axis_ids,lendim
44  INTEGER                    :: latid, lonid, vertid, timeid
45  INTEGER                    :: varid, cmorvarid
46  INTEGER                    :: ilat, ilon, ivert, itime
47  INTEGER                    :: lunout      ! device de sortie
48
49  logical                    :: found = .false.
50
51  REAL, ALLOCATABLE, DIMENSION(:) :: lon, lat, vert, time
52  REAL, ALLOCATABLE, DIMENSION(:) :: lon_bounds, lat_bounds
53  REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: donnees
54  DOUBLE PRECISION, DIMENSION(1)                :: rdate
55  real                            :: missing_value
56
57  external iargc
58
59!
60! quelques initialisations
61  lunout = 6
62  varname = 'xxxxxxxx'
63
64!
65! On vérifie que l'appel au programme a bien un argument:
66  CALL getarg(1, orig_file)
67  IF (iargc() == 0 .OR.  orig_file == '-h') then
68    WRITE(lunout,*)' '
69    WRITE(lunout,*)' Utilisation de ce programme: '
70    WRITE(lunout,*)' ./ts2IPCC nom_de_fichier [variable]'
71    WRITE(lunout,*)'        ou nom_de_fichier est le nom du fichier a traiter'
72    WRITE(lunout,*)'        et variable la variable a traiter [optionel]'
73    WRITE(lunout,*)' '
74    WRITE(lunout,*)' ./ts2IPCC -h sort ce message'
75    WRITE(lunout,*)' '
76    stop
77  ENDIF
78  if (iargc() == 2) then
79    CALL getarg(2, varname)
80  endif
81
82!
83! Lecture du fichier de configuration
84  OPEN (20, IOSTAT=iostat, file='config.def',form='formatted')
85  IF (iostat /= 0) then
86    WRITE(lunout,*)'Erreur ouverture du fichier config.def'
87    stop
88  endif
89
90  config: do
91    READ(20,'(A)',iostat=iostat)line_read
92    IF (iostat /= 0) exit
93    line_read = TRIM(line_read)
94    IF (INDEX(line_read, '#') /= 1) THEN
95      first_part = trim(line_read(1:INDEX(line_read, '=')-1))
96      second_part = trim(line_read(INDEX(line_read, '=')+1:))
97      selectcase(first_part)
98        case('inpath')
99         inpath = trim(second_part)
100        case('file_action')
101         action = trim(second_part)
102        case('verbosity')
103          READ(second_part,'(i)') verbos
104        case('exit_control')
105          READ(second_part,'(i)') exit_ctl
106        case('repertoire')
107         repert = trim(second_part)
108        case('experiment_ID')
109         expt_id = trim(second_part)
110        case('institut')
111         instit = trim(second_part)
112        case('source')
113         source = trim(second_part)
114        case('realisation')
115          READ(second_part,'(i)') realis
116        case('hist_gen')
117         hist_gen = trim(second_part)
118        case('comment')
119         comment = trim(second_part)
120        case('refs')
121         refs = trim(second_part)
122        case('hist_var')
123         hist_var = trim(second_part)
124        case('contact')
125          contact = trim(second_part)
126      end select
127    endif
128  enddo config
129  if (iostat > 0) then
130    WRITE(lunout,*)'Probleme de lecture du fichier config.def, iostat = ',iostat
131    stop
132  endif
133  close(20)
134
135!
136! Lecture du tableau de correspondance nom IPSL <=> nom IPCC
137  OPEN (20, IOSTAT=iostat, file='table.def',form='formatted')
138  IF (iostat /= 0) then
139    WRITE(lunout,*)'Erreur ouverture du fichier table.def'
140    stop
141  endif
142  indice = 0
143  table: do
144    READ(20,'(A)',iostat=iostat)line_read
145    IF (iostat /= 0) exit
146    line_read = TRIM(line_read)
147    IF (INDEX(line_read, '#') /= 1) THEN
148      indice = indice + 1
149      ipsl_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
150      line_read = trim(line_read(INDEX(line_read, '|')+1:))
151      ipsl_units(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
152      line_read = trim(line_read(INDEX(line_read, '|')+1:))
153      ipsl_pos(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
154      line_read = trim(line_read(INDEX(line_read, '|')+1:))
155      ipcc_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
156      ipcc_table(indice) = trim(line_read(INDEX(line_read, '|')+1:))
157    endif
158  enddo table
159  indice = indice - 1
160  close(20)
161!  DO i = 1, indice
162!    WRITE(lunout,*)ipsl_name(i),ipsl_units(i),ipcc_name(i),ipcc_table(i)
163!  enddo
164
165!
166! Ouverture du fichier a traiter
167  ierr = nf_open(orig_file, NF_NOWRITE, orig_file_id)
168  IF (ierr /= NF_NOERR) then
169    WRITE(lunout,*)NF_STRERROR(ierr)
170    stop
171  endif
172!
173! trouver la variable a traiter, c'est une variable a 3 ou 4 dimensions
174  ierr = nf_inq_nvars(orig_file_id, nvars)
175  IF (ierr /= NF_NOERR) then
176    WRITE(lunout,*)NF_STRERROR(ierr)
177    stop
178  endif
179
180  i = 0
181  if (varname == 'xxxxxxxx') then
182    DO while (.not.found)
183      i = i + 1
184      if (i > nvars) then
185        WRITE(lunout,*)' pas de variable 3d ou 4d trouvee'
186        stop
187      endif
188      ierr = nf_inq_varname(orig_file_id, i, varname)
189      IF (ierr /= NF_NOERR) then
190        WRITE(lunout,*)NF_STRERROR(ierr)
191        stop
192      endif
193      ierr =  nf_inq_varndims(orig_file_id, i, ndims)
194      IF (ierr /= NF_NOERR) then
195        WRITE(lunout,*)NF_STRERROR(ierr)
196        stop
197      endif
198      if (ndims > 2) found = .true.
199    enddo
200  else
201    ierr = nf_inq_varid(orig_file_id, varname, varid)
202      IF (ierr /= NF_NOERR) then
203        WRITE(lunout,*)NF_STRERROR(ierr)
204        stop
205      endif
206  endif
207   
208!
209! recherche de la correspondance nom IPSL <=> nom IPCC
210  found = .false.
211  i = 0
212  do while (.not. found)
213    i = i + 1
214    if (i > indice) then
215      WRITE(lunout,*)'La variable ',trim(varname),' n''est pas dans le tableau de correspondance table.def'
216      stop
217    endif 
218    IF (varname == ipsl_name(i)) THEN
219      index_table = i
220      found = .true.
221    endif
222  enddo
223
224  WRITE(lunout,*)' found variable = ', trim(varname)
225  WRITE(lunout,*)' ipcc_name = ', trim(ipcc_name(index_table))
226
227
228!
229! Initialisation CMOR
230  ierr = cmor_setup(inpath=inpath, netcdf_file_action=action,set_verbosity=verbos,&
231 &                  exit_control=exit_ctl)
232  IF (ierr /= 0) then
233    WRITE(lunout,*)'Probleme dans cmor_setup, ierr = ', ierr
234  endif
235
236!
237! Initialisation dataset
238  ierr = cmor_dataset(outpath=repert,        &
239 &                    experiment_id=expt_id, &
240 &                    institution=instit,    &
241 &                    source=source,         &
242 &                    calendar='360_day',    &
243 &                    realization=realis,    &
244 &                    contact=contact,       &
245 &                    history=hist_gen,      &
246 &                    comment=comment,       &
247 &                    references=refs)
248  IF (ierr /= 0) then
249    WRITE(lunout,*)'Probleme dans cmor_dataset, ierr = ', ierr
250  endif
251
252!
253! Definition des axes
254
255  ierr = nf90_inq_varid(orig_file_id,TRIM(varname), varid)
256  if (ierr /= 0) call handle_err(ierr)
257  ierr = nf90_Inquire_Variable(orig_file_id, varid, ndims = ndims)
258  if (ierr /= 0) call handle_err(ierr)
259  allocate (dimids(ndims))
260  allocate (axis_ids(ndims))
261  allocate (lendim(ndims))
262  ierr = nf90_Inquire_Variable(orig_file_id, varid, dimids = dimids)
263  if (ierr /= 0) call handle_err(ierr)
264
265  do idim = 1, ndims
266    ierr = nf90_Inquire_Dimension(orig_file_id, dimids(idim), &
267                          name = namedim, len = lendim(idim))
268    if (ierr /= 0) call handle_err(ierr)
269    units=' '
270    selectcase(trim(namedim))
271      case('lat')
272!     lecture de la latitude:
273        allocate(lat(lendim(idim)))
274        ierr = nf_inq_varid(orig_file_id, namedim, latid)
275        if (ierr /= 0) call handle_err(ierr)
276        ierr = nf_get_var_real(orig_file_id, latid, lat)
277        if (ierr /= 0) call handle_err(ierr)
278        ierr = nf_get_att_text(orig_file_id, latid, 'units', units)
279        if (ierr /= 0) call handle_err(ierr)
280        allocate(lat_bounds(lendim(idim)+1))
281        DO i = 2, lendim(idim)
282          lat_bounds(i) = lat(i-1) - (lat(i-1) - lat(i))/2
283        enddo
284        lat_bounds(1) = lat(1)
285        lat_bounds(lendim(idim)+1) = lat(lendim(idim))
286!       definition de la latitude
287        axis_ids(idim) = cmor_axis(                           &
288                    table=trim(ipcc_table(index_table)),      &
289                    table_entry='latitude',                   &
290                    units=trim(units),                        & 
291                    length=lendim(idim),                      &
292                    coord_vals=lat,                           &
293                    cell_bounds=lat_bounds)
294!
295!       
296      case('lon')
297!       lecture de la longitude:
298        allocate(lon(lendim(idim)))
299        ierr = nf_inq_varid(orig_file_id, namedim, lonid)
300        if (ierr /= 0) call handle_err(ierr)
301        ierr = nf_get_var_real(orig_file_id, lonid, lon)
302        if (ierr /= 0) call handle_err(ierr)
303        ierr = nf_get_att_text(orig_file_id, lonid, 'units', units)
304        if (ierr /= 0) call handle_err(ierr)
305        ALLOCATE(lon_bounds(lendim(idim)+1))
306        DO i = 2, lendim(idim)
307         lon_bounds(i) = lon(i-1) - (lon(i-1) - lon(i))/2
308        enddo
309        lon_bounds(1) = lon(1) - (lon_bounds(3) -lon_bounds(2))/2.
310        lon_bounds(lendim(idim)+1) = lon(lendim(idim)) + (lon_bounds(lendim(idim))-lon_bounds(lendim(idim)-1))/2.
311
312!       definition de la longitude
313        axis_ids(idim) = cmor_axis(                           &
314                    table=trim(ipcc_table(index_table)),      &
315                    table_entry='longitude',                  &
316                    units=trim(units),                        & 
317                    length=lendim(idim),                            &
318                    coord_vals=lon,                           &
319                    cell_bounds=lon_bounds)       
320!
321!
322      case('presnivs')
323!     lecture de la verticale:
324        allocate(vert(lendim(idim)))
325        ierr = nf_inq_varid(orig_file_id, namedim, vertid)
326        if (ierr /= 0) call handle_err(ierr)
327        ierr = nf_get_var_real(orig_file_id, vertid, vert)
328        if (ierr /= 0) call handle_err(ierr)
329        ierr = nf_get_att_text(orig_file_id, vertid, 'units', units)
330        if (ierr /= 0) call handle_err(ierr)
331!
332!       definition de la verticale
333        if (units == 'mb') units='Pa'
334        axis_ids(idim) = cmor_axis(                           &
335                    table=trim(ipcc_table(index_table)),      &
336                    table_entry='pressure',                   &
337                    units=trim(units),                        & 
338                    length=lendim(idim),                      &
339                    coord_vals=vert)
340!
341!       
342      case('time_counter')         
343!     definition du temps
344      if (idim /= ndims) then
345        write(lunout,*)'la dimension temps doit etre la derniere dimension'
346        stop
347      endif
348      allocate(time(lendim(idim)))
349      ierr = nf_inq_varid(orig_file_id,namedim,timeid)
350      if (ierr /=0) call handle_err(ierr)
351      ierr = nf_get_var_real(orig_file_id, timeid, time)
352      if (ierr /=0) call handle_err(ierr)
353      ierr = nf_get_att_text(orig_file_id,timeid, 'units', units)
354      if (ierr /=0) call handle_err(ierr)
355      axis_ids(idim) = cmor_axis(                          &
356           table=trim(ipcc_table(index_table)),            &
357           table_entry='time',                             &
358           units=trim(units),                              &
359           length=lendim(idim),                            &
360           interval='30 minutes')
361      itime= lendim(idim)
362    case default
363      write(lunout,*)'Dimension: ', trim(namedim),' non reconnue'
364      stop
365   endselect
366  enddo
367 
368!
369! Definition de la variable a ecrire
370  units=' '
371  ierr = nf_inq_varid(orig_file_id,TRIM(varname), varid)
372  if (ierr /= 0) call handle_err(ierr)
373  ierr = nf_get_att_text(orig_file_id, varid, 'units', units)
374  if (ierr /= 0) call handle_err(ierr)
375  ierr = nf_get_att_real(orig_file_id, varid, 'missing_value', missing_value)
376  if (ierr /= 0) call handle_err(ierr)
377  cmorvarid = cmor_variable(                         &
378       table=trim(ipcc_table(index_table)),          &
379       table_entry=trim(ipcc_name(index_table)),     &
380       units=trim(units),                            &
381       axis_ids= axis_ids,                           &
382       missing_value=real(missing_value),            &
383       positive = trim(ipsl_pos(index_table)),       &
384       original_name=trim(varname))
385!
386! Lecture de la variable
387  if (ndims == 3) then
388    ALLOCATE (donnees(lendim(1), lendim(2), 1, lendim(3) ))
389  else if (ndims ==4) then
390    allocate (donnees(lendim(1), lendim(2), lendim(3), lendim(4) ))
391  endif
392  ierr = nf_get_var_real(orig_file_id, varid, donnees)
393!
394! Ecriture de la variable
395 
396  DO i = 1, itime
397    rdate(1) = dble(time(i))
398    ierr = cmor_write(                                     &
399             var_id        = cmorvarid,                    &
400             DATA          = REAL(donnees(:,:,:,i)),       &
401             ntimes_passed = 1,                            &
402             time_vals     = rdate)
403  enddo
404!
405! Fin CMOR
406  ierr = cmor_close()
407  IF (ierr /= 0) then
408    WRITE(lunout,*)'Probleme dans cmor_close, ierr = ', ierr
409  endif
410
411!
412! fermeture fichier originel
413  ierr = nf_close(orig_file_id)
414  IF (ierr /= NF_NOERR) then
415    WRITE(lunout,*)NF_STRERROR(ierr)
416    stop
417  endif
418
419contains
420
421subroutine handle_err(status)
422  use netcdf
423
424  implicit none
425  integer, intent(in) :: status
426
427  write(lunout,*)nf90_strerror(status)
428  stop
429
430end subroutine handle_err
431 
432END PROGRAM ts2IPCC
433
Note: See TracBrowser for help on using the repository browser.