lfi2cdf: is now able to create new variables from read variables.
[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
11   TYPE workfield
12      CHARACTER(LEN=FM_FIELD_SIZE)            :: name   ! nom du champ
13      INTEGER                                 :: TYPE   ! type (entier ou reel)    
14      CHARACTER(LEN=:), POINTER               :: comment
15      TYPE(dimCDF),                   POINTER :: dim
16      INTEGER                                 :: id
17      INTEGER                                 :: grid
18      LOGICAL                                 :: found  ! T if found in the input file
19      LOGICAL                                 :: calc   ! T if computed from other variables
20      LOGICAL                                 :: tbw    ! to be written or not
21      LOGICAL                                 :: tbr    ! to be read or not
22      INTEGER,DIMENSION(MAXRAW)               :: src    ! List of variables used to compute the variable (needed only if calc=.true.)
23      INTEGER                                 :: tgt    ! Target: id of the variable that use it (calc variable)
24   END TYPE workfield
25
26 #ifndef LOWMEM
27   TYPE lfidata
28      INTEGER(KIND=8), DIMENSION(:), POINTER :: iwtab
29   END TYPE lfidata
30   TYPE(lfidata), DIMENSION(:), ALLOCATABLE :: lfiart
31 #endif
32
33   LOGICAL(KIND=LFI_INT), PARAMETER :: ltrue  = .TRUE.
34   LOGICAL(KIND=LFI_INT), PARAMETER :: lfalse = .FALSE.
35
36 CONTAINS 
37   FUNCTION str_replace(hstr, hold, hnew)
38     CHARACTER(LEN=*) :: hstr, hold, hnew
39     CHARACTER(LEN=LEN_TRIM(hstr)+MAX(0,LEN(hnew)-LEN(hold))) :: str_replace
40     
41     INTEGER :: pos
42     
43     pos = INDEX(hstr,hold)
44     IF (pos /= 0) THEN
45        str_replace = hstr(1:pos-1)//hnew//hstr(pos+LEN(hold):)
46     ELSE 
47        str_replace = hstr 
48     END IF
49
50   END FUNCTION str_replace
51
52   SUBROUTINE FMREADLFIN1(klu,hrecfm,kval,kresp)
53   INTEGER(KIND=LFI_INT), INTENT(IN) :: klu ! logical fortran unit au lfi file
54   CHARACTER(LEN=*),INTENT(IN)       :: hrecfm ! article name to be read
55   INTEGER, INTENT(OUT)        :: kval ! integer value for hrecfm article
56   INTEGER(KIND=LFI_INT), INTENT(OUT):: kresp! return code null if OK
57   !
58   INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE::iwork
59   INTEGER :: icomlen
60   INTEGER(KIND=LFI_INT) :: iresp,ilenga,iposex
61   !
62   CALL LFINFO(iresp,klu,hrecfm,ilenga,iposex)
63   IF (iresp /=0 .OR. ilenga == 0) THEN
64     kresp = -1
65     kval = 0
66   ELSE
67     ALLOCATE(IWORK(ilenga))
68     CALL LFILEC(iresp,klu,hrecfm,iwork,ilenga)
69     icomlen = iwork(2)
70     kval = iwork(3+icomlen)
71     kresp = iresp
72     DEALLOCATE(IWORK)
73   END IF
74   END SUBROUTINE FMREADLFIN1
75
76   SUBROUTINE parse_lfi(klu, hvarlist, nbvar_lfi, nbvar_tbr, nbvar_calc, nbvar_tbw, tpreclist, kbuflen, current_level)
77     INTEGER, INTENT(IN)                    :: klu
78     INTEGER, INTENT(IN)                    :: nbvar_lfi, nbvar_tbr, nbvar_calc, nbvar_tbw
79     CHARACTER(LEN=*), intent(IN)           :: hvarlist
80     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist    
81     INTEGER, INTENT(OUT)                   :: kbuflen
82     INTEGER, INTENT(IN), OPTIONAL          :: current_level
83
84     INTEGER                                  :: ji,jj
85     INTEGER                                  :: ndb, nde, ndey, idx, idx_var, maxvar
86     LOGICAL                                  :: ladvan
87     INTEGER                                  :: ich
88     INTEGER                                  :: fsize,sizemax
89     CHARACTER(LEN=FM_FIELD_SIZE)             :: yrecfm
90     CHARACTER(LEN=4)                         :: suffix
91 #ifdef LOWMEM
92     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
93 #endif
94     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
95     !JUAN CYCCL3
96     INTEGER                        :: JPHEXT
97
98     ilu = klu
99
100     CALL FMREADLFIN1(klu,'JPHEXT',JPHEXT,iresp)
101     IF (iresp /= 0) JPHEXT=1
102     ! First check if IMAX,JMAX,KMAX exist in LFI file
103     ! to handle 3D, 2D variables -> update IDIMX,IDIMY,IDIMZ
104     CALL FMREADLFIN1(klu,'IMAX',IDIMX,iresp)
105     IF (iresp == 0) IDIMX = IDIMX+2*JPHEXT  ! IMAX + 2*JPHEXT
106     !
107     CALL FMREADLFIN1(klu,'JMAX',IDIMY,iresp)
108     IF (iresp == 0) IDIMY = IDIMY+2*JPHEXT  ! JMAX + 2*JPHEXT
109     !
110     CALL FMREADLFIN1(ilu,'KMAX',IDIMZ,iresp)
111     IF (iresp == 0) IDIMZ = IDIMZ+2  ! KMAX + 2*JPVEXT
112     GUSEDIM = (IDIMX*IDIMY > 0)
113     IF (GUSEDIM) THEN
114       PRINT *,'MESONH 3D, 2D articles DIMENSIONS used :'
115       PRINT *,'DIMX =',IDIMX
116       PRINT *,'DIMY =',IDIMY
117       PRINT *,'DIMZ =',IDIMZ ! IDIMZ may be equal to 0 (PGD files)
118     ELSE
119       PRINT *,'BEWARE : ALL MesoNH arrays are handled as 1D arrays !'
120     END IF
121
122     sizemax = 0
123
124     IF (present(current_level)) THEN
125       write(suffix,'(I4.4)') current_level
126     ElSE
127       suffix=''
128     END IF
129
130     ! Phase 1 : build articles list to convert.
131     !
132     !    Pour l'instant tous les articles du fichier LFI sont
133     !    convertis. On peut modifier cette phase pour prendre en
134     !    compte un sous-ensemble d'article (liste definie par
135     !    l'utilisateur par exemple)  
136     !
137     IF (LEN_TRIM(hvarlist) > 0) THEN
138 #ifndef LOWMEM
139       IF(.NOT.ALLOCATED(lfiart)) ALLOCATE(lfiart(nbvar_tbr+nbvar_calc))
140 #endif
141       ALLOCATE(tpreclist(nbvar_tbr+nbvar_calc))
142       DO ji=1,nbvar_tbr+nbvar_calc
143         tpreclist(ji)%found  = .FALSE.
144         tpreclist(ji)%calc   = .FALSE. !By default variables are not computed from others
145         tpreclist(ji)%tbw    = .TRUE.  !By default variables are written
146         tpreclist(ji)%tbr    = .TRUE.  !By default variables are written
147         tpreclist(ji)%src(:) = -1
148         tpreclist(ji)%tgt    = -1
149       END DO
150
151        ! A variable list is provided with -v var1,...
152        ndb  = 1
153        idx_var = 1
154        DO ji=1,nbvar_tbw
155           nde = INDEX(TRIM(hvarlist(ndb:)),',')
156           yrecfm = hvarlist(ndb:ndb+nde-2)
157
158           !Detect operations on variables (only + is supported now)
159           ndey = INDEX(TRIM(yrecfm),'=')
160           idx = 1
161           IF (ndey /= 0) THEN
162             var_calc = yrecfm(1:ndey-1)
163             DO WHILE (ndey /= 0)
164               IF (idx>MAXRAW) THEN
165                 print *,'Error: MAXRAW exceeded (too many raw variables for 1 computed one)'
166                 STOP
167               END IF
168               yrecfm = yrecfm(ndey+1:)
169               ndey = INDEX(TRIM(yrecfm),'+')
170               IF (ndey /= 0) THEN
171                 var_raw(idx) = yrecfm(1:ndey-1)
172               ELSE
173                 var_raw(idx) = TRIM(yrecfm)
174               END IF
175               idx = idx + 1
176             END DO
177
178             tpreclist(idx_var)%name = trim(var_calc)
179             tpreclist(idx_var)%calc = .TRUE.
180             tpreclist(idx_var)%tbw  = .TRUE.
181             tpreclist(idx_var)%tbr  = .FALSE.
182             idx_var=idx_var+1
183             DO jj = 1, idx-1
184               tpreclist(idx_var-jj)%src(jj) = idx_var
185               tpreclist(idx_var)%name = trim(var_raw(jj))
186               tpreclist(idx_var)%calc = .FALSE.
187               tpreclist(idx_var)%tbw  = .FALSE.
188               tpreclist(idx_var)%tbr  = .TRUE.
189               tpreclist(idx_var)%tgt  = idx_var-jj
190               idx_var=idx_var+1
191             END DO
192
193           ELSE
194             tpreclist(idx_var)%name = trim(yrecfm)
195             tpreclist(idx_var)%calc = .FALSE.
196             tpreclist(idx_var)%tbw  = .TRUE.
197             idx_var=idx_var+1
198
199           END IF
200
201           ndb = nde+ndb
202        END DO
203
204 !TODO: merge loop?
205        DO ji=1,nbvar_tbr+nbvar_calc
206           IF (tpreclist(ji)%calc) CYCLE
207           yrecfm = TRIM(tpreclist(ji)%name)
208           CALL LFINFO(iresp,ilu,trim(yrecfm)//trim(suffix),ileng,ipos)
209           
210           IF (iresp /= 0 .OR. ileng == 0) THEN
211              PRINT *,'Article ',TRIM(yrecfm), ' not found!'
212              tpreclist(ji)%found = .FAlSE.
213              tpreclist(ji)%tbw   = .FAlSE.
214              tpreclist(ji)%tbr   = .FAlSE.
215           ELSE
216              tpreclist(ji)%found = .TRUE.
217              ! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
218              IF (ileng > sizemax) sizemax = ileng        
219 #ifndef LOWMEM
220              ALLOCATE(lfiart(ji)%iwtab(ileng))
221 #endif
222           end IF
223        END DO
224
225        maxvar = nbvar_tbr+nbvar_calc
226
227 DO ji=1,nbvar_tbr+nbvar_calc
228   print *,ji,'name=',trim(tpreclist(ji)%name),' calc=',tpreclist(ji)%calc,' tbw=',tpreclist(ji)%tbw,&
229           ' tbr=',tpreclist(ji)%tbr,' found=',tpreclist(ji)%found
230 END DO
231
232     ELSE
233        ! Entire file is converted
234 #ifndef LOWMEM
235        IF(.NOT.ALLOCATED(lfiart)) ALLOCATE(lfiart(nbvar_lfi))
236 #endif
237        ALLOCATE(tpreclist(nbvar_lfi))
238        DO ji=1,nbvar_lfi
239          tpreclist(ji)%calc   = .FALSE. !By default variables are not computed from others
240          tpreclist(ji)%tbw    = .TRUE.  !By default variables are written
241          tpreclist(ji)%src(:) = -1
242        END DO
243
244        CALL LFIPOS(iresp,ilu)
245        ladvan = .TRUE.
246        
247        DO ji=1,nbvar_lfi
248           CALL LFICAS(iresp,ilu,yrecfm,ileng,ipos,ladvan)
249           ! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
250           tpreclist(ji)%name = trim(yrecfm)
251           tpreclist(ji)%found  = .TRUE.
252           IF (ileng > sizemax) sizemax = ileng        
253 #ifndef LOWMEM       
254           ALLOCATE(lfiart(ji)%iwtab(ileng))
255 #endif
256        END DO
257        maxvar = nbvar_lfi
258     END IF
259
260     kbuflen = sizemax
261
262 #ifdef LOWMEM
263     WRITE(*,'("Taille maximale du buffer :",f10.3," Mio")') sizemax*8./1048576.
264     ALLOCATE(iwork(sizemax))
265 #endif
266     
267     ! Phase 2 : Extract comments and dimensions for valid articles.
268     !           Infos are put in tpreclist.
269     CALL init_dimCDF()
270     DO ji=1,maxvar
271        IF (tpreclist(ji)%calc .OR. .NOT.tpreclist(ji)%found) CYCLE
272
273        yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
274        CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
275 #ifdef LOWMEM
276        CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
277        tpreclist(ji)%TYPE = get_ftype(yrecfm)               
278        tpreclist(ji)%grid = iwork(1)
279
280        ALLOCATE(character(len=iwork(2)) :: tpreclist(ji)%comment)
281        DO jj=1,iwork(2)
282           ich = iwork(2+jj)
283           tpreclist(ji)%comment(jj:jj) = CHAR(ich)
284        END DO
285        fsize = ileng-(2+iwork(2))
286 #else
287        CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
288        tpreclist(ji)%TYPE = get_ftype(yrecfm)               
289        tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
290
291        ALLOCATE(character(len=lfiart(ji)%iwtab(2)) :: tpreclist(ji)%comment)
292        DO jj=1,lfiart(ji)%iwtab(2)
293           ich = lfiart(ji)%iwtab(2+jj)
294           tpreclist(ji)%comment(jj:jj) = CHAR(ich)
295        END DO
296        fsize = ileng-(2+lfiart(ji)%iwtab(2))
297 #endif
298        tpreclist(ji)%dim=>get_dimCDF(fsize)
299     END DO
300
301     !Complete info for calculated variables
302     IF (nbvar_calc>0) THEN
303     DO ji=1,maxvar
304        IF (.NOT.tpreclist(ji)%calc) CYCLE
305        tpreclist(ji)%TYPE = tpreclist(tpreclist(ji)%src(1))%TYPE
306        tpreclist(ji)%grid = tpreclist(tpreclist(ji)%src(1))%grid
307        tpreclist(ji)%dim  => tpreclist(tpreclist(ji)%src(1))%dim
308
309 !TODO: cleaner length!
310        ALLOCATE(character(len=256) :: tpreclist(ji)%comment)
311        tpreclist(ji)%comment='Constructed from'
312        jj = 1
313        DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
314          tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' '//trim(tpreclist(tpreclist(ji)%src(jj))%name)
315          IF (jj<MAXRAW .AND. tpreclist(ji)%src(jj+1)>0) THEN
316            tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' +'
317          END IF
318          jj=jj+1
319        END DO
320     END DO
321     END IF
322
323   
324     PRINT *,'Nombre de dimensions = ', size_dimCDF()
325 #ifdef LOWMEM
326     DEALLOCATE(iwork)
327 #endif
328   END SUBROUTINE parse_lfi
329   
330   SUBROUTINE read_data_lfi(klu, hvarlist, nbvar, tpreclist, kbuflen, current_level)
331     INTEGER, INTENT(IN)                    :: klu
332     INTEGER, INTENT(INOUT)                 :: nbvar
333     CHARACTER(LEN=*), intent(IN)           :: hvarlist
334     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
335     INTEGER, INTENT(IN)                    :: kbuflen
336     INTEGER, INTENT(IN), OPTIONAL          :: current_level
337
338     INTEGER                                  :: ji,jj
339     INTEGER                                  :: ndb, nde
340     LOGICAL                                  :: ladvan
341     INTEGER                                  :: ich
342     INTEGER                                  :: fsize,sizemax
343     CHARACTER(LEN=FM_FIELD_SIZE)             :: yrecfm
344     CHARACTER(LEN=4)                         :: suffix
345 #ifdef LOWMEM
346     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
347 #endif
348     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
349     CHARACTER(LEN=FM_FIELD_SIZE)             :: var_calc
350     CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
351
352     ilu = klu
353
354     IF (present(current_level)) THEN
355       write(suffix,'(I4.4)') current_level
356     ElSE
357       suffix=''
358     END IF
359
360 #ifdef LOWMEM
361     ALLOCATE(iwork(kbuflen))
362 #endif
363
364     DO ji=1,nbvar
365        IF (.NOT.tpreclist(ji)%tbr) CYCLE
366        yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
367        CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
368 #ifdef LOWMEM
369        CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
370        tpreclist(ji)%grid = iwork(1)
371 #else
372        CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
373        tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
374 #endif
375     END DO
376
377 #ifdef LOWMEM
378     DEALLOCATE(iwork)
379 #endif
380   END SUBROUTINE read_data_lfi
381
382   SUBROUTINE HANDLE_ERR(status,line)
383     INTEGER :: status,line
384
385     IF (status /= NF90_NOERR) THEN
386        PRINT *, 'line ',line,': ',NF90_STRERROR(status)
387        STOP
388     END IF
389   END SUBROUTINE HANDLE_ERR
390
391   SUBROUTINE def_ncdf(tpreclist,nbvar,oreduceprecision,kcdf_id,omerge,ocompress,compress_level)
392     TYPE(workfield),DIMENSION(:),INTENT(INOUT) :: tpreclist
393     INTEGER,                     INTENT(IN) :: nbvar
394     LOGICAL,                     INTENT(IN) :: oreduceprecision
395     INTEGER,                     INTENT(OUT):: kcdf_id
396     LOGICAL,                     INTENT(IN) :: omerge
397     LOGICAL,                     INTENT(IN) :: ocompress
398     INTEGER,                     INTENT(IN) :: compress_level
399
400     INTEGER :: status
401     INTEGER :: ji
402     TYPE(dimCDF), POINTER :: tzdim
403     INTEGER               :: invdims
404     INTEGER               :: type_float
405     INTEGER, DIMENSION(10) :: ivdims
406     CHARACTER(LEN=20)     :: ycdfvar
407
408
409     IF (oreduceprecision) THEN
410       type_float = NF90_REAL
411     ELSE
412       type_float = NF90_DOUBLE
413     END IF
414
415       ! global attributes
416       status = NF90_PUT_ATT(kcdf_id,NF90_GLOBAL,'Title',VERSION_ID)
417       IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
418
419     ! define DIMENSIONS
420     tzdim=>first_DimCDF()
421     DO WHILE(ASSOCIATED(tzdim))
422       IF (tzdim%create) THEN
423         status = NF90_DEF_DIM(kcdf_id,tzdim%name,tzdim%len,tzdim%id)
424         IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
425       END IF
426       tzdim=>tzdim%next
427     END DO
428
429     PRINT *,'------------- NetCDF DEFINITION ---------------'
430
431     ! define VARIABLES and ATTRIBUTES
432     DO ji=1,nbvar
433        IF (.NOT.tpreclist(ji)%tbw) CYCLE
434
435        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
436          IF (tpreclist(ji)%dim%create) THEN
437            invdims   = 1
438            ivdims(1) = tpreclist(ji)%dim%id
439          ELSE
440            invdims = tpreclist(ji)%dim%ndims
441            IF(omerge) invdims=invdims+1 !when merging variables from LFI splitted files
442            SELECT CASE(invdims)
443            CASE(2)
444               ivdims(1)=ptdimx%id
445               ivdims(2)=ptdimy%id
446            CASE(3)
447               ivdims(1)=ptdimx%id
448               ivdims(2)=ptdimy%id
449               ivdims(3)=ptdimz%id
450            CASE(12)
451               ivdims(1)=ptdimx%id
452               ivdims(2)=ptdimz%id
453               invdims = 2 ! on retablit la bonne valeur du nbre de dimension
454            CASE default
455              PRINT *,'Fatal error in NetCDF dimension definition'
456              STOP
457            END SELECT
458          END IF
459        ELSE
460          ! scalar variables
461           invdims   = 0
462           ivdims(1) = 0 ! ignore dans ce cas
463        END IF
464        
465        ! Variables definition
466
467        !! NetCDF n'aime pas les '%' dans le nom des variables
468        !! "%" remplaces par '__' 
469        ycdfvar = str_replace(tpreclist(ji)%name,'%','__')
470        !! ni les '.' remplaces par '--'
471        ycdfvar = str_replace(ycdfvar,'.','--')
472
473        SELECT CASE(tpreclist(ji)%TYPE)
474        CASE (TEXT)
475 !          PRINT *,'TEXT : ',tpreclist(ji)%name
476           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_CHAR,&
477                    ivdims(:invdims),tpreclist(ji)%id)
478           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
479
480        CASE (INT,BOOL)
481 !          PRINT *,'INT,BOOL : ',tpreclist(ji)%name
482           status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_INT,&
483                    ivdims(:invdims),tpreclist(ji)%id)
484           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
485
486        CASE(FLOAT)
487 !          PRINT *,'FLOAT : ',tpreclist(ji)%name
488           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
489                    ivdims(:invdims),tpreclist(ji)%id)
490           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
491
492           
493        CASE default
494           PRINT *,'ATTENTION : ',TRIM(tpreclist(ji)%name),' est de&
495                & TYPE inconnu --> force a REAL'
496           status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
497                    ivdims(:invdims),tpreclist(ji)%id)
498           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
499           
500
501        END SELECT
502
503        ! Compress data (costly operation for the CPU)
504        IF (ocompress .AND. invdims>0) THEN
505          status = NF90_DEF_VAR_DEFLATE(kcdf_id,tpreclist(ji)%id,1,1,compress_level)
506          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
507        END IF
508
509        ! GRID attribute definition
510        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id,'GRID',tpreclist(ji)%grid)
511        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
512
513        ! COMMENT attribute definition
514        status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id,'COMMENT',trim(tpreclist(ji)%comment))
515        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
516
517     END DO
518     
519     status = NF90_ENDDEF(kcdf_id)
520     IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
521
522     
523   END SUBROUTINE def_ncdf
524
525   SUBROUTINE fill_ncdf(klu,kcdf_id,tpreclist,knaf,kbuflen,current_level)
526     INTEGER,                      INTENT(IN):: klu
527     INTEGER,                      INTENT(IN):: kcdf_id
528     TYPE(workfield), DIMENSION(:),INTENT(IN):: tpreclist
529     INTEGER,                      INTENT(IN):: knaf
530     INTEGER,                      INTENT(IN):: kbuflen
531     INTEGER, INTENT(IN), OPTIONAL           :: current_level
532
533 #ifdef LOWMEM
534     INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
535 #endif
536     INTEGER                                  :: ji,jj
537     INTEGER,DIMENSION(:),ALLOCATABLE         :: itab
538     REAL   (KIND=8),DIMENSION(:),ALLOCATABLE :: xtab
539     CHARACTER, DIMENSION(:), ALLOCATABLE     :: ytab
540     INTEGER                                  :: status
541     INTEGER                                  :: extent, ndims
542     INTEGER                                  :: ich
543     INTEGER                                  :: src
544     INTEGER                                  :: level
545     INTEGER(KIND=LFI_INT)                    :: iresp,ilu,ileng,ipos
546     CHARACTER(LEN=4)                         :: suffix
547
548     !
549     ilu = klu
550     !
551
552     IF (present(current_level)) THEN
553       write(suffix,'(I4.4)') current_level
554       level = current_level
555     ElSE
556       suffix=''
557       level = 1
558     END IF
559
560 #if LOWMEM
561     ALLOCATE(iwork(kbuflen))
562 #endif
563     ALLOCATE(itab(kbuflen))
564     ALLOCATE(xtab(kbuflen))
565
566     DO ji=1,knaf
567        IF (.NOT.tpreclist(ji)%tbw) CYCLE
568
569 #if LOWMEM
570        CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
571        CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
572 #endif
573        IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
574           extent = tpreclist(ji)%dim%len
575           ndims = tpreclist(ji)%dim%ndims
576        ELSE
577           extent = 1
578           ndims = 0
579        END IF
580
581        SELECT CASE(tpreclist(ji)%TYPE)
582        CASE (INT,BOOL)
583 #if LOWMEM
584 ***
585 print *,'lowmem: not tested!!!!!' (to be compared with no low mem version)
586          itab(1:extent) = iwork(3+iwork(2):)
587 #else
588          IF (.NOT.tpreclist(ji)%calc) THEN
589            itab(1:extent) = lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):)
590          ELSE
591            src=tpreclist(ji)%src(1)
592            xtab(1:extent) = lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
593            jj = 2
594            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
595              src=tpreclist(ji)%src(jj)
596              xtab(1:extent) = xtab(1:extent) + lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
597              jj=jj+1
598            END DO
599          END IF
600          itab(1:extent) = lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):)
601 #endif
602 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
603          SELECT CASE(ndims)
604          CASE (0)
605            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,itab(1))
606          CASE (1)
607            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,itab(1:extent),count=(/extent/))
608          CASE (2)
609            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(itab,(/ptdimx%len,ptdimy%len/)), &
610                                  start = (/1,1,level/) )
611          CASE (3)
612            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(itab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
613          CASE DEFAULT
614            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
615          END SELECT
616          
617        CASE (FLOAT)
618 #if LOWMEM
619 ***
620 print *,'lowmem: not tested!!!!!' (to be compared with no low mem version)
621          xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
622 #else
623          IF (.NOT.tpreclist(ji)%calc) THEN
624            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
625          ELSE
626            src=tpreclist(ji)%src(1)
627            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
628            jj = 2
629            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
630              src=tpreclist(ji)%src(jj)
631              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
632              jj=jj+1
633            END DO
634          END IF
635 #endif
636 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
637          SELECT CASE(ndims)
638          CASE (0)
639            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,xtab(1))
640          CASE (1)
641            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,xtab(1:extent),count=(/extent/))
642          CASE (2)
643            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
644                                  start = (/1,1,level/) )
645          CASE (3)
646            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
647          CASE DEFAULT
648            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
649          END SELECT
650
651        CASE (TEXT)
652          ALLOCATE(ytab(extent))
653          DO jj=1,extent
654 #if LOWMEM
655 ***
656 print *,'lowmem: not tested!!!!!' (to be compared with no low mem version)
657            ich = iwork(2+iwork(2)+jj)
658 #else
659            ich = lfiart(ji)%iwtab(2+lfiart(ji)%iwtab(2)+jj)
660 #endif
661            ytab(jj) = CHAR(ich)
662          END DO
663          status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,ytab,count=(/extent/))
664          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
665          DEALLOCATE(ytab)
666
667        CASE default
668 #if LOWMEM
669 ***
670 print *,'lowmem: not tested!!!!!' (to be compared with no low mem version)
671          xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
672 #else         
673          IF (.NOT.tpreclist(ji)%calc) THEN
674            xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
675          ELSE
676            src=tpreclist(ji)%src(1)
677            xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
678            jj = 2
679            DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
680              src=tpreclist(ji)%src(jj)
681              xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
682              jj=jj+1
683            END DO
684          END IF
685 #endif
686 !TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
687          SELECT CASE(ndims)
688          CASE (0)
689            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,xtab(1))
690          CASE (1)
691            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,xtab(1:extent),count=(/extent/))
692          CASE (2)
693            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
694                                  start = (/1,1,level/) )
695          CASE (3)
696            status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
697          CASE DEFAULT
698            print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
699          END SELECT
700          IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
701
702        END SELECT
703
704     END DO
705     DEALLOCATE(itab,xtab)
706 #if LOWMEM
707     DEALLOCATE(iwork)
708 #endif 
709   END SUBROUTINE fill_ncdf
710
711   SUBROUTINE parse_cdf(kcdf_id,tpreclist,kbuflen)
712     INTEGER, INTENT(IN)                    :: kcdf_id
713     TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
714     INTEGER, INTENT(OUT)                   :: kbuflen
715
716
717     INTEGER :: status
718     INTEGER :: nvars, var_id
719     INTEGER :: jdim
720     INTEGER :: sizemax
721     INTEGER :: itype
722     INTEGER, DIMENSION(10) :: idim_id
723     INTEGER :: icomlen,idimlen,idims,idimtmp
724     
725     status = NF90_INQUIRE(kcdf_id, nvariables = nvars)
726     IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
727     ALLOCATE(tpreclist(nvars))
728
729     sizemax = 0
730
731     CALL init_dimCDF()
732     
733     ! Parcours de toutes les variables et extraction des infos
734     !      - nom de dimension
735     !      - dimension, etendue
736     !      - attributs
737     DO var_id = 1, nvars
738        ! Pour la forme
739        tpreclist(var_id)%id = var_id  
740        
741        ! Nom, type et dimensions de la variable
742        status = NF90_INQUIRE_VARIABLE(kcdf_id, var_id, name = tpreclist(var_id)%name, xtype = itype, ndims = idims, &
743                                       dimids = idim_id)
744        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
745        
746        SELECT CASE(itype)
747        CASE(NF90_CHAR)
748           tpreclist(var_id)%TYPE = TEXT
749        CASE(NF90_INT)
750           tpreclist(var_id)%TYPE = INT
751        CASE(NF90_FLOAT,NF90_DOUBLE)
752           tpreclist(var_id)%TYPE = FLOAT
753        CASE default 
754           PRINT *, 'Attention : variable ',TRIM(tpreclist(var_id)&
755                & %name), ' a un TYPE non reconnu par le convertisseur.'
756           PRINT *, '--> TYPE force a REAL(KIND 8) dans LFI !'
757        END SELECT
758       
759        IF (idims == 0) THEN
760           ! variable scalaire
761           NULLIFY(tpreclist(var_id)%dim)
762           idimlen = 1
763        ELSE
764           ! infos sur dimensions
765           idimlen = 1
766           DO jdim=1,idims
767             status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
768             IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
769             idimlen = idimlen*idimtmp
770           END DO
771           
772           tpreclist(var_id)%dim=>get_dimCDF(idimlen)
773           ! seul le champ 'len' de dimCDF sera utilise par la suite
774        END IF
775        
776        ! GRID et COMMENT attributes
777        status = NF90_GET_ATT(kcdf_id,var_id,'GRID',tpreclist(var_id)%grid)
778        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
779
780        status = NF90_INQUIRE_ATTRIBUTE(kcdf_id,var_id,'COMMENT',len = icomlen)
781        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
782        
783        ALLOCATE(character(len=icomlen) :: tpreclist(var_id)%comment)
784        status = NF90_GET_ATT(kcdf_id,var_id,'COMMENT',tpreclist(var_id)%comment)
785        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
786
787        
788        IF (sizemax < icomlen+idimlen) sizemax = icomlen+idimlen 
789
790     END DO
791     
792     kbuflen = sizemax
793
794   END SUBROUTINE parse_cdf
795
796   SUBROUTINE build_lfi(kcdf_id,klu,tpreclist,kbuflen)
797     INTEGER,                       INTENT(IN) :: kcdf_id 
798     INTEGER,                       INTENT(IN) :: klu
799     TYPE(workfield), DIMENSION(:), INTENT(IN) :: tpreclist
800     INTEGER,                       INTENT(IN) :: kbuflen
801     
802     INTEGER :: status
803     INTEGER :: ivar,jj
804     INTEGER(KIND=8), DIMENSION(:), POINTER  :: iwork
805     INTEGER(KIND=8), DIMENSION(:), POINTER  :: idata
806     REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: xtab
807     INTEGER,      DIMENSION(:), ALLOCATABLE :: itab
808     CHARACTER,    DIMENSION(:), ALLOCATABLE :: ytab
809     CHARACTER(LEN=FM_FIELD_SIZE)            :: yrecfm
810
811     INTEGER :: iartlen, idlen, icomlen
812     INTEGER(KIND=LFI_INT) :: iresp,ilu,iartlen8
813
814     ! Un article LFI est compose de :
815     !   - 1 entier identifiant le numero de grille
816     !   - 1 entier contenant la taille du commentaire
817     !   - le commentaire code en entier 64 bits
818     !   - les donnees proprement dites
819
820     PRINT *,'Taille buffer = ',2+kbuflen
821
822     ALLOCATE(iwork(2+kbuflen))
823     ALLOCATE(itab(2+kbuflen))
824     ALLOCATE(xtab(2+kbuflen))
825
826     DO ivar=1,SIZE(tpreclist)
827        icomlen = LEN(tpreclist(ivar)%comment)
828
829        ! traitement Grille et Commentaire
830        iwork(1) = tpreclist(ivar)%grid
831        iwork(2) = icomlen
832        DO jj=1,iwork(2)
833           iwork(2+jj)=ICHAR(tpreclist(ivar)%comment(jj:jj))
834        END DO
835
836        IF (ASSOCIATED(tpreclist(ivar)%dim)) THEN
837           idlen = tpreclist(ivar)%dim%len
838        ELSE 
839           idlen = 1
840        END IF
841        
842        iartlen = 2+icomlen+idlen
843        idata=>iwork(3+icomlen:iartlen)
844
845
846        SELECT CASE(tpreclist(ivar)%TYPE)
847        CASE(INT,BOOL)
848           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id,itab)
849           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
850
851 !          PRINT *,'INT,BOOL --> ',tpreclist(ivar)%name,',len = ',idlen
852           idata(1:idlen) = itab(1:idlen)
853
854        CASE(FLOAT)
855           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id,xtab)
856           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
857           
858 !          PRINT *,'FLOAT    --> ',tpreclist(ivar)%name,',len = ',idlen
859           ! La ligne suivante ne pose aucun pb sur Cray alors que sur
860           ! fuji, elle genere une erreur d'execution
861 !          idata(1:idlen) = TRANSFER(xtab(1:idlen),(/ 0_8 /))
862           
863           ! la correction pour Fuji (valable sur CRAY) est :
864           idata(1:idlen) = TRANSFER(xtab,(/ 0_8 /),idlen)
865
866 !          IF (idlen < 10) PRINT *,'xtab = ',xtab(1:idlen)
867
868        CASE(TEXT)
869           ALLOCATE(ytab(idlen))
870           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id,ytab)
871           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
872
873 !          PRINT *,'TEXT -->     ',tpreclist(ivar)%name,',len = ',idlen
874
875           DO jj=1,idlen
876              idata(jj) = ICHAR(ytab(jj))
877           END DO
878           
879           DEALLOCATE(ytab)
880
881        CASE default
882           status = NF90_GET_VAR(kcdf_id,tpreclist(ivar)%id,xtab)
883           IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
884
885           PRINT *,'Default (ERROR) -->',tpreclist(ivar)%name,',len = ',idlen
886           idata(1:idlen) = TRANSFER(xtab,(/ 0_8 /),idlen)
887
888        END SELECT
889        
890        ! Attention restoration des '%' dans le nom des champs LFI
891        yrecfm = str_replace(tpreclist(ivar)%name,'__','%')
892        ! et des '.'
893        yrecfm = str_replace(yrecfm,'--','.')
894        ilu = klu
895        iartlen8 = iartlen
896        CALL LFIECR(iresp,ilu,yrecfm,iwork,iartlen8)
897
898     END DO
899     DEALLOCATE(iwork,itab,xtab)
900
901   END SUBROUTINE build_lfi
902
903   SUBROUTINE OPEN_FILES(hinfile,houtfile,olfi2cdf,olfilist,ohdf5,kcdf_id,klu,knaf)
904     LOGICAL,          INTENT(IN)  :: olfi2cdf, olfilist, ohdf5
905     CHARACTER(LEN=*), INTENT(IN)  :: hinfile
906     CHARACTER(LEN=*), INTENT(IN)  :: houtfile
907     INTEGER         , INTENT(OUT) :: kcdf_id,klu,knaf
908
909     INTEGER                     :: extindex
910     INTEGER(KIND=LFI_INT)       :: ilu,iresp,iverb,inap,inaf
911     INTEGER                     :: status
912     CHARACTER(LEN=4)            :: ypextsrc, ypextdest
913     LOGICAL                     :: fexist
914     INTEGER                     :: omode
915
916     iverb = 0
917     ilu   = 11
918
919     CALL init_sysfield()
920
921     IF (olfi2cdf) THEN 
922        ! Cas LFI -> NetCDF
923        CALL LFIOUV(iresp,ilu,ltrue,hinfile,'OLD',lfalse&
924             & ,lfalse,iverb,inap,inaf)
925
926        IF (olfilist) THEN
927           CALL LFILAF(iresp,ilu,lfalse)
928           CALL LFIFER(iresp,ilu,'KEEP')
929           return
930        end IF
931
932        IF (ohdf5) THEN
933           status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_NETCDF4), kcdf_id)
934        ELSE
935           status = NF90_CREATE(houtfile, IOR(NF90_CLOBBER,NF90_64BIT_OFFSET), kcdf_id) 
936        end IF
937        
938        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
939
940        status = NF90_SET_FILL(kcdf_id,NF90_NOFILL,omode)
941        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
942 !!$       SELECT CASE(omode)
943 !!$       CASE (NF90_FILL)
944 !!$          PRINT *,'Ancien mode : NF90_FILL'
945 !!$       CASE (NF90_NOFILL)
946 !!$          PRINT *,'Ancien mode : NF90_NOFILL'
947 !!$       CASE default
948 !!$          PRINT *, 'Ancien mode : inconnu'
949 !!$       END SELECT
950        
951     ELSE
952        ! Cas NetCDF -> LFI
953        status = NF90_OPEN(hinfile,NF90_NOWRITE,kcdf_id)
954        IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
955        
956        inap = 100
957        CALL LFIOUV(iresp,ilu,ltrue,houtfile,'NEW'&
958             & ,lfalse,lfalse,iverb,inap,inaf)
959     END IF
960
961     klu  = ilu
962     knaf = inaf
963
964     PRINT *,'--> Fichier converti : ', houtfile
965
966   END SUBROUTINE OPEN_FILES
967   
968   SUBROUTINE CLOSE_FILES(klu,kcdf_id)
969     INTEGER, INTENT(IN) :: klu, kcdf_id
970     
971     INTEGER(KIND=LFI_INT) :: iresp,ilu
972     INTEGER               :: status
973
974     ilu = klu
975     ! close LFI file
976     CALL LFIFER(iresp,ilu,'KEEP')
977
978     ! close NetCDF file
979     status = NF90_CLOSE(kcdf_id)
980     IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
981     
982   END SUBROUTINE CLOSE_files
983   
984 END MODULE mode_util