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