lfi2cdf: options are now stored in a structure, command line is read in Fortran
[MNH-git_open_source-lfs.git] / tools / lfi2cdf / src / mode_util.f90
1 MODULE mode_util
2   USE MODE_FIELDTYPE
3   USE mode_dimlist
4   USE mode_options
5   USE MODD_PARAM
6   USE netcdf
7
8   IMPLICIT NONE 
9
10   INTEGER,PARAMETER :: MAXRAW=10
11   INTEGER,PARAMETER :: MAXLEN=512
12   INTEGER,PARAMETER :: MAXFILES=100
13
14   INTEGER,PARAMETER :: UNDEFINED = -1, READING = 1, WRITING = 2
15   INTEGER,PARAMETER :: UNKNOWN_FORMAT = -1, NETCDF_FORMAT = 1, LFI_FORMAT = 2
16
17   TYPE filestruct
18     INTEGER :: lun_id                  ! Logical ID of file
19     INTEGER :: format = UNKNOWN_FORMAT ! NETCDF, LFI
20     INTEGER :: status = UNDEFINED      ! Opened for reading or writing
21     INTEGER :: var_id                  ! Position of the variable in the workfield structure
22     LOGICAL :: opened = .false.
23   END TYPE filestruct
24
25   TYPE filelist_struct
26     INTEGER :: nbfiles = 0
27 !    TYPE(filestruct),DIMENSION(:),ALLOCATABLE :: files
28     TYPE(filestruct),DIMENSION(MAXFILES) :: files
29   END TYPE filelist_struct
30
31
32   TYPE workfield
33      CHARACTER(LEN=FM_FIELD_SIZE)            :: name   ! nom du champ
34      INTEGER                                 :: TYPE   ! type (entier ou reel)    
35      CHARACTER(LEN=:), POINTER               :: comment
36      TYPE(dimCDF),                   POINTER :: dim
37      INTEGER                                 :: id_in = -1, id_out = -1
38      INTEGER                                 :: grid
39      LOGICAL                                 :: found  ! T if found in the input file
40      LOGICAL                                 :: calc   ! T if computed from other variables
41      LOGICAL                                 :: tbw    ! to be written or not
42      LOGICAL                                 :: tbr    ! to be read or not
43      INTEGER,DIMENSION(MAXRAW)               :: src    ! List of variables used to compute the variable (needed only if calc=.true.)
44      INTEGER                                 :: tgt    ! Target: id of the variable that use it (calc variable)
45   END TYPE workfield
46
47 #ifndef LOWMEM
48   TYPE lfidata
49      INTEGER(KIND=8), DIMENSION(:), POINTER :: iwtab
50   END TYPE lfidata
51   TYPE(lfidata), DIMENSION(:), ALLOCATABLE :: lfiart
52 #endif
53
54   LOGICAL(KIND=LFI_INT), PARAMETER :: ltrue  = .TRUE.
55   LOGICAL(KIND=LFI_INT), PARAMETER :: lfalse = .FALSE.
56
57 CONTAINS 
58   FUNCTION str_replace(hstr, hold, hnew)
59     CHARACTER(LEN=*) :: hstr, hold, hnew
60     CHARACTER(LEN=LEN_TRIM(hstr)+MAX(0,LEN(hnew)-LEN(hold))) :: str_replace
61     
62     INTEGER :: pos
63     
64     pos = INDEX(hstr,hold)
65     IF (pos /= 0) THEN
66        str_replace = hstr(1:pos-1)//hnew//hstr(pos+LEN(hold):)
67     ELSE 
68        str_replace = hstr 
69     END IF
70
71   END FUNCTION str_replace
72
73   SUBROUTINE FMREADLFIN1(klu,hrecfm,kval,kresp)
74   INTEGER(KIND=LFI_INT), INTENT(IN) :: klu ! logical fortran unit au lfi file
75   CHARACTER(LEN=*),INTENT(IN)       :: hrecfm ! article name to be read
76   INTEGER, INTENT(OUT)        :: kval ! integer value for hrecfm article
77   INTEGER(KIND=LFI_INT), INTENT(OUT):: kresp! return code null if OK
78   !
79   INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE::iwork
80   INTEGER :: icomlen
81   INTEGER(KIND=LFI_INT) :: iresp,ilenga,iposex
82   !
83   CALL LFINFO(iresp,klu,hrecfm,ilenga,iposex)
84   IF (iresp /=0 .OR. ilenga == 0) THEN
85     kresp = -1
86     kval = 0
87   ELSE
88     ALLOCATE(IWORK(ilenga))
89     CALL LFILEC(iresp,klu,hrecfm,iwork,ilenga)
90     icomlen = iwork(2)
91     kval = iwork(3+icomlen)
92     kresp = iresp
93     DEALLOCATE(IWORK)
94   END IF
95   END SUBROUTINE FMREADLFIN1
96
97   SUBROUTINE parse_infiles(infiles, nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw, tpreclist, kbuflen, options, icurrent_level)
98     TYPE(filelist_struct),      INTENT(IN) :: infiles
99     INTEGER,                    INTENT(IN) :: nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw
100     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
101     INTEGER,                   INTENT(OUT) :: kbuflen
102     TYPE(option),DIMENSION(:), INTENT(IN)  :: options
103     INTEGER,          INTENT(IN), OPTIONAL :: icurrent_level
104
105     INTEGER                                  :: ji,jj, kcdf_id, itype
106     INTEGER                                  :: ndb, nde, ndey, idx, idx_var, maxvar
107     INTEGER                                  :: idims, idimtmp, jdim, status, var_id
108     LOGICAL                                  :: ladvan
109     INTEGER                                  :: ich, current_level, leng
110     INTEGER                                  :: comment_size, fsize, sizemax
111     CHARACTER(LEN=FM_FIELD_SIZE)             :: yrecfm
112     CHARACTER(LEN=4)                         :: suffix
113 #ifdef LOWMEM
114     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
115 #endif
116     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
117     CHARACTER(LEN=FM_FIELD_SIZE)             :: var_calc
118     CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
119     INTEGER, DIMENSION(10)                   :: idim_id
120     !JUAN CYCCL3
121     INTEGER                        :: JPHEXT
122
123     IF (infiles%files(1)%format == LFI_FORMAT) THEN
124       ilu = infiles%files(1)%lun_id
125
126       CALL FMREADLFIN1(ilu,'JPHEXT',JPHEXT,iresp)
127       IF (iresp /= 0) JPHEXT=1
128
129       ! First check if IMAX,JMAX,KMAX exist in LFI file
130       ! to handle 3D, 2D variables -> update IDIMX,IDIMY,IDIMZ
131       CALL FMREADLFIN1(ilu,'IMAX',IDIMX,iresp)
132       IF (iresp == 0) IDIMX = IDIMX+2*JPHEXT  ! IMAX + 2*JPHEXT
133        !
134       CALL FMREADLFIN1(ilu,'JMAX',IDIMY,iresp)
135       IF (iresp == 0) IDIMY = IDIMY+2*JPHEXT  ! JMAX + 2*JPHEXT
136       !
137       CALL FMREADLFIN1(ilu,'KMAX',IDIMZ,iresp)
138       IF (iresp == 0) IDIMZ = IDIMZ+2  ! KMAX + 2*JPVEXT
139     ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
140       kcdf_id = infiles%files(1)%lun_id
141
142       status = NF90_INQ_DIMID(kcdf_id, "DIMX", idim_id(1))
143       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
144       status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(1),len = IDIMX)
145
146       status = NF90_INQ_DIMID(kcdf_id, "DIMY", idim_id(2))
147       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
148       status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(2),len = IDIMY)
149
150       status = NF90_INQ_DIMID(kcdf_id, "DIMZ", idim_id(3))
151       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
152       status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(3),len = IDIMZ)
153     END IF
154
155     GUSEDIM = (IDIMX*IDIMY > 0)
156     IF (GUSEDIM) THEN
157       PRINT *,'MESONH 3D, 2D articles DIMENSIONS used :'
158       PRINT *,'DIMX =',IDIMX
159       PRINT *,'DIMY =',IDIMY
160       PRINT *,'DIMZ =',IDIMZ ! IDIMZ may be equal to 0 (PGD files)
161     ELSE
162       PRINT *,'BEWARE : ALL MesoNH arrays are handled as 1D arrays !'
163     END IF
164
165     sizemax = 0
166
167     IF (present(icurrent_level)) THEN
168       write(suffix,'(I4.4)') icurrent_level
169       current_level = icurrent_level
170     ElSE
171       suffix=''
172       current_level = -1
173     END IF
174
175     ! Phase 1 : build articles list to convert.
176     !
177     !    Pour l'instant tous les articles du fichier LFI sont
178     !    convertis. On peut modifier cette phase pour prendre en
179     !    compte un sous-ensemble d'article (liste definie par
180     !    l'utilisateur par exemple)
181     !
182     IF (options(OPTVAR)%set) THEN
183 #ifndef LOWMEM
184       IF(.NOT.ALLOCATED(lfiart) .AND. infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(nbvar_tbr+nbvar_calc))
185 #endif
186       ALLOCATE(tpreclist(nbvar_tbr+nbvar_calc))
187       DO ji=1,nbvar_tbr+nbvar_calc
188         tpreclist(ji)%found  = .FALSE.
189         tpreclist(ji)%calc   = .FALSE. !By default variables are not computed from others
190         tpreclist(ji)%tbw    = .TRUE.  !By default variables are written
191         tpreclist(ji)%tbr    = .TRUE.  !By default variables are written
192         tpreclist(ji)%src(:) = -1
193         tpreclist(ji)%tgt    = -1
194       END DO
195
196        ! A variable list is provided with -v var1,...
197        ndb  = 1
198        idx_var = 1
199        DO ji=1,nbvar_tbw
200           !crash compiler GCC 4.2.0: nde = INDEX(TRIM(options(OPTVAR)%cvalue(ndb:)),',')
201           nde = INDEX(TRIM(options(OPTVAR)%cvalue(ndb:len(trim(options(OPTVAR)%cvalue)))),',')
202           yrecfm = options(OPTVAR)%cvalue(ndb:ndb+nde-2)
203
204           !Detect operations on variables (only + is supported now)
205           ndey = INDEX(TRIM(yrecfm),'=')
206           idx = 1
207           IF (ndey /= 0) THEN
208             var_calc = yrecfm(1:ndey-1)
209             DO WHILE (ndey /= 0)
210               IF (idx>MAXRAW) THEN
211                 print *,'Error: MAXRAW exceeded (too many raw variables for 1 computed one)'
212                 STOP
213               END IF
214               yrecfm = yrecfm(ndey+1:)
215               ndey = INDEX(TRIM(yrecfm),'+')
216               IF (ndey /= 0) THEN
217                 var_raw(idx) = yrecfm(1:ndey-1)
218               ELSE
219                 var_raw(idx) = TRIM(yrecfm)
220               END IF
221               idx = idx + 1
222             END DO
223
224             tpreclist(idx_var)%name = trim(var_calc)
225             tpreclist(idx_var)%calc = .TRUE.
226             tpreclist(idx_var)%tbw  = .TRUE.
227             tpreclist(idx_var)%tbr  = .FALSE.
228             idx_var=idx_var+1
229             DO jj = 1, idx-1
230               tpreclist(idx_var-jj)%src(jj) = idx_var
231               tpreclist(idx_var)%name = trim(var_raw(jj))
232               tpreclist(idx_var)%calc = .FALSE.
233               tpreclist(idx_var)%tbw  = .FALSE.
234               tpreclist(idx_var)%tbr  = .TRUE.
235               tpreclist(idx_var)%tgt  = idx_var-jj
236               idx_var=idx_var+1
237             END DO
238
239           ELSE
240             tpreclist(idx_var)%name = trim(yrecfm)
241             tpreclist(idx_var)%calc = .FALSE.
242             tpreclist(idx_var)%tbw  = .TRUE.
243             idx_var=idx_var+1
244
245           END IF
246
247           ndb = nde+ndb
248        END DO
249
250        DO ji=1,nbvar_tbr+nbvar_calc
251           IF (tpreclist(ji)%calc) CYCLE
252
253           yrecfm = TRIM(tpreclist(ji)%name)
254           IF (infiles%files(1)%format == LFI_FORMAT) THEN
255             CALL LFINFO(iresp,ilu,trim(yrecfm)//trim(suffix),ileng,ipos)
256             IF (iresp == 0 .AND. ileng /= 0) tpreclist(ji)%found = .true.
257             leng = ileng
258           ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
259             status = NF90_INQ_VARID(kcdf_id,trim(yrecfm)//trim(suffix),tpreclist(ji)%id_in)
260             IF (status == NF90_NOERR) THEN
261               tpreclist(ji)%found = .true.
262               status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in,ndims = idims,dimids = idim_id)
263               IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
264
265 !TODO:useful?
266 !DUPLICATED
267               IF (idims == 0) THEN
268                  ! variable scalaire
269                  leng = 1
270               ELSE
271                  ! infos sur dimensions
272                  leng = 1
273                  DO jdim=1,idims
274                    status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
275                    IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
276                    leng = leng*idimtmp
277                 END DO
278               END IF
279             END IF
280           END IF
281
282           IF (.NOT.tpreclist(ji)%found) THEN
283              PRINT *,'Article ',TRIM(yrecfm), ' not found!'
284              tpreclist(ji)%tbw   = .FAlSE.
285              tpreclist(ji)%tbr   = .FAlSE.
286           ELSE
287              ! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
288              IF (leng > sizemax) sizemax = leng
289 #ifndef LOWMEM
290 !TODO:useful for netcdf?
291              IF (infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(ji)%iwtab(leng))
292 #endif
293           END IF
294        END DO
295
296        maxvar = nbvar_tbr+nbvar_calc
297
298 DO ji=1,nbvar_tbr+nbvar_calc
299   print *,ji,'name=',trim(tpreclist(ji)%name),' calc=',tpreclist(ji)%calc,' tbw=',tpreclist(ji)%tbw,&
300           ' tbr=',tpreclist(ji)%tbr,' found=',tpreclist(ji)%found
301 END DO
302
303     ELSE
304        ! Entire file is converted
305 #ifndef LOWMEM
306        IF(.NOT.ALLOCATED(lfiart) .AND. infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(nbvar_infile))
307 #endif
308        ALLOCATE(tpreclist(nbvar_infile))
309        DO ji=1,nbvar_infile
310          tpreclist(ji)%calc   = .FALSE. !By default variables are not computed from others
311          tpreclist(ji)%tbw    = .TRUE.  !By default variables are written
312          tpreclist(ji)%src(:) = -1
313        END DO
314
315        IF (infiles%files(1)%format == LFI_FORMAT) THEN
316          CALL LFIPOS(iresp,ilu)
317          ladvan = .TRUE.
318
319          DO ji=1,nbvar_infile
320            CALL LFICAS(iresp,ilu,yrecfm,ileng,ipos,ladvan)
321            ! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
322            tpreclist(ji)%name = trim(yrecfm)
323            tpreclist(ji)%found  = .TRUE.
324            IF (ileng > sizemax) sizemax = ileng
325 #ifndef LOWMEM
326            ALLOCATE(lfiart(ji)%iwtab(ileng))
327 #endif
328          END DO
329        ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
330          DO ji=1,nbvar_infile
331            tpreclist(ji)%id_in = ji
332            status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in, name = tpreclist(ji)%name, ndims = idims, &
333                                           dimids = idim_id)
334            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
335            ! PRINT *,'Article ',ji,' : ',TRIM(tpreclist(ji)%name),', longueur = ',ileng
336            tpreclist(ji)%found  = .TRUE.
337 !TODO:useful?
338 !DUPLICATED
339            IF (idims == 0) THEN
340              ! variable scalaire
341              leng = 1
342            ELSE
343              ! infos sur dimensions
344              leng = 1
345              DO jdim=1,idims
346                status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
347                IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
348                leng = leng*idimtmp
349              END DO
350            END IF
351            IF (leng > sizemax) sizemax = leng
352          END DO
353        END IF
354
355        maxvar = nbvar_infile
356     END IF
357
358     kbuflen = sizemax
359
360 #ifdef LOWMEM
361     WRITE(*,'("Taille maximale du buffer :",f10.3," Mio")') sizemax*8./1048576.
362     ALLOCATE(iwork(sizemax))
363 #endif
364
365     ! Phase 2 : Extract comments and dimensions for valid articles.
366     !           Infos are put in tpreclist.
367     CALL init_dimCDF()
368     DO ji=1,maxvar
369        IF (tpreclist(ji)%calc .OR. .NOT.tpreclist(ji)%found) CYCLE
370
371        IF (infiles%files(1)%format == LFI_FORMAT) THEN
372          yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
373          CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
374 #ifdef LOWMEM
375          CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
376          tpreclist(ji)%grid = iwork(1)
377          comment_size = iwork(2)
378 #else
379          CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
380          tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
381          comment_size = lfiart(ji)%iwtab(2)
382 #endif
383          tpreclist(ji)%TYPE = get_ftype(yrecfm,current_level)
384
385          ALLOCATE(character(len=comment_size) :: tpreclist(ji)%comment)
386          DO jj=1,comment_size
387 #ifdef LOWMEM
388            ich = iwork(2+jj)
389 #else
390            ich = lfiart(ji)%iwtab(2+jj)
391 #endif
392            tpreclist(ji)%comment(jj:jj) = CHAR(ich)
393          END DO
394
395          fsize = ileng-(2+comment_size)
396
397        ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
398          ! GRID attribute definition
399          status = NF90_GET_ATT(kcdf_id,tpreclist(ji)%id_in,'GRID',tpreclist(ji)%grid)
400          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
401
402          ! COMMENT attribute definition
403          status = NF90_INQUIRE_ATTRIBUTE(kcdf_id,tpreclist(ji)%id_in,'COMMENT',len=comment_size)
404          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
405          ALLOCATE(character(len=comment_size) :: tpreclist(ji)%comment)
406          status = NF90_GET_ATT(kcdf_id,tpreclist(ji)%id_in,'COMMENT',tpreclist(ji)%comment)
407          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
408
409          status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in, xtype = itype, ndims = idims, &
410                                         dimids = idim_id)
411          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
412
413          SELECT CASE(itype)
414          CASE(NF90_CHAR)
415            tpreclist(ji)%TYPE = TEXT
416          CASE(NF90_INT)
417            tpreclist(ji)%TYPE = INT
418          CASE(NF90_FLOAT,NF90_DOUBLE)
419            tpreclist(ji)%TYPE = FLOAT
420          CASE default
421            PRINT *, 'Attention : variable ',TRIM(tpreclist(ji)%name), ' a un TYPE non reconnu par le convertisseur.'
422            PRINT *, '--> TYPE force a REAL(KIND 8) dans LFI !'
423          END SELECT
424
425 !DUPLICATED
426          IF (idims == 0) THEN
427            ! variable scalaire
428            leng = 1
429          ELSE
430            ! infos sur dimensions
431            leng = 1
432            DO jdim=1,idims
433              status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
434              IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
435              leng = leng*idimtmp
436            END DO
437          END IF
438
439          fsize = leng
440        END IF
441
442        tpreclist(ji)%dim=>get_dimCDF(fsize)
443     END DO
444
445     !Complete info for calculated variables
446     IF (nbvar_calc>0) THEN
447     DO ji=1,maxvar
448        IF (.NOT.tpreclist(ji)%calc) CYCLE
449        tpreclist(ji)%TYPE = tpreclist(tpreclist(ji)%src(1))%TYPE
450        tpreclist(ji)%grid = tpreclist(tpreclist(ji)%src(1))%grid
451        tpreclist(ji)%dim  => tpreclist(tpreclist(ji)%src(1))%dim
452
453 !TODO: cleaner length!
454        ALLOCATE(character(len=256) :: tpreclist(ji)%comment)
455        tpreclist(ji)%comment='Constructed from'
456        jj = 1
457        DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
458          tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' '//trim(tpreclist(tpreclist(ji)%src(jj))%name)
459          IF (jj<MAXRAW .AND. tpreclist(ji)%src(jj+1)>0) THEN
460            tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' +'
461          END IF
462          jj=jj+1
463        END DO
464     END DO
465     END IF
466
467
468     PRINT *,'Nombre de dimensions = ', size_dimCDF()
469 #ifdef LOWMEM
470     DEALLOCATE(iwork)
471 #endif
472   END SUBROUTINE parse_infiles
473   
474   SUBROUTINE read_data_lfi(infiles, nbvar, tpreclist, kbuflen, current_level)
475     TYPE(filelist_struct),      INTENT(IN) :: infiles
476     INTEGER, INTENT(INOUT)                 :: nbvar
477     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
478     INTEGER, INTENT(IN)                    :: kbuflen
479     INTEGER, INTENT(IN), OPTIONAL          :: current_level
480
481     INTEGER                                  :: ji,jj
482     INTEGER                                  :: ndb, nde
483     LOGICAL                                  :: ladvan
484     INTEGER                                  :: ich
485     INTEGER                                  :: fsize,sizemax
486     CHARACTER(LEN=FM_FIELD_SIZE)             :: yrecfm
487     CHARACTER(LEN=4)                         :: suffix
488 #ifdef LOWMEM
489     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
490 #endif
491     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
492     CHARACTER(LEN=FM_FIELD_SIZE)             :: var_calc
493     CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
494
495
496     ilu = infiles%files(1)%lun_id
497
498     IF (present(current_level)) THEN
499       write(suffix,'(I4.4)') current_level
500     ElSE
501       suffix=''
502     END IF
503
504 #ifdef LOWMEM
505     ALLOCATE(iwork(kbuflen))
506 #endif
507
508     DO ji=1,nbvar
509        IF (.NOT.tpreclist(ji)%tbr) CYCLE
510        yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
511        CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
512 #ifdef LOWMEM
513        CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
514        tpreclist(ji)%grid = iwork(1)
515 #else
516        CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
517        tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
518 #endif
519     END DO
520
521 #ifdef LOWMEM
522     DEALLOCATE(iwork)
523 #endif
524   END SUBROUTINE read_data_lfi
525
526   SUBROUTINE HANDLE_ERR(status,line)
527     INTEGER :: status,line
528
529     IF (status /= NF90_NOERR) THEN
530        PRINT *, 'line ',line,': ',NF90_STRERROR(status)
531            STOP
532     END IF
533   END SUBROUTINE HANDLE_ERR
534
535   SUBROUTINE def_ncdf(outfiles,tpreclist,nbvar,options)
536     TYPE(filelist_struct),       INTENT(IN) :: outfiles
537     TYPE(workfield),DIMENSION(:),INTENT(INOUT) :: tpreclist
538     INTEGER,                     INTENT(IN) :: nbvar
539     TYPE(option),DIMENSION(:),   INTENT(IN) :: options
540
541     INTEGER :: compress_level, status
542     INTEGER :: idx, ji, nbfiles
543     INTEGER:: kcdf_id
544     TYPE(dimCDF), POINTER :: tzdim
545     INTEGER               :: invdims
546     INTEGER               :: type_float
547     INTEGER, DIMENSION(10) :: ivdims
548     CHARACTER(LEN=20)     :: ycdfvar
549
550
551     nbfiles = outfiles%nbfiles
552
553     IF (options(OPTREDUCE)%set) THEN
554       type_float = NF90_REAL
555     ELSE
556       type_float = NF90_DOUBLE
557     END IF
558
559     DO ji = 1,nbfiles
560       kcdf_id = outfiles%files(ji)%lun_id
561
562       ! global attributes
563       status = NF90_PUT_ATT(kcdf_id,NF90_GLOBAL,'Title',VERSION_ID)
564       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
565
566       ! define DIMENSIONS
567       tzdim=>first_DimCDF()
568       DO WHILE(ASSOCIATED(tzdim))
569         IF (tzdim%create) THEN
570           status = NF90_DEF_DIM(kcdf_id,tzdim%name,tzdim%len,tzdim%id)
571           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
572         END IF
573         tzdim=>tzdim%next
574       END DO
575     END DO
576
577     PRINT *,'------------- NetCDF DEFINITION ---------------'
578
579     ! define VARIABLES and ATTRIBUTES
580     idx = 1
581     DO ji=1,nbvar
582        IF (.NOT.tpreclist(ji)%tbw) CYCLE
583
584        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
585          IF (tpreclist(ji)%dim%create) THEN
586            invdims   = 1
587            ivdims(1) = tpreclist(ji)%dim%id
588          ELSE
589            invdims = tpreclist(ji)%dim%ndims
590            IF(options(OPTMERGE)%set) invdims=invdims+1 !when merging variables from LFI splitted files
591            SELECT CASE(invdims)
592            CASE(2)
593               ivdims(1)=ptdimx%id
594               ivdims(2)=ptdimy%id
595            CASE(3)
596               ivdims(1)=ptdimx%id
597               ivdims(2)=ptdimy%id
598               ivdims(3)=ptdimz%id
599            CASE(12)
600               ivdims(1)=ptdimx%id
601               ivdims(2)=ptdimz%id
602               invdims = 2 ! on retablit la bonne valeur du nbre de dimension
603            CASE default
604              PRINT *,'Fatal error in NetCDF dimension definition'
605              STOP
606            END SELECT
607          END IF
608        ELSE
609          ! scalar variables
610           invdims   = 0
611           ivdims(1) = 0 ! ignore dans ce cas
612        END IF
613        
614        ! Variables definition
615
616        !! NetCDF n'aime pas les '%' dans le nom des variables
617        !! "%" remplaces par '__' 
618        ycdfvar = str_replace(tpreclist(ji)%name,'%','__')
619        !! ni les '.' remplaces par '--'
620        ycdfvar = str_replace(ycdfvar,'.','--')
621
622        kcdf_id = outfiles%files(idx)%lun_id
623
624        SELECT CASE(tpreclist(ji)%TYPE)
625        CASE (TEXT)
626 !          PRINT *,'TEXT : ',tpreclist(ji)%name
627           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_CHAR,&
628                    ivdims(:invdims),tpreclist(ji)%id_out)
629           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
630
631        CASE (INT,BOOL)
632 !          PRINT *,'INT,BOOL : ',tpreclist(ji)%name
633           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_INT,&
634                    ivdims(:invdims),tpreclist(ji)%id_out)
635           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
636
637        CASE(FLOAT)
638 !          PRINT *,'FLOAT : ',tpreclist(ji)%name
639           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
640                    ivdims(:invdims),tpreclist(ji)%id_out)
641           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
642
643           
644        CASE default
645           PRINT *,'ATTENTION : ',TRIM(tpreclist(ji)%name),' est de&
646                & TYPE inconnu --> force a REAL'
647           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
648                    ivdims(:invdims),tpreclist(ji)%id_out)
649           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
650           
651
652        END SELECT
653
654        ! Compress data (costly operation for the CPU)
655        IF (options(OPTCOMPRESS)%set .AND. invdims>0) THEN
656          compress_level = options(OPTCOMPRESS)%ivalue
657          status = NF90_DEF_VAR_DEFLATE(kcdf_id,tpreclist(ji)%id_out,1,1,compress_level)
658          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
659        END IF
660
661        ! GRID attribute definition
662        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'GRID',tpreclist(ji)%grid)
663        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
664
665        ! COMMENT attribute definition
666        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'COMMENT',trim(tpreclist(ji)%comment))
667        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
668
669        IF (options(OPTSPLIT)%set) idx = idx + 1
670     END DO
671     
672     DO ji = 1,nbfiles
673       kcdf_id = outfiles%files(ji)%lun_id
674       status = NF90_ENDDEF(kcdf_id)
675       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
676     END DO
677
678   END SUBROUTINE def_ncdf
679
680   SUBROUTINE fill_ncdf(infiles,outfiles,tpreclist,knaf,kbuflen,options,current_level)
681     TYPE(filelist_struct),        INTENT(IN):: infiles, outfiles
682     TYPE(workfield), DIMENSION(:),INTENT(IN):: tpreclist
683     INTEGER,                      INTENT(IN):: knaf
684     INTEGER,                      INTENT(IN):: kbuflen
685     TYPE(option),DIMENSION(:),    INTENT(IN):: options
686     INTEGER, INTENT(IN), OPTIONAL           :: current_level
687
688 #ifdef LOWMEM
689     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
690 #endif
691     INTEGER                                  :: idx, ji,jj
692     INTEGER                                  :: kcdf_id
693     INTEGER                                  :: status
694     INTEGER                                  :: extent, ndims
695     INTEGER                                  :: ich
696     INTEGER                                  :: src
697     INTEGER                                  :: level
698     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
699     CHARACTER(LEN=4)                         :: suffix
700     INTEGER,DIMENSION(3)                     :: idims, start
701     INTEGER,DIMENSION(:),ALLOCATABLE         :: itab
702     REAL(KIND=8),DIMENSION(:),ALLOCATABLE    :: xtab
703     CHARACTER, DIMENSION(:), ALLOCATABLE     :: ytab
704     REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: xtab3d, xtab3d2
705     INTEGER,      DIMENSION(:,:,:), ALLOCATABLE :: itab3d, itab3d2
706
707
708     !
709     IF (infiles%files(1)%format == LFI_FORMAT) ilu = infiles%files(1)%lun_id
710     !
711
712     IF (present(current_level)) THEN
713       write(suffix,'(I4.4)') current_level
714       level = current_level
715     ElSE
716       suffix=''
717       level = 1
718     END IF
719
720 #if LOWMEM
721     ALLOCATE(iwork(kbuflen))
722 #endif
723     ALLOCATE(itab(kbuflen))
724     ALLOCATE(xtab(kbuflen))
725
726     idx = 1
727     DO ji=1,knaf
728        IF (.NOT.tpreclist(ji)%tbw) CYCLE
729
730        kcdf_id = outfiles%files(idx)%lun_id
731
732        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
733           extent = tpreclist(ji)%dim%len
734           ndims = tpreclist(ji)%dim%ndims
735        ELSE
736           extent = 1
737           ndims = 0
738        END IF
739
740        idims(:) = 1
741        if(ndims>0) idims(1) = ptdimx%len
742        if(ndims>1) idims(2) = ptdimy%len
743        if(ndims>2) idims(3) = ptdimz%len
744        if(ndims>3) then
745          PRINT *,'Too many dimensions'
746          STOP
747        endif
748
749        SELECT CASE(tpreclist(ji)%TYPE)
750        CASE (INT,BOOL)
751         IF (infiles%files(1)%format == LFI_FORMAT) THEN
752 #if LOWMEM
753          IF (.NOT.tpreclist(ji)%calc) THEN
754            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
755            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
756            itab(1:extent) = iwork(3+iwork(2):)
757          ELSE
758            src=tpreclist(ji)%src(1)
759            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
760            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
761            itab(1:extent) = iwork(3+iwork(2):)
762            jj = 2
763            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
764              src=tpreclist(ji)%src(jj)
765              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
766              itab(1:extent) = itab(1:extent) + iwork(3+iwork(2):)
767              jj=jj+1
768            END DO
769          ENDIF
770 #else
771          IF (.NOT.tpreclist(ji)%calc) THEN
772            itab(1:extent) = lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):)
773          ELSE
774            src=tpreclist(ji)%src(1)
775            itab(1:extent) = lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
776            jj = 2
777            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
778              src=tpreclist(ji)%src(jj)
779              itab(1:extent) = xtab(1:extent) + lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
780              jj=jj+1
781            END DO
782          END IF
783 #endif
784
785 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
786          SELECT CASE(ndims)
787          CASE (0)
788            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1))
789          CASE (1)
790            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1:extent),count=(/extent/))
791          CASE (2)
792            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len/)), &
793                                  start = (/1,1,level/) )
794          CASE (3)
795            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
796          CASE DEFAULT
797            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
798          END SELECT
799
800         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
801          ALLOCATE( itab3d(idims(1),idims(2),idims(3)) )
802          IF (.NOT.tpreclist(ji)%calc) THEN
803            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,itab3d)
804            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
805          ELSE
806            ALLOCATE( itab3d2(idims(1),idims(2),idims(3)) )
807            src=tpreclist(ji)%src(1)
808            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d)
809            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
810            jj = 2
811            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
812              src=tpreclist(ji)%src(jj)
813              status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d2)
814              IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
815              itab3d(:,:,:) = itab3d(:,:,:) + itab3d2(:,:,:)
816              jj=jj+1
817            END DO
818            DEALLOCATE(itab3d2)
819          END IF
820
821 !TODO: not clean, should be done only if merging z-levels
822          IF (ndims == 2) THEN
823            start = (/1,1,level/)
824          ELSE
825            start = (/1,1,1/)
826          ENDIF
827          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab3d,start=start)
828          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
829
830          DEALLOCATE(itab3d)
831         END IF
832
833          
834        CASE (FLOAT)
835         IF (infiles%files(1)%format == LFI_FORMAT) THEN
836 #if LOWMEM
837          IF (.NOT.tpreclist(ji)%calc) THEN
838            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
839            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
840            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
841          ELSE
842            src=tpreclist(ji)%src(1)
843            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
844            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
845            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
846            jj = 2
847            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
848              src=tpreclist(ji)%src(jj)
849              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
850              xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
851              jj=jj+1
852            END DO
853          ENDIF
854 #else
855          IF (.NOT.tpreclist(ji)%calc) THEN
856            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
857          ELSE
858            src=tpreclist(ji)%src(1)
859            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
860            jj = 2
861            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
862              src=tpreclist(ji)%src(jj)
863              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
864              jj=jj+1
865            END DO
866          END IF
867 #endif
868 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
869          SELECT CASE(ndims)
870          CASE (0)
871            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
872          CASE (1)
873            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
874          CASE (2)
875            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
876                                  start = (/1,1,level/) )
877          CASE (3)
878            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
879          CASE DEFAULT
880            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
881          END SELECT
882
883         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
884          ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
885          IF (.NOT.tpreclist(ji)%calc) THEN
886            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,xtab3d)
887            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
888          ELSE
889            ALLOCATE( xtab3d2(idims(1),idims(2),idims(3)) )
890            src=tpreclist(ji)%src(1)
891            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d)
892            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
893            jj = 2
894            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
895              src=tpreclist(ji)%src(jj)
896              status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d2)
897              IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
898              xtab3d(:,:,:) = xtab3d(:,:,:) + xtab3d2(:,:,:)
899              jj=jj+1
900            END DO
901            DEALLOCATE(xtab3d2)
902          END IF
903
904 !TODO: not clean, should be done only if merging z-levels
905          IF (ndims == 2) THEN
906            start = (/1,1,level/)
907          ELSE
908            start = (/1,1,1/)
909          ENDIF
910          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab3d,start=start)
911          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
912
913          DEALLOCATE(xtab3d)
914         END IF
915
916        CASE (TEXT)
917         IF (infiles%files(1)%format == LFI_FORMAT) THEN
918          ALLOCATE(ytab(extent))
919          DO jj=1,extent
920 #if LOWMEM
921            ich = iwork(2+iwork(2)+jj)
922 #else
923            ich = lfiart(ji)%iwtab(2+lfiart(ji)%iwtab(2)+jj)
924 #endif
925            ytab(jj) = CHAR(ich)
926          END DO
927          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
928          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
929          DEALLOCATE(ytab)
930         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
931          status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,ytab,count=(/extent/))
932          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
933
934          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
935          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
936         END IF
937
938        CASE default
939         IF (infiles%files(1)%format == LFI_FORMAT) THEN
940 #if LOWMEM
941          IF (.NOT.tpreclist(ji)%calc) THEN
942            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
943            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
944            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
945          ELSE
946            src=tpreclist(ji)%src(1)
947            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
948            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
949            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
950            jj = 2
951            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
952              src=tpreclist(ji)%src(jj)
953              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
954              xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
955              jj=jj+1
956            END DO
957          ENDIF
958 #else         
959          IF (.NOT.tpreclist(ji)%calc) THEN
960            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
961          ELSE
962            src=tpreclist(ji)%src(1)
963            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
964            jj = 2
965            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
966              src=tpreclist(ji)%src(jj)
967              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
968              jj=jj+1
969            END DO
970          END IF
971 #endif
972 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
973          SELECT CASE(ndims)
974          CASE (0)
975            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
976          CASE (1)
977            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
978          CASE (2)
979            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
980                                  start = (/1,1,level/) )
981          CASE (3)
982            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
983          CASE DEFAULT
984            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
985          END SELECT
986          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
987         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
988          print *,'Error: unknown datatype'
989          STOP
990         END IF
991
992        END SELECT
993
994        if (options(OPTSPLIT)%set) idx = idx + 1
995     END DO
996     DEALLOCATE(itab,xtab)
997 #if LOWMEM
998     DEALLOCATE(iwork)
999 #endif 
1000   END SUBROUTINE fill_ncdf
1001
1002   SUBROUTINE build_lfi(infiles,outfiles,tpreclist,kbuflen)
1003     TYPE(filelist_struct),         INTENT(IN) :: infiles, outfiles
1004     TYPE(workfield), DIMENSION(:), INTENT(IN) :: tpreclist
1005     INTEGER,                       INTENT(IN) :: kbuflen
1006     
1007     INTEGER :: kcdf_id, status
1008     INTEGER :: ivar,ji,jj,ndims
1009     INTEGER,DIMENSION(3) :: idims
1010     INTEGER(KIND=8), DIMENSION(:), POINTER  :: iwork
1011     INTEGER(KIND=8), DIMENSION(:), POINTER  :: idata
1012     REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: xtab3d
1013     INTEGER,      DIMENSION(:,:,:), ALLOCATABLE :: itab3d
1014     CHARACTER,    DIMENSION(:), ALLOCATABLE :: ytab
1015     CHARACTER(LEN=FM_FIELD_SIZE)            :: yrecfm
1016
1017     INTEGER :: iartlen, idlen, icomlen
1018     INTEGER(KIND=LFI_INT) :: iresp,ilu,iartlen8
1019
1020
1021     ilu = outfiles%files(1)%lun_id
1022     kcdf_id = infiles%files(1)%lun_id
1023
1024     ! Un article LFI est compose de :
1025     !   - 1 entier identifiant le numero de grille
1026     !   - 1 entier contenant la taille du commentaire
1027     !   - le commentaire code en entier 64 bits
1028     !   - les donnees proprement dites
1029
1030     PRINT *,'Taille buffer = ',2+kbuflen
1031
1032     ALLOCATE(iwork(2+kbuflen))
1033
1034     DO ivar=1,SIZE(tpreclist)
1035        icomlen = LEN(tpreclist(ivar)%comment)
1036
1037        ! traitement Grille et Commentaire
1038        iwork(1) = tpreclist(ivar)%grid
1039        iwork(2) = icomlen
1040        DO jj=1,iwork(2)
1041           iwork(2+jj)=ICHAR(tpreclist(ivar)%comment(jj:jj))
1042        END DO
1043
1044        IF (ASSOCIATED(tpreclist(ivar)%dim)) THEN
1045           idlen = tpreclist(ivar)%dim%len
1046           ndims = tpreclist(ivar)%dim%ndims
1047        ELSE 
1048           idlen = 1
1049           ndims = 0
1050        END IF
1051        
1052        idims(:) = 1
1053        if(ndims>0) idims(1) = ptdimx%len
1054        if(ndims>1) idims(2) = ptdimy%len
1055        if(ndims>2) idims(3) = ptdimz%len
1056        if(ndims>3) then
1057          PRINT *,'Too many dimensions'
1058          STOP
1059        endif
1060
1061        iartlen = 2+icomlen+idlen
1062        idata=>iwork(3+icomlen:iartlen)
1063
1064
1065        SELECT CASE(tpreclist(ivar)%TYPE)
1066        CASE(INT,BOOL)
1067           ALLOCATE( itab3d(idims(1),idims(2),idims(3)) )
1068           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,itab3d)
1069           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1070
1071 !          PRINT *,'INT,BOOL --> ',tpreclist(ivar)%name,',len = ',idlen
1072           idata(1:idlen) = RESHAPE( itab3d , (/ idims(1)*idims(2)*idims(3) /) )
1073
1074           DEALLOCATE(itab3d)
1075
1076        CASE(FLOAT)
1077           ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
1078           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,xtab3d)
1079           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1080
1081 !          PRINT *,'FLOAT -->    ',tpreclist(ivar)%name,',len = ',idlen
1082           idata(1:idlen) = RESHAPE( TRANSFER(xtab3d,(/ 0_8 /),idlen) , (/ idims(1)*idims(2)*idims(3) /) )
1083
1084           DEALLOCATE(xtab3d)
1085
1086        CASE(TEXT)
1087           ALLOCATE(ytab(idlen))
1088           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,ytab)
1089           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1090
1091 !          PRINT *,'TEXT -->     ',tpreclist(ivar)%name,',len = ',idlen
1092           DO jj=1,idlen
1093              idata(jj) = ICHAR(ytab(jj))
1094           END DO
1095           
1096           DEALLOCATE(ytab)
1097
1098        CASE default
1099           ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
1100           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,xtab3d)
1101           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1102
1103           PRINT *,'Default (ERROR) -->',tpreclist(ivar)%name,',len = ',idlen
1104           idata(1:idlen) = RESHAPE( TRANSFER(xtab3d,(/ 0_8 /),idlen) , (/ idims(1)*idims(2)*idims(3) /) )
1105
1106           DEALLOCATE(xtab3d)
1107
1108        END SELECT
1109        
1110        ! Attention restoration des '%' dans le nom des champs LFI
1111        yrecfm = str_replace(tpreclist(ivar)%name,'__','%')
1112        ! et des '.'
1113        yrecfm = str_replace(yrecfm,'--','.')
1114        iartlen8 = iartlen
1115        CALL LFIECR(iresp,ilu,yrecfm,iwork,iartlen8)
1116
1117     END DO
1118     DEALLOCATE(iwork)
1119
1120   END SUBROUTINE build_lfi
1121
1122   SUBROUTINE UPDATE_VARID_IN(infiles,hinfile,tpreclist,nbvar,current_level)
1123     !Update the id_in for netCDF files (could change from one file to the other)
1124     TYPE(filelist_struct),         INTENT(IN)    :: infiles
1125     CHARACTER(LEN=*),              INTENT(IN)    :: hinfile
1126     TYPE(workfield), DIMENSION(:), INTENT(INOUT) :: tpreclist
1127     INTEGER,                       INTENT(IN)    :: nbvar
1128     INTEGER,                       INTENT(IN)    :: current_level
1129
1130     INTEGER :: ji, status
1131     CHARACTER(len=4) :: suffix
1132
1133
1134     if (infiles%files(1)%format /= NETCDF_FORMAT) return
1135
1136     write(suffix,'(I4.4)') current_level
1137
1138     DO ji=1,nbvar
1139       IF (.NOT.tpreclist(ji)%tbr) CYCLE
1140       status = NF90_INQ_VARID(infiles%files(1)%lun_id,trim(tpreclist(ji)%name)//trim(suffix),tpreclist(ji)%id_in)
1141       IF (status /= NF90_NOERR .AND. tpreclist(ji)%found) THEN
1142         tpreclist(ji)%found=.false.
1143         tpreclist(ji)%tbr=.false.
1144         tpreclist(ji)%tbw=.false.
1145         print *,'Error: variable ',trim(tpreclist(ji)%name),' not found anymore in split file'
1146       END IF
1147     END DO
1148   END SUBROUTINE UPDATE_VARID_IN
1149
1150   SUBROUTINE OPEN_FILES(infiles,outfiles,hinfile,houtfile,nbvar_infile,options,runmode)
1151     TYPE(filelist_struct),INTENT(OUT) :: infiles, outfiles
1152     CHARACTER(LEN=*), INTENT(IN)  :: hinfile
1153     CHARACTER(LEN=*), INTENT(IN)  :: houtfile
1154     INTEGER         , INTENT(OUT) :: nbvar_infile
1155     TYPE(option),DIMENSION(:),INTENT(IN) :: options
1156     INTEGER         , INTENT(IN)  :: runmode
1157
1158     INTEGER                     :: extindex
1159     INTEGER(KIND=LFI_INT)       :: iresp,iverb,inap,inaf
1160     INTEGER                     :: idx,status
1161     CHARACTER(LEN=4)            :: ypextsrc, ypextdest
1162     LOGICAL                     :: fexist
1163     INTEGER                     :: omode
1164
1165     iverb = 0
1166
1167     CALL init_sysfield()
1168
1169     IF (runmode == MODELFI2CDF) THEN
1170        ! Cas LFI -> NetCDF
1171        infiles%nbfiles = infiles%nbfiles + 1
1172        idx = infiles%nbfiles
1173        infiles%files(idx)%lun_id = 11
1174        infiles%files(idx)%format = LFI_FORMAT
1175        infiles%files(idx)%status = READING
1176        CALL LFIOUV(iresp,infiles%files(idx)%lun_id,ltrue,hinfile,'OLD',lfalse&
1177             & ,lfalse,iverb,inap,inaf)
1178        infiles%files(idx)%opened  = .TRUE.
1179
1180        nbvar_infile = inaf
1181
1182        IF (options(OPTLIST)%set) THEN
1183           CALL LFILAF(iresp,infiles%files(idx)%lun_id,lfalse)
1184           CALL LFIFER(iresp,infiles%files(idx)%lun_id,'KEEP')
1185           return
1186        END IF
1187
1188        IF (.NOT.options(OPTSPLIT)%set) THEN
1189          outfiles%nbfiles = outfiles%nbfiles + 1
1190
1191          idx = outfiles%nbfiles
1192          outfiles%files(idx)%format = NETCDF_FORMAT
1193          outfiles%files(idx)%status = WRITING
1194          IF (options(OPTCDF4)%set) THEN
1195             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1196          ELSE
1197             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1198          END IF
1199        
1200          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1201          outfiles%files(idx)%opened  = .TRUE.
1202
1203          status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1204          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1205 !!$       SELECT CASE(omode)
1206 !!$       CASE (NF90_FILL)
1207 !!$          PRINT *,'Ancien mode : NF90_FILL'
1208 !!$       CASE (NF90_NOFILL)
1209 !!$          PRINT *,'Ancien mode : NF90_NOFILL'
1210 !!$       CASE default
1211 !!$          PRINT *, 'Ancien mode : inconnu'
1212 !!$       END SELECT
1213          END IF ! .NOT.osplit
1214        
1215     ELSE IF (runmode == MODECDF2CDF) THEN
1216        ! Cas netCDF -> netCDF
1217
1218        infiles%nbfiles = infiles%nbfiles + 1
1219        idx = infiles%nbfiles
1220        status = NF90_OPEN(hinfile,NF90_NOWRITE,infiles%files(idx)%lun_id)
1221        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1222        infiles%files(idx)%opened  = .TRUE.
1223        infiles%files(idx)%format = NETCDF_FORMAT
1224        infiles%files(idx)%status = READING
1225
1226        status = NF90_INQUIRE(infiles%files(idx)%lun_id, nvariables = nbvar_infile)
1227        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1228
1229
1230        IF (.NOT.options(OPTSPLIT)%set) THEN
1231          outfiles%nbfiles = outfiles%nbfiles + 1
1232          idx = outfiles%nbfiles
1233
1234          IF (options(OPTCDF4)%set) THEN
1235             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1236          ELSE
1237             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1238          END IF
1239
1240          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1241          outfiles%files(idx)%opened  = .TRUE.
1242          outfiles%files(idx)%format = NETCDF_FORMAT
1243          outfiles%files(idx)%status = WRITING
1244
1245          status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1246          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1247        END IF ! .NOT.osplit
1248
1249     ELSE
1250        ! Cas NetCDF -> LFI
1251        infiles%nbfiles = infiles%nbfiles + 1
1252        idx = infiles%nbfiles
1253        status = NF90_OPEN(hinfile,NF90_NOWRITE,infiles%files(idx)%lun_id)
1254        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1255        infiles%files(idx)%opened  = .TRUE.
1256        infiles%files(idx)%format = NETCDF_FORMAT
1257        infiles%files(idx)%status = READING
1258        
1259        status = NF90_INQUIRE(infiles%files(idx)%lun_id, nvariables = nbvar_infile)
1260        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1261
1262        inap = 100
1263        outfiles%nbfiles = outfiles%nbfiles + 1
1264        idx = outfiles%nbfiles
1265        outfiles%files(idx)%lun_id = 11
1266        outfiles%files(idx)%format = LFI_FORMAT
1267        outfiles%files(idx)%status = WRITING
1268        CALL LFIOUV(iresp,outfiles%files(idx)%lun_id,ltrue,houtfile,'NEW'&
1269             & ,lfalse,lfalse,iverb,inap,inaf)
1270        outfiles%files(idx)%opened  = .TRUE.
1271     END IF
1272
1273     PRINT *,'--> Fichier converti : ', houtfile
1274
1275   END SUBROUTINE OPEN_FILES
1276
1277   SUBROUTINE OPEN_SPLIT_LFIFILE_IN(infiles,hinfile,current_level)
1278     TYPE(filelist_struct), INTENT(INOUT) :: infiles
1279     CHARACTER(LEN=*), INTENT(IN) :: hinfile
1280     INTEGER,          INTENT(IN) :: current_level
1281
1282     INTEGER(KIND=LFI_INT) :: ilu,iresp,iverb,inap,nbvar
1283
1284     CHARACTER(LEN=3)      :: suffix
1285     CHARACTER(LEN=:),ALLOCATABLE :: filename
1286
1287
1288     iverb = 0 !Verbosity level for LFI
1289
1290     ALLOCATE(character(len=len(hinfile)) :: filename)
1291
1292     ilu = infiles%files(1)%lun_id !We assume only 1 infile
1293
1294     write(suffix,'(I3.3)') current_level
1295     filename=hinfile(1:len(hinfile)-7)//suffix//'.lfi'
1296     CALL LFIOUV(iresp,ilu,ltrue,filename,'OLD',lfalse,lfalse,iverb,inap,nbvar)
1297     infiles%files(1)%opened = .TRUE.
1298
1299     DEALLOCATE(filename)
1300   END SUBROUTINE OPEN_SPLIT_LFIFILE_IN
1301
1302   SUBROUTINE OPEN_SPLIT_NCFILE_IN(infiles,hinfile,current_level)
1303     TYPE(filelist_struct), INTENT(INOUT) :: infiles
1304     CHARACTER(LEN=*), INTENT(IN) :: hinfile
1305     INTEGER,          INTENT(IN) :: current_level
1306
1307     INTEGER :: status
1308     CHARACTER(LEN=3)      :: suffix
1309     CHARACTER(LEN=:),ALLOCATABLE :: filename
1310
1311
1312     ALLOCATE(character(len=len(hinfile)) :: filename)
1313
1314     write(suffix,'(I3.3)') current_level
1315     filename=hinfile(1:len(hinfile)-6)//suffix//'.nc'
1316     status = NF90_OPEN(filename,NF90_NOWRITE,infiles%files(1)%lun_id)
1317     IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1318     infiles%files(1)%opened  = .TRUE.
1319
1320     DEALLOCATE(filename)
1321   END SUBROUTINE OPEN_SPLIT_NCFILE_IN
1322
1323   SUBROUTINE OPEN_SPLIT_NCFILES_OUT(outfiles,houtfile,nbvar,tpreclist,options)
1324     TYPE(filelist_struct),         INTENT(INOUT) :: outfiles
1325     CHARACTER(LEN=*),              INTENT(IN)    :: houtfile
1326     INTEGER,                       INTENT(IN)    :: nbvar
1327     TYPE(workfield), DIMENSION(:), INTENT(IN)    :: tpreclist
1328     TYPE(option),DIMENSION(:),     INTENT(IN)    :: options
1329
1330     INTEGER :: ji, idx
1331     INTEGER :: status
1332     INTEGER :: omode
1333     CHARACTER(LEN=MAXLEN) :: filename
1334
1335
1336     DO ji = 1,nbvar
1337       IF (tpreclist(ji)%tbw) outfiles%nbfiles = outfiles%nbfiles + 1
1338     END DO
1339
1340     idx = 1
1341     DO ji = 1,nbvar
1342       IF (.NOT.tpreclist(ji)%tbw) CYCLE
1343       outfiles%files(idx)%var_id = ji
1344
1345       IF (options(OPTCDF4)%set) THEN
1346         filename = trim(houtfile)//'.'//trim(tpreclist(ji)%name)//'.nc4'
1347         status = NF90_CREATE(trim(filename), IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1348       ELSE
1349         filename = trim(houtfile)//'.'//trim(tpreclist(ji)%name)//'.nc'
1350         status = NF90_CREATE(trim(filename), IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1351       END IF
1352
1353       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1354
1355       status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1356       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1357
1358       outfiles%files(idx)%opened  = .TRUE.
1359       outfiles%files(idx)%format = NETCDF_FORMAT
1360       outfiles%files(idx)%status = WRITING
1361
1362       idx = idx + 1
1363     END DO
1364
1365   END SUBROUTINE OPEN_SPLIT_NCFILES_OUT
1366   
1367   SUBROUTINE CLOSE_FILES(filelist)
1368     TYPE(filelist_struct),INTENT(INOUT) :: filelist
1369     
1370     INTEGER(KIND=LFI_INT) :: iresp
1371     INTEGER               :: ji,status
1372
1373     DO ji=1,filelist%nbfiles
1374       IF ( .NOT.filelist%files(ji)%opened ) CYCLE
1375
1376       IF ( filelist%files(ji)%format == LFI_FORMAT ) THEN
1377         CALL LFIFER(iresp,filelist%files(ji)%lun_id,'KEEP')
1378       ELSE IF ( filelist%files(ji)%format == NETCDF_FORMAT ) THEN
1379         status = NF90_CLOSE(filelist%files(ji)%lun_id)
1380         IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1381       END IF
1382       filelist%files(ji)%opened=.false.
1383     END DO
1384
1385   END SUBROUTINE CLOSE_FILES
1386
1387 END MODULE mode_util