31cab4c7efd5cc1e72a5e067d0dac5d981981016
[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 MODD_PARAM
5   USE netcdf
6
7   IMPLICIT NONE 
8
9   INTEGER,PARAMETER :: MAXRAW=10
10   INTEGER,PARAMETER :: MAXLEN=512
11   INTEGER,PARAMETER :: MAXFILES=100
12
13   INTEGER,PARAMETER :: UNDEFINED = -1, READING = 1, WRITING = 2
14   INTEGER,PARAMETER :: UNKNOWN_FORMAT = -1, NETCDF_FORMAT = 1, LFI_FORMAT = 2
15
16   TYPE filestruct
17     INTEGER :: lun_id                  ! Logical ID of file
18     INTEGER :: format = UNKNOWN_FORMAT ! NETCDF, LFI
19     INTEGER :: status = UNDEFINED      ! Opened for reading or writing
20     INTEGER :: var_id                  ! Position of the variable in the workfield structure
21     LOGICAL :: opened = .false.
22   END TYPE filestruct
23
24   TYPE filelist_struct
25     INTEGER :: nbfiles = 0
26 !    TYPE(filestruct),DIMENSION(:),ALLOCATABLE :: files
27     TYPE(filestruct),DIMENSION(MAXFILES) :: files
28   END TYPE filelist_struct
29
30
31   TYPE workfield
32      CHARACTER(LEN=FM_FIELD_SIZE)            :: name   ! nom du champ
33      INTEGER                                 :: TYPE   ! type (entier ou reel)    
34      CHARACTER(LEN=:), POINTER               :: comment
35      TYPE(dimCDF),                   POINTER :: dim
36      INTEGER                                 :: id_in = -1, id_out = -1
37      INTEGER                                 :: grid
38      LOGICAL                                 :: found  ! T if found in the input file
39      LOGICAL                                 :: calc   ! T if computed from other variables
40      LOGICAL                                 :: tbw    ! to be written or not
41      LOGICAL                                 :: tbr    ! to be read or not
42      INTEGER,DIMENSION(MAXRAW)               :: src    ! List of variables used to compute the variable (needed only if calc=.true.)
43      INTEGER                                 :: tgt    ! Target: id of the variable that use it (calc variable)
44   END TYPE workfield
45
46 #ifndef LOWMEM
47   TYPE lfidata
48      INTEGER(KIND=8), DIMENSION(:), POINTER :: iwtab
49   END TYPE lfidata
50   TYPE(lfidata), DIMENSION(:), ALLOCATABLE :: lfiart
51 #endif
52
53   LOGICAL(KIND=LFI_INT), PARAMETER :: ltrue  = .TRUE.
54   LOGICAL(KIND=LFI_INT), PARAMETER :: lfalse = .FALSE.
55
56 CONTAINS 
57   FUNCTION str_replace(hstr, hold, hnew)
58     CHARACTER(LEN=*) :: hstr, hold, hnew
59     CHARACTER(LEN=LEN_TRIM(hstr)+MAX(0,LEN(hnew)-LEN(hold))) :: str_replace
60     
61     INTEGER :: pos
62     
63     pos = INDEX(hstr,hold)
64     IF (pos /= 0) THEN
65        str_replace = hstr(1:pos-1)//hnew//hstr(pos+LEN(hold):)
66     ELSE 
67        str_replace = hstr 
68     END IF
69
70   END FUNCTION str_replace
71
72   SUBROUTINE FMREADLFIN1(klu,hrecfm,kval,kresp)
73   INTEGER(KIND=LFI_INT), INTENT(IN) :: klu ! logical fortran unit au lfi file
74   CHARACTER(LEN=*),INTENT(IN)       :: hrecfm ! article name to be read
75   INTEGER, INTENT(OUT)        :: kval ! integer value for hrecfm article
76   INTEGER(KIND=LFI_INT), INTENT(OUT):: kresp! return code null if OK
77   !
78   INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE::iwork
79   INTEGER :: icomlen
80   INTEGER(KIND=LFI_INT) :: iresp,ilenga,iposex
81   !
82   CALL LFINFO(iresp,klu,hrecfm,ilenga,iposex)
83   IF (iresp /=0 .OR. ilenga == 0) THEN
84     kresp = -1
85     kval = 0
86   ELSE
87     ALLOCATE(IWORK(ilenga))
88     CALL LFILEC(iresp,klu,hrecfm,iwork,ilenga)
89     icomlen = iwork(2)
90     kval = iwork(3+icomlen)
91     kresp = iresp
92     DEALLOCATE(IWORK)
93   END IF
94   END SUBROUTINE FMREADLFIN1
95
96   SUBROUTINE parse_infiles(infiles, hvarlist, nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw, tpreclist, kbuflen, icurrent_level)
97     TYPE(filelist_struct),      INTENT(IN) :: infiles
98     INTEGER,                    INTENT(IN) :: nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw
99     CHARACTER(LEN=*),           INTENT(IN) :: hvarlist
100     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
101     INTEGER,                   INTENT(OUT) :: kbuflen
102     INTEGER,          INTENT(IN), OPTIONAL :: icurrent_level
103
104     INTEGER                                  :: ji,jj, kcdf_id, itype
105     INTEGER                                  :: ndb, nde, ndey, idx, idx_var, maxvar
106     INTEGER                                  :: idims, idimtmp, jdim, status, var_id
107     LOGICAL                                  :: ladvan
108     INTEGER                                  :: ich, current_level, leng
109     INTEGER                                  :: comment_size, fsize, sizemax
110     CHARACTER(LEN=FM_FIELD_SIZE)             :: yrecfm
111     CHARACTER(LEN=4)                         :: suffix
112 #ifdef LOWMEM
113     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
114 #endif
115     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
116     CHARACTER(LEN=FM_FIELD_SIZE)             :: var_calc
117     CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
118     INTEGER, DIMENSION(10)                   :: idim_id
119     !JUAN CYCCL3
120     INTEGER                        :: JPHEXT
121
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 (LEN_TRIM(hvarlist) > 0) 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           nde = INDEX(TRIM(hvarlist(ndb:)),',')
201           yrecfm = hvarlist(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, hvarlist, nbvar, tpreclist, kbuflen, current_level)
474     TYPE(filelist_struct),      INTENT(IN) :: infiles
475     INTEGER, INTENT(INOUT)                 :: nbvar
476     CHARACTER(LEN=*), intent(IN)           :: hvarlist
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,oreduceprecision,omerge,osplit,ocompress,compress_level)
536     TYPE(filelist_struct),       INTENT(IN) :: outfiles
537     TYPE(workfield),DIMENSION(:),INTENT(INOUT) :: tpreclist
538     INTEGER,                     INTENT(IN) :: nbvar
539     LOGICAL,                     INTENT(IN) :: oreduceprecision
540     LOGICAL,                     INTENT(IN) :: omerge
541     LOGICAl,                     INTENT(IN) :: osplit
542     LOGICAL,                     INTENT(IN) :: ocompress
543     INTEGER,                     INTENT(IN) :: compress_level
544
545     INTEGER :: status
546     INTEGER :: idx, ji, nbfiles
547     INTEGER:: kcdf_id
548     TYPE(dimCDF), POINTER :: tzdim
549     INTEGER               :: invdims
550     INTEGER               :: type_float
551     INTEGER, DIMENSION(10) :: ivdims
552     CHARACTER(LEN=20)     :: ycdfvar
553
554
555     nbfiles = outfiles%nbfiles
556
557     IF (oreduceprecision) THEN
558       type_float = NF90_REAL
559     ELSE
560       type_float = NF90_DOUBLE
561     END IF
562
563     DO ji = 1,nbfiles
564       kcdf_id = outfiles%files(ji)%lun_id
565
566       ! global attributes
567       status = NF90_PUT_ATT(kcdf_id,NF90_GLOBAL,'Title',VERSION_ID)
568       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
569
570       ! define DIMENSIONS
571       tzdim=>first_DimCDF()
572       DO WHILE(ASSOCIATED(tzdim))
573         IF (tzdim%create) THEN
574           status = NF90_DEF_DIM(kcdf_id,tzdim%name,tzdim%len,tzdim%id)
575           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
576         END IF
577         tzdim=>tzdim%next
578       END DO
579     END DO
580
581     PRINT *,'------------- NetCDF DEFINITION ---------------'
582
583     ! define VARIABLES and ATTRIBUTES
584     idx = 1
585     DO ji=1,nbvar
586        IF (.NOT.tpreclist(ji)%tbw) CYCLE
587
588        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
589          IF (tpreclist(ji)%dim%create) THEN
590            invdims   = 1
591            ivdims(1) = tpreclist(ji)%dim%id
592          ELSE
593            invdims = tpreclist(ji)%dim%ndims
594            IF(omerge) invdims=invdims+1 !when merging variables from LFI splitted files
595            SELECT CASE(invdims)
596            CASE(2)
597               ivdims(1)=ptdimx%id
598               ivdims(2)=ptdimy%id
599            CASE(3)
600               ivdims(1)=ptdimx%id
601               ivdims(2)=ptdimy%id
602               ivdims(3)=ptdimz%id
603            CASE(12)
604               ivdims(1)=ptdimx%id
605               ivdims(2)=ptdimz%id
606               invdims = 2 ! on retablit la bonne valeur du nbre de dimension
607            CASE default
608              PRINT *,'Fatal error in NetCDF dimension definition'
609              STOP
610            END SELECT
611          END IF
612        ELSE
613          ! scalar variables
614           invdims   = 0
615           ivdims(1) = 0 ! ignore dans ce cas
616        END IF
617        
618        ! Variables definition
619
620        !! NetCDF n'aime pas les '%' dans le nom des variables
621        !! "%" remplaces par '__' 
622        ycdfvar = str_replace(tpreclist(ji)%name,'%','__')
623        !! ni les '.' remplaces par '--'
624        ycdfvar = str_replace(ycdfvar,'.','--')
625
626        kcdf_id = outfiles%files(idx)%lun_id
627
628        SELECT CASE(tpreclist(ji)%TYPE)
629        CASE (TEXT)
630 !          PRINT *,'TEXT : ',tpreclist(ji)%name
631           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_CHAR,&
632                    ivdims(:invdims),tpreclist(ji)%id_out)
633           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
634
635        CASE (INT,BOOL)
636 !          PRINT *,'INT,BOOL : ',tpreclist(ji)%name
637           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_INT,&
638                    ivdims(:invdims),tpreclist(ji)%id_out)
639           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
640
641        CASE(FLOAT)
642 !          PRINT *,'FLOAT : ',tpreclist(ji)%name
643           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
644                    ivdims(:invdims),tpreclist(ji)%id_out)
645           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
646
647           
648        CASE default
649           PRINT *,'ATTENTION : ',TRIM(tpreclist(ji)%name),' est de&
650                & TYPE inconnu --> force a REAL'
651           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
652                    ivdims(:invdims),tpreclist(ji)%id_out)
653           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
654           
655
656        END SELECT
657
658        ! Compress data (costly operation for the CPU)
659        IF (ocompress .AND. invdims>0) THEN
660          status = NF90_DEF_VAR_DEFLATE(kcdf_id,tpreclist(ji)%id_out,1,1,compress_level)
661          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
662        END IF
663
664        ! GRID attribute definition
665        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'GRID',tpreclist(ji)%grid)
666        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
667
668        ! COMMENT attribute definition
669        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'COMMENT',trim(tpreclist(ji)%comment))
670        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
671
672        IF (osplit) idx = idx + 1
673     END DO
674     
675     DO ji = 1,nbfiles
676       kcdf_id = outfiles%files(ji)%lun_id
677       status = NF90_ENDDEF(kcdf_id)
678       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
679     END DO
680
681   END SUBROUTINE def_ncdf
682
683   SUBROUTINE fill_ncdf(infiles,outfiles,tpreclist,knaf,kbuflen,osplit,current_level)
684     TYPE(filelist_struct),        INTENT(IN):: infiles, outfiles
685     TYPE(workfield), DIMENSION(:),INTENT(IN):: tpreclist
686     INTEGER,                      INTENT(IN):: knaf
687     INTEGER,                      INTENT(IN):: kbuflen
688     LOGICAl,                      INTENT(IN):: osplit
689     INTEGER, INTENT(IN), OPTIONAL           :: current_level
690
691 #ifdef LOWMEM
692     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
693 #endif
694     INTEGER                                  :: idx, ji,jj
695     INTEGER                                  :: kcdf_id
696     INTEGER                                  :: status
697     INTEGER                                  :: extent, ndims
698     INTEGER                                  :: ich
699     INTEGER                                  :: src
700     INTEGER                                  :: level
701     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
702     CHARACTER(LEN=4)                         :: suffix
703     INTEGER,DIMENSION(3)                     :: idims, start
704     INTEGER,DIMENSION(:),ALLOCATABLE         :: itab
705     REAL(KIND=8),DIMENSION(:),ALLOCATABLE    :: xtab
706     CHARACTER, DIMENSION(:), ALLOCATABLE     :: ytab
707     REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: xtab3d, xtab3d2
708     INTEGER,      DIMENSION(:,:,:), ALLOCATABLE :: itab3d, itab3d2
709
710
711     !
712     IF (infiles%files(1)%format == LFI_FORMAT) ilu = infiles%files(1)%lun_id
713     !
714
715     IF (present(current_level)) THEN
716       write(suffix,'(I4.4)') current_level
717       level = current_level
718     ElSE
719       suffix=''
720       level = 1
721     END IF
722
723 #if LOWMEM
724     ALLOCATE(iwork(kbuflen))
725 #endif
726     ALLOCATE(itab(kbuflen))
727     ALLOCATE(xtab(kbuflen))
728
729     idx = 1
730     DO ji=1,knaf
731        IF (.NOT.tpreclist(ji)%tbw) CYCLE
732
733        kcdf_id = outfiles%files(idx)%lun_id
734
735        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
736           extent = tpreclist(ji)%dim%len
737           ndims = tpreclist(ji)%dim%ndims
738        ELSE
739           extent = 1
740           ndims = 0
741        END IF
742
743        idims(:) = 1
744        if(ndims>0) idims(1) = ptdimx%len
745        if(ndims>1) idims(2) = ptdimy%len
746        if(ndims>2) idims(3) = ptdimz%len
747        if(ndims>3) then
748          PRINT *,'Too many dimensions'
749          STOP
750        endif
751
752        SELECT CASE(tpreclist(ji)%TYPE)
753        CASE (INT,BOOL)
754         IF (infiles%files(1)%format == LFI_FORMAT) THEN
755 #if LOWMEM
756          IF (.NOT.tpreclist(ji)%calc) THEN
757            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
758            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
759            itab(1:extent) = iwork(3+iwork(2):)
760          ELSE
761            src=tpreclist(ji)%src(1)
762            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
763            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
764            itab(1:extent) = iwork(3+iwork(2):)
765            jj = 2
766            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
767              src=tpreclist(ji)%src(jj)
768              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
769              itab(1:extent) = itab(1:extent) + iwork(3+iwork(2):)
770              jj=jj+1
771            END DO
772          ENDIF
773 #else
774          IF (.NOT.tpreclist(ji)%calc) THEN
775            itab(1:extent) = lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):)
776          ELSE
777            src=tpreclist(ji)%src(1)
778            itab(1:extent) = lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
779            jj = 2
780            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
781              src=tpreclist(ji)%src(jj)
782              itab(1:extent) = xtab(1:extent) + lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
783              jj=jj+1
784            END DO
785          END IF
786 #endif
787
788 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
789          SELECT CASE(ndims)
790          CASE (0)
791            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1))
792          CASE (1)
793            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1:extent),count=(/extent/))
794          CASE (2)
795            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len/)), &
796                                  start = (/1,1,level/) )
797          CASE (3)
798            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
799          CASE DEFAULT
800            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
801          END SELECT
802
803         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
804          ALLOCATE( itab3d(idims(1),idims(2),idims(3)) )
805          IF (.NOT.tpreclist(ji)%calc) THEN
806            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,itab3d)
807            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
808          ELSE
809            ALLOCATE( itab3d2(idims(1),idims(2),idims(3)) )
810            src=tpreclist(ji)%src(1)
811            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d)
812            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
813            jj = 2
814            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
815              src=tpreclist(ji)%src(jj)
816              status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d2)
817              IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
818              itab3d(:,:,:) = itab3d(:,:,:) + itab3d2(:,:,:)
819              jj=jj+1
820            END DO
821            DEALLOCATE(itab3d2)
822          END IF
823
824 !TODO: not clean, should be done only if merging z-levels
825          IF (ndims == 2) THEN
826            start = (/1,1,level/)
827          ELSE
828            start = (/1,1,1/)
829          ENDIF
830          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab3d,start=start)
831          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
832
833          DEALLOCATE(itab3d)
834         END IF
835
836          
837        CASE (FLOAT)
838         IF (infiles%files(1)%format == LFI_FORMAT) THEN
839 #if LOWMEM
840          IF (.NOT.tpreclist(ji)%calc) THEN
841            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
842            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
843            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
844          ELSE
845            src=tpreclist(ji)%src(1)
846            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
847            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
848            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
849            jj = 2
850            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
851              src=tpreclist(ji)%src(jj)
852              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
853              xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
854              jj=jj+1
855            END DO
856          ENDIF
857 #else
858          IF (.NOT.tpreclist(ji)%calc) THEN
859            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
860          ELSE
861            src=tpreclist(ji)%src(1)
862            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
863            jj = 2
864            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
865              src=tpreclist(ji)%src(jj)
866              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
867              jj=jj+1
868            END DO
869          END IF
870 #endif
871 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
872          SELECT CASE(ndims)
873          CASE (0)
874            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
875          CASE (1)
876            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
877          CASE (2)
878            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
879                                  start = (/1,1,level/) )
880          CASE (3)
881            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
882          CASE DEFAULT
883            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
884          END SELECT
885
886         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
887          ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
888          IF (.NOT.tpreclist(ji)%calc) THEN
889            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,xtab3d)
890            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
891          ELSE
892            ALLOCATE( xtab3d2(idims(1),idims(2),idims(3)) )
893            src=tpreclist(ji)%src(1)
894            status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d)
895            IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
896            jj = 2
897            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
898              src=tpreclist(ji)%src(jj)
899              status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d2)
900              IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
901              xtab3d(:,:,:) = xtab3d(:,:,:) + xtab3d2(:,:,:)
902              jj=jj+1
903            END DO
904            DEALLOCATE(xtab3d2)
905          END IF
906
907 !TODO: not clean, should be done only if merging z-levels
908          IF (ndims == 2) THEN
909            start = (/1,1,level/)
910          ELSE
911            start = (/1,1,1/)
912          ENDIF
913          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab3d,start=start)
914          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
915
916          DEALLOCATE(xtab3d)
917         END IF
918
919        CASE (TEXT)
920         IF (infiles%files(1)%format == LFI_FORMAT) THEN
921          ALLOCATE(ytab(extent))
922          DO jj=1,extent
923 #if LOWMEM
924            ich = iwork(2+iwork(2)+jj)
925 #else
926            ich = lfiart(ji)%iwtab(2+lfiart(ji)%iwtab(2)+jj)
927 #endif
928            ytab(jj) = CHAR(ich)
929          END DO
930          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
931          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
932          DEALLOCATE(ytab)
933         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
934          status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,ytab,count=(/extent/))
935          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
936
937          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
938          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
939         END IF
940
941        CASE default
942         IF (infiles%files(1)%format == LFI_FORMAT) THEN
943 #if LOWMEM
944          IF (.NOT.tpreclist(ji)%calc) THEN
945            CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
946            CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
947            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
948          ELSE
949            src=tpreclist(ji)%src(1)
950            CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
951            CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
952            xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
953            jj = 2
954            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
955              src=tpreclist(ji)%src(jj)
956              CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
957              xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
958              jj=jj+1
959            END DO
960          ENDIF
961 #else         
962          IF (.NOT.tpreclist(ji)%calc) THEN
963            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
964          ELSE
965            src=tpreclist(ji)%src(1)
966            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
967            jj = 2
968            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
969              src=tpreclist(ji)%src(jj)
970              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
971              jj=jj+1
972            END DO
973          END IF
974 #endif
975 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
976          SELECT CASE(ndims)
977          CASE (0)
978            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
979          CASE (1)
980            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
981          CASE (2)
982            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
983                                  start = (/1,1,level/) )
984          CASE (3)
985            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
986          CASE DEFAULT
987            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
988          END SELECT
989          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
990         ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
991          print *,'Error: unknown datatype'
992          STOP
993         END IF
994
995        END SELECT
996
997        if (osplit) idx = idx + 1
998     END DO
999     DEALLOCATE(itab,xtab)
1000 #if LOWMEM
1001     DEALLOCATE(iwork)
1002 #endif 
1003   END SUBROUTINE fill_ncdf
1004
1005   SUBROUTINE build_lfi(infiles,outfiles,tpreclist,kbuflen)
1006     TYPE(filelist_struct),         INTENT(IN) :: infiles, outfiles
1007     TYPE(workfield), DIMENSION(:), INTENT(IN) :: tpreclist
1008     INTEGER,                       INTENT(IN) :: kbuflen
1009     
1010     INTEGER :: kcdf_id, status
1011     INTEGER :: ivar,ji,jj,ndims
1012     INTEGER,DIMENSION(3) :: idims
1013     INTEGER(KIND=8), DIMENSION(:), POINTER  :: iwork
1014     INTEGER(KIND=8), DIMENSION(:), POINTER  :: idata
1015     REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: xtab3d
1016     INTEGER,      DIMENSION(:,:,:), ALLOCATABLE :: itab3d
1017     CHARACTER,    DIMENSION(:), ALLOCATABLE :: ytab
1018     CHARACTER(LEN=FM_FIELD_SIZE)            :: yrecfm
1019
1020     INTEGER :: iartlen, idlen, icomlen
1021     INTEGER(KIND=LFI_INT) :: iresp,ilu,iartlen8
1022
1023
1024     ilu = outfiles%files(1)%lun_id
1025     kcdf_id = infiles%files(1)%lun_id
1026
1027     ! Un article LFI est compose de :
1028     !   - 1 entier identifiant le numero de grille
1029     !   - 1 entier contenant la taille du commentaire
1030     !   - le commentaire code en entier 64 bits
1031     !   - les donnees proprement dites
1032
1033     PRINT *,'Taille buffer = ',2+kbuflen
1034
1035     ALLOCATE(iwork(2+kbuflen))
1036
1037     DO ivar=1,SIZE(tpreclist)
1038        icomlen = LEN(tpreclist(ivar)%comment)
1039
1040        ! traitement Grille et Commentaire
1041        iwork(1) = tpreclist(ivar)%grid
1042        iwork(2) = icomlen
1043        DO jj=1,iwork(2)
1044           iwork(2+jj)=ICHAR(tpreclist(ivar)%comment(jj:jj))
1045        END DO
1046
1047        IF (ASSOCIATED(tpreclist(ivar)%dim)) THEN
1048           idlen = tpreclist(ivar)%dim%len
1049           ndims = tpreclist(ivar)%dim%ndims
1050        ELSE 
1051           idlen = 1
1052           ndims = 0
1053        END IF
1054        
1055        idims(:) = 1
1056        if(ndims>0) idims(1) = ptdimx%len
1057        if(ndims>1) idims(2) = ptdimy%len
1058        if(ndims>2) idims(3) = ptdimz%len
1059        if(ndims>3) then
1060          PRINT *,'Too many dimensions'
1061          STOP
1062        endif
1063
1064        iartlen = 2+icomlen+idlen
1065        idata=>iwork(3+icomlen:iartlen)
1066
1067
1068        SELECT CASE(tpreclist(ivar)%TYPE)
1069        CASE(INT,BOOL)
1070           ALLOCATE( itab3d(idims(1),idims(2),idims(3)) )
1071           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,itab3d)
1072           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1073
1074 !          PRINT *,'INT,BOOL --> ',tpreclist(ivar)%name,',len = ',idlen
1075           idata(1:idlen) = RESHAPE( itab3d , (/ idims(1)*idims(2)*idims(3) /) )
1076
1077           DEALLOCATE(itab3d)
1078
1079        CASE(FLOAT)
1080           ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
1081           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,xtab3d)
1082           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1083
1084 !          PRINT *,'FLOAT -->    ',tpreclist(ivar)%name,',len = ',idlen
1085           idata(1:idlen) = RESHAPE( TRANSFER(xtab3d,(/ 0_8 /),idlen) , (/ idims(1)*idims(2)*idims(3) /) )
1086
1087           DEALLOCATE(xtab3d)
1088
1089        CASE(TEXT)
1090           ALLOCATE(ytab(idlen))
1091           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,ytab)
1092           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1093
1094 !          PRINT *,'TEXT -->     ',tpreclist(ivar)%name,',len = ',idlen
1095           DO jj=1,idlen
1096              idata(jj) = ICHAR(ytab(jj))
1097           END DO
1098           
1099           DEALLOCATE(ytab)
1100
1101        CASE default
1102           ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
1103           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id_in,xtab3d)
1104           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1105
1106           PRINT *,'Default (ERROR) -->',tpreclist(ivar)%name,',len = ',idlen
1107           idata(1:idlen) = RESHAPE( TRANSFER(xtab3d,(/ 0_8 /),idlen) , (/ idims(1)*idims(2)*idims(3) /) )
1108
1109           DEALLOCATE(xtab3d)
1110
1111        END SELECT
1112        
1113        ! Attention restoration des '%' dans le nom des champs LFI
1114        yrecfm = str_replace(tpreclist(ivar)%name,'__','%')
1115        ! et des '.'
1116        yrecfm = str_replace(yrecfm,'--','.')
1117        iartlen8 = iartlen
1118        CALL LFIECR(iresp,ilu,yrecfm,iwork,iartlen8)
1119
1120     END DO
1121     DEALLOCATE(iwork)
1122
1123   END SUBROUTINE build_lfi
1124
1125   SUBROUTINE UPDATE_VARID_IN(infiles,hinfile,tpreclist,nbvar,current_level)
1126     !Update the id_in for netCDF files (could change from one file to the other)
1127     TYPE(filelist_struct),         INTENT(IN)    :: infiles
1128     CHARACTER(LEN=*),              INTENT(IN)    :: hinfile
1129     TYPE(workfield), DIMENSION(:), INTENT(INOUT) :: tpreclist
1130     INTEGER,                       INTENT(IN)    :: nbvar
1131     INTEGER,                       INTENT(IN)    :: current_level
1132
1133     INTEGER :: ji, status
1134     CHARACTER(len=4) :: suffix
1135
1136
1137     if (infiles%files(1)%format /= NETCDF_FORMAT) return
1138
1139     write(suffix,'(I4.4)') current_level
1140
1141     DO ji=1,nbvar
1142       IF (.NOT.tpreclist(ji)%tbr) CYCLE
1143       status = NF90_INQ_VARID(infiles%files(1)%lun_id,trim(tpreclist(ji)%name)//trim(suffix),tpreclist(ji)%id_in)
1144       IF (status /= NF90_NOERR .AND. tpreclist(ji)%found) THEN
1145         tpreclist(ji)%found=.false.
1146         tpreclist(ji)%tbr=.false.
1147         tpreclist(ji)%tbw=.false.
1148         print *,'Error: variable ',trim(tpreclist(ji)%name),' not found anymore in split file'
1149       END IF
1150     END DO
1151   END SUBROUTINE UPDATE_VARID_IN
1152
1153   SUBROUTINE OPEN_FILES(infiles,outfiles,hinfile,houtfile,ocdf2cdf,olfi2cdf,olfilist,ohdf5,nbvar_infile,osplit)
1154     TYPE(filelist_struct),INTENT(OUT) :: infiles, outfiles
1155     LOGICAL,          INTENT(IN)  :: ocdf2cdf, olfi2cdf, olfilist, ohdf5, osplit
1156     CHARACTER(LEN=*), INTENT(IN)  :: hinfile
1157     CHARACTER(LEN=*), INTENT(IN)  :: houtfile
1158     INTEGER         , INTENT(OUT) :: nbvar_infile
1159
1160     INTEGER                     :: extindex
1161     INTEGER(KIND=LFI_INT)       :: iresp,iverb,inap,inaf
1162     INTEGER                     :: idx,status
1163     CHARACTER(LEN=4)            :: ypextsrc, ypextdest
1164     LOGICAL                     :: fexist
1165     INTEGER                     :: omode
1166
1167     iverb = 0
1168
1169     CALL init_sysfield()
1170
1171     IF (olfi2cdf) THEN 
1172        ! Cas LFI -> NetCDF
1173        infiles%nbfiles = infiles%nbfiles + 1
1174        idx = infiles%nbfiles
1175        infiles%files(idx)%lun_id = 11
1176        infiles%files(idx)%format = LFI_FORMAT
1177        infiles%files(idx)%status = READING
1178        CALL LFIOUV(iresp,infiles%files(idx)%lun_id,ltrue,hinfile,'OLD',lfalse&
1179             & ,lfalse,iverb,inap,inaf)
1180        infiles%files(idx)%opened  = .TRUE.
1181
1182        nbvar_infile = inaf
1183
1184        IF (olfilist) THEN
1185           CALL LFILAF(iresp,infiles%files(idx)%lun_id,lfalse)
1186           CALL LFIFER(iresp,infiles%files(idx)%lun_id,'KEEP')
1187           return
1188        END IF
1189
1190        IF (.NOT.osplit) THEN
1191          outfiles%nbfiles = outfiles%nbfiles + 1
1192
1193          idx = outfiles%nbfiles
1194          outfiles%files(idx)%format = NETCDF_FORMAT
1195          outfiles%files(idx)%status = WRITING
1196          IF (ohdf5) THEN
1197             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1198          ELSE
1199             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1200          END IF
1201        
1202          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1203          outfiles%files(idx)%opened  = .TRUE.
1204
1205          status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1206          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1207 !!$       SELECT CASE(omode)
1208 !!$       CASE (NF90_FILL)
1209 !!$          PRINT *,'Ancien mode : NF90_FILL'
1210 !!$       CASE (NF90_NOFILL)
1211 !!$          PRINT *,'Ancien mode : NF90_NOFILL'
1212 !!$       CASE default
1213 !!$          PRINT *, 'Ancien mode : inconnu'
1214 !!$       END SELECT
1215          END IF ! .NOT.osplit
1216        
1217     ELSE IF (ocdf2cdf) THEN
1218        ! Cas netCDF -> netCDF
1219
1220        infiles%nbfiles = infiles%nbfiles + 1
1221        idx = infiles%nbfiles
1222        status = NF90_OPEN(hinfile,NF90_NOWRITE,infiles%files(idx)%lun_id)
1223        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1224        infiles%files(idx)%opened  = .TRUE.
1225        infiles%files(idx)%format = NETCDF_FORMAT
1226        infiles%files(idx)%status = READING
1227
1228        status = NF90_INQUIRE(infiles%files(idx)%lun_id, nvariables = nbvar_infile)
1229        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1230
1231
1232        IF (.NOT.osplit) THEN
1233          outfiles%nbfiles = outfiles%nbfiles + 1
1234          idx = outfiles%nbfiles
1235
1236          IF (ohdf5) THEN
1237             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1238          ELSE
1239             status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1240          END IF
1241
1242          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1243          outfiles%files(idx)%opened  = .TRUE.
1244          outfiles%files(idx)%format = NETCDF_FORMAT
1245          outfiles%files(idx)%status = WRITING
1246
1247          status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1248          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1249        END IF ! .NOT.osplit
1250
1251     ELSE
1252        ! Cas NetCDF -> LFI
1253        infiles%nbfiles = infiles%nbfiles + 1
1254        idx = infiles%nbfiles
1255        status = NF90_OPEN(hinfile,NF90_NOWRITE,infiles%files(idx)%lun_id)
1256        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1257        infiles%files(idx)%opened  = .TRUE.
1258        infiles%files(idx)%format = NETCDF_FORMAT
1259        infiles%files(idx)%status = READING
1260        
1261        status = NF90_INQUIRE(infiles%files(idx)%lun_id, nvariables = nbvar_infile)
1262        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1263
1264        inap = 100
1265        outfiles%nbfiles = outfiles%nbfiles + 1
1266        idx = outfiles%nbfiles
1267        outfiles%files(idx)%lun_id = 11
1268        outfiles%files(idx)%format = LFI_FORMAT
1269        outfiles%files(idx)%status = WRITING
1270        CALL LFIOUV(iresp,outfiles%files(idx)%lun_id,ltrue,houtfile,'NEW'&
1271             & ,lfalse,lfalse,iverb,inap,inaf)
1272        outfiles%files(idx)%opened  = .TRUE.
1273     END IF
1274
1275     PRINT *,'--> Fichier converti : ', houtfile
1276
1277   END SUBROUTINE OPEN_FILES
1278
1279   SUBROUTINE OPEN_SPLIT_LFIFILE_IN(infiles,hinfile,current_level)
1280     TYPE(filelist_struct), INTENT(INOUT) :: infiles
1281     CHARACTER(LEN=*), INTENT(IN) :: hinfile
1282     INTEGER,          INTENT(IN) :: current_level
1283
1284     INTEGER(KIND=LFI_INT) :: ilu,iresp,iverb,inap,nbvar
1285
1286     CHARACTER(LEN=3)      :: suffix
1287     CHARACTER(LEN=:),ALLOCATABLE :: filename
1288
1289
1290     iverb = 0 !Verbosity level for LFI
1291
1292     ALLOCATE(character(len=len(hinfile)) :: filename)
1293
1294     ilu = infiles%files(1)%lun_id !We assume only 1 infile
1295
1296     write(suffix,'(I3.3)') current_level
1297     filename=hinfile(1:len(hinfile)-7)//suffix//'.lfi'
1298     CALL LFIOUV(iresp,ilu,ltrue,filename,'OLD',lfalse,lfalse,iverb,inap,nbvar)
1299     infiles%files(1)%opened = .TRUE.
1300
1301     DEALLOCATE(filename)
1302   END SUBROUTINE OPEN_SPLIT_LFIFILE_IN
1303
1304   SUBROUTINE OPEN_SPLIT_NCFILE_IN(infiles,hinfile,current_level)
1305     TYPE(filelist_struct), INTENT(INOUT) :: infiles
1306     CHARACTER(LEN=*), INTENT(IN) :: hinfile
1307     INTEGER,          INTENT(IN) :: current_level
1308
1309     INTEGER :: status
1310     CHARACTER(LEN=3)      :: suffix
1311     CHARACTER(LEN=:),ALLOCATABLE :: filename
1312
1313
1314     ALLOCATE(character(len=len(hinfile)) :: filename)
1315
1316     write(suffix,'(I3.3)') current_level
1317     filename=hinfile(1:len(hinfile)-6)//suffix//'.nc'
1318     status = NF90_OPEN(filename,NF90_NOWRITE,infiles%files(1)%lun_id)
1319     IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1320     infiles%files(1)%opened  = .TRUE.
1321
1322     DEALLOCATE(filename)
1323   END SUBROUTINE OPEN_SPLIT_NCFILE_IN
1324
1325   SUBROUTINE OPEN_SPLIT_NCFILES_OUT(outfiles,houtfile,nbvar,tpreclist,ohdf5)
1326     TYPE(filelist_struct),         INTENT(INOUT) :: outfiles
1327     CHARACTER(LEN=*),              INTENT(IN)    :: houtfile
1328     INTEGER,                       INTENT(IN)    :: nbvar
1329     TYPE(workfield), DIMENSION(:), INTENT(IN)    :: tpreclist
1330     LOGICAL,                       INTENT(IN)    :: ohdf5
1331
1332     INTEGER :: ji, idx
1333     INTEGER :: status
1334     INTEGER :: omode
1335     CHARACTER(LEN=MAXLEN) :: filename
1336
1337
1338     DO ji = 1,nbvar
1339       IF (tpreclist(ji)%tbw) outfiles%nbfiles = outfiles%nbfiles + 1
1340     END DO
1341
1342     idx = 1
1343     DO ji = 1,nbvar
1344       IF (.NOT.tpreclist(ji)%tbw) CYCLE
1345       outfiles%files(idx)%var_id = ji
1346
1347       IF (ohdf5) THEN
1348         filename = trim(houtfile)//'.'//trim(tpreclist(ji)%name)//'.nc4'
1349         status = NF90_CREATE(trim(filename), IOR(NF90_CLOBBER,NF90_NETCDF4), outfiles%files(idx)%lun_id)
1350       ELSE
1351         filename = trim(houtfile)//'.'//trim(tpreclist(ji)%name)//'.nc'
1352         status = NF90_CREATE(trim(filename), IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), outfiles%files(idx)%lun_id)
1353       END IF
1354
1355       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1356
1357       status = NF90_SET_FILL(outfiles%files(idx)%lun_id,NF90_NOFILL,omode)
1358       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1359
1360       outfiles%files(idx)%opened  = .TRUE.
1361       outfiles%files(idx)%format = NETCDF_FORMAT
1362       outfiles%files(idx)%status = WRITING
1363
1364       idx = idx + 1
1365     END DO
1366
1367   END SUBROUTINE OPEN_SPLIT_NCFILES_OUT
1368   
1369   SUBROUTINE CLOSE_FILES(filelist)
1370     TYPE(filelist_struct),INTENT(INOUT) :: filelist
1371     
1372     INTEGER(KIND=LFI_INT) :: iresp
1373     INTEGER               :: ji,status
1374
1375     DO ji=1,filelist%nbfiles
1376       IF ( .NOT.filelist%files(ji)%opened ) CYCLE
1377
1378       IF ( filelist%files(ji)%format == LFI_FORMAT ) THEN
1379         CALL LFIFER(iresp,filelist%files(ji)%lun_id,'KEEP')
1380       ELSE IF ( filelist%files(ji)%format == NETCDF_FORMAT ) THEN
1381         status = NF90_CLOSE(filelist%files(ji)%lun_id)
1382         IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
1383       END IF
1384       filelist%files(ji)%opened=.false.
1385     END DO
1386
1387   END SUBROUTINE CLOSE_FILES
1388
1389 END MODULE mode_util