source: BOL/trunk/IPCC/ts2IPCC.F90 @ 570

Last change on this file since 570 was 561, checked in by lmdzadmin, 20 years ago

Rajout d'une colonne indiquant dans quelle direction (up/down) les
variables sont positives (pour les variables pour lesquelles cela a un sens)
IM/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  do while (iostat == 0)
91    READ(20,'(A)',iostat=iostat)line_read
92    line_read = trim(line_read)
93    IF (INDEX(line_read, '#') /= 1) THEN
94      first_part = trim(line_read(1:INDEX(line_read, '=')-1))
95      second_part = trim(line_read(INDEX(line_read, '=')+1:))
96      selectcase(first_part)
97        case('inpath')
98         inpath = trim(second_part)
99        case('file_action')
100         action = trim(second_part)
101        case('verbosity')
102          READ(second_part,'(i)') verbos
103        case('exit_control')
104          READ(second_part,'(i)') exit_ctl
105        case('repertoire')
106         repert = trim(second_part)
107        case('experiment_ID')
108         expt_id = trim(second_part)
109        case('institut')
110         instit = trim(second_part)
111        case('source')
112         source = trim(second_part)
113        case('realisation')
114          READ(second_part,'(i)') realis
115        case('hist_gen')
116         hist_gen = trim(second_part)
117        case('comment')
118         comment = trim(second_part)
119        case('refs')
120         refs = trim(second_part)
121        case('hist_var')
122         hist_var = trim(second_part)
123        case('contact')
124          contact = trim(second_part)
125      end select
126    endif
127  enddo
128  if (iostat > 0) then
129    WRITE(lunout,*)'Probleme de lecture du fichier config.def, iostat = ',iostat
130    stop
131  endif
132  close(20)
133
134!
135! Lecture du tableau de correspondance nom IPSL <=> nom IPCC
136  OPEN (20, IOSTAT=iostat, file='table.def',form='formatted')
137  IF (iostat /= 0) then
138    WRITE(lunout,*)'Erreur ouverture du fichier table.def'
139    stop
140  endif
141  indice = 0
142  do while (iostat == 0)
143    READ(20,'(A)',iostat=iostat)line_read
144    line_read = trim(line_read)
145    IF (INDEX(line_read, '#') /= 1) THEN
146      indice = indice + 1
147      ipsl_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
148      line_read = trim(line_read(INDEX(line_read, '|')+1:))
149      ipsl_units(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
150      line_read = trim(line_read(INDEX(line_read, '|')+1:))
151      ipsl_pos(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
152      line_read = trim(line_read(INDEX(line_read, '|')+1:))
153      ipcc_name(indice) = trim(line_read(1:INDEX(line_read, '|')-1))
154      ipcc_table(indice) = trim(line_read(INDEX(line_read, '|')+1:))
155    endif
156  enddo
157  indice = indice - 1
158  close(20)
159!  DO i = 1, indice
160!    WRITE(lunout,*)ipsl_name(i),ipsl_units(i),ipcc_name(i),ipcc_table(i)
161!  enddo
162
163!
164! Ouverture du fichier a traiter
165  ierr = nf_open(orig_file, NF_NOWRITE, orig_file_id)
166  IF (ierr /= NF_NOERR) then
167    WRITE(lunout,*)NF_STRERROR(ierr)
168    stop
169  endif
170!
171! trouver la variable a traiter, c'est une variable a 3 ou 4 dimensions
172  ierr = nf_inq_nvars(orig_file_id, nvars)
173  IF (ierr /= NF_NOERR) then
174    WRITE(lunout,*)NF_STRERROR(ierr)
175    stop
176  endif
177
178  i = 0
179  if (varname == 'xxxxxxxx') then
180    DO while (.not.found)
181      i = i + 1
182      if (i > nvars) then
183        WRITE(lunout,*)' pas de variable 3d ou 4d trouvee'
184        stop
185      endif
186      ierr = nf_inq_varname(orig_file_id, i, varname)
187      IF (ierr /= NF_NOERR) then
188        WRITE(lunout,*)NF_STRERROR(ierr)
189        stop
190      endif
191      ierr =  nf_inq_varndims(orig_file_id, i, ndims)
192      IF (ierr /= NF_NOERR) then
193        WRITE(lunout,*)NF_STRERROR(ierr)
194        stop
195      endif
196      if (ndims > 2) found = .true.
197    enddo
198  else
199    ierr = nf_inq_varid(orig_file_id, varname, varid)
200      IF (ierr /= NF_NOERR) then
201        WRITE(lunout,*)NF_STRERROR(ierr)
202        stop
203      endif
204  endif
205   
206!
207! recherche de la correspondance nom IPSL <=> nom IPCC
208  found = .false.
209  i = 0
210  do while (.not. found)
211    i = i + 1
212    if (i > indice) then
213      WRITE(lunout,*)'La variable ',trim(varname),' n''est pas dans le tableau de correspondance table.def'
214      stop
215    endif 
216    IF (varname == ipsl_name(i)) THEN
217      index_table = i
218      found = .true.
219    endif
220  enddo
221
222  WRITE(lunout,*)' found variable = ', trim(varname)
223  WRITE(lunout,*)' ipcc_name = ', trim(ipcc_name(index_table))
224
225
226!
227! Initialisation CMOR
228  ierr = cmor_setup(inpath=inpath, netcdf_file_action=action,set_verbosity=verbos,&
229 &                  exit_control=exit_ctl)
230  IF (ierr /= 0) then
231    WRITE(lunout,*)'Probleme dans cmor_setup, ierr = ', ierr
232  endif
233
234!
235! Initialisation dataset
236  ierr = cmor_dataset(outpath=repert,        &
237 &                    experiment_id=expt_id, &
238 &                    institution=instit,    &
239 &                    source=source,         &
240 &                    calendar='360_day',    &
241 &                    realization=realis,    &
242 &                    contact=contact,       &
243 &                    history=hist_gen,      &
244 &                    comment=comment,       &
245 &                    references=refs)
246  IF (ierr /= 0) then
247    WRITE(lunout,*)'Probleme dans cmor_dataset, ierr = ', ierr
248  endif
249
250!
251! Definition des axes
252
253  ierr = nf90_inq_varid(orig_file_id,TRIM(varname), varid)
254  if (ierr /= 0) call handle_err(ierr)
255  ierr = nf90_Inquire_Variable(orig_file_id, varid, ndims = ndims)
256  if (ierr /= 0) call handle_err(ierr)
257  allocate (dimids(ndims))
258  allocate (axis_ids(ndims))
259  allocate (lendim(ndims))
260  ierr = nf90_Inquire_Variable(orig_file_id, varid, dimids = dimids)
261  if (ierr /= 0) call handle_err(ierr)
262
263  do idim = 1, ndims
264    ierr = nf90_Inquire_Dimension(orig_file_id, dimids(idim), &
265                          name = namedim, len = lendim(idim))
266    if (ierr /= 0) call handle_err(ierr)
267    units=' '
268    selectcase(trim(namedim))
269      case('lat')
270!     lecture de la latitude:
271        allocate(lat(lendim(idim)))
272        ierr = nf_inq_varid(orig_file_id, namedim, latid)
273        if (ierr /= 0) call handle_err(ierr)
274        ierr = nf_get_var_real(orig_file_id, latid, lat)
275        if (ierr /= 0) call handle_err(ierr)
276        ierr = nf_get_att_text(orig_file_id, latid, 'units', units)
277        if (ierr /= 0) call handle_err(ierr)
278        allocate(lat_bounds(lendim(idim)+1))
279        DO i = 2, lendim(idim)
280          lat_bounds(i) = lat(i-1) - (lat(i-1) - lat(i))/2
281        enddo
282        lat_bounds(1) = lat(1)
283        lat_bounds(lendim(idim)+1) = lat(lendim(idim))
284!       definition de la latitude
285        axis_ids(idim) = cmor_axis(                           &
286                    table=trim(ipcc_table(index_table)),      &
287                    table_entry='latitude',                   &
288                    units=units,                              & 
289                    length=lendim(idim),                      &
290                    coord_vals=lat,                           &
291                    cell_bounds=lat_bounds)
292!
293!       
294      case('lon')
295!       lecture de la longitude:
296        allocate(lon(lendim(idim)))
297        ierr = nf_inq_varid(orig_file_id, namedim, lonid)
298        if (ierr /= 0) call handle_err(ierr)
299        ierr = nf_get_var_real(orig_file_id, lonid, lon)
300        if (ierr /= 0) call handle_err(ierr)
301        ierr = nf_get_att_text(orig_file_id, lonid, 'units', units)
302        if (ierr /= 0) call handle_err(ierr)
303        ALLOCATE(lon_bounds(lendim(idim)+1))
304        DO i = 2, lendim(idim)
305         lon_bounds(i) = lon(i-1) - (lon(i-1) - lon(i))/2
306        enddo
307        lon_bounds(1) = lon(1) - (lon_bounds(3) -lon_bounds(2))/2.
308        lon_bounds(lendim(idim)+1) = lon(lendim(idim)) + (lon_bounds(lendim(idim))-lon_bounds(lendim(idim)-1))/2.
309
310!       definition de la longitude
311        axis_ids(idim) = cmor_axis(                           &
312                    table=trim(ipcc_table(index_table)),      &
313                    table_entry='longitude',                  &
314                    units=units,                              & 
315                    length=lendim(idim),                            &
316                    coord_vals=lon,                           &
317                    cell_bounds=lon_bounds)       
318!
319!
320      case('presnivs')
321!     lecture de la verticale:
322        allocate(vert(lendim(idim)))
323        ierr = nf_inq_varid(orig_file_id, namedim, vertid)
324        if (ierr /= 0) call handle_err(ierr)
325        ierr = nf_get_var_real(orig_file_id, vertid, vert)
326        if (ierr /= 0) call handle_err(ierr)
327        ierr = nf_get_att_text(orig_file_id, vertid, 'units', units)
328        if (ierr /= 0) call handle_err(ierr)
329!
330!       definition de la verticale
331        if (units == 'mb') units='Pa'
332        axis_ids(idim) = cmor_axis(                           &
333                    table=trim(ipcc_table(index_table)),      &
334                    table_entry='pressure',                   &
335                    units=units,                              & 
336                    length=lendim(idim),                      &
337                    coord_vals=vert)
338!
339!       
340      case('time_counter')         
341!     definition du temps
342      if (idim /= ndims) then
343        write(lunout,*)'la dimension temps doit etre la derniere dimension'
344        stop
345      endif
346      allocate(time(lendim(idim)))
347      ierr = nf_inq_varid(orig_file_id,namedim,timeid)
348      if (ierr /=0) call handle_err(ierr)
349      ierr = nf_get_var_real(orig_file_id, timeid, time)
350      if (ierr /=0) call handle_err(ierr)
351      ierr = nf_get_att_text(orig_file_id,timeid, 'units', units)
352      if (ierr /=0) call handle_err(ierr)
353      axis_ids(idim) = cmor_axis(                          &
354           table=trim(ipcc_table(index_table)),            &
355           table_entry='time',                             &
356           units=units,                                    &
357           length=lendim(idim),                            &
358           interval='30 minutes')
359      itime= lendim(idim)
360    case default
361      write(lunout,*)'Dimension: ', trim(namedim),' non reconnue'
362      stop
363   endselect
364  enddo
365 
366!
367! Definition de la variable a ecrire
368  units=' '
369  ierr = nf_inq_varid(orig_file_id,TRIM(varname), varid)
370  if (ierr /= 0) call handle_err(ierr)
371  ierr = nf_get_att_text(orig_file_id, varid, 'units', units)
372  if (ierr /= 0) call handle_err(ierr)
373  ierr = nf_get_att_real(orig_file_id, varid, 'missing_value', missing_value)
374  if (ierr /= 0) call handle_err(ierr)
375  cmorvarid = cmor_variable(                         &
376       table=trim(ipcc_table(index_table)),          &
377       table_entry=trim(ipcc_name(index_table)),     &
378       units=units,                                  &
379       axis_ids= axis_ids,                           &
380       missing_value=real(missing_value),            &
381       positive = trim(ipsl_pos(index_table)),       &
382       original_name=varname)
383!
384! Lecture de la variable
385  if (ndims == 3) then
386    ALLOCATE (donnees(lendim(1), lendim(2), 1, lendim(3) ))
387  else if (ndims ==4) then
388    allocate (donnees(lendim(1), lendim(2), lendim(3), lendim(4) ))
389  endif
390  ierr = nf_get_var_real(orig_file_id, varid, donnees)
391!
392! Ecriture de la variable
393 
394  DO i = 1, itime
395    rdate(1) = dble(time(i))
396    ierr = cmor_write(                                     &
397             var_id        = cmorvarid,                    &
398             DATA          = REAL(donnees(:,:,:,i)),       &
399             ntimes_passed = 1,                            &
400             time_vals     = rdate)
401  enddo
402!
403! Fin CMOR
404  ierr = cmor_close()
405  IF (ierr /= 0) then
406    WRITE(lunout,*)'Probleme dans cmor_close, ierr = ', ierr
407  endif
408
409!
410! fermeture fichier originel
411  ierr = nf_close(orig_file_id)
412  IF (ierr /= NF_NOERR) then
413    WRITE(lunout,*)NF_STRERROR(ierr)
414    stop
415  endif
416
417contains
418
419subroutine handle_err(status)
420  use netcdf
421
422  implicit none
423  integer, intent(in) :: status
424
425  write(lunout,*)nf90_strerror(status)
426  stop
427
428end subroutine handle_err
429 
430END PROGRAM ts2IPCC
431
Note: See TracBrowser for help on using the repository browser.