Newer
Older
MODULE mode_util
USE MODE_FIELDTYPE
USE mode_dimlist
WAUTELET Philippe
committed
USE mode_options
USE netcdf
INTEGER,PARAMETER :: MAXRAW=10
WAUTELET Philippe
committed
INTEGER,PARAMETER :: MAXLEN=512
WAUTELET Philippe
committed
INTEGER,PARAMETER :: UNDEFINED = -1, READING = 1, WRITING = 2
INTEGER,PARAMETER :: UNKNOWN_FORMAT = -1, NETCDF_FORMAT = 1, LFI_FORMAT = 2
TYPE filestruct
INTEGER :: lun_id ! Logical ID of file
INTEGER :: format = UNKNOWN_FORMAT ! NETCDF, LFI
INTEGER :: status = UNDEFINED ! Opened for reading or writing
INTEGER :: var_id ! Position of the variable in the workfield structure
LOGICAL :: opened = .false.
END TYPE filestruct
TYPE filelist_struct
INTEGER :: nbfiles = 0
! TYPE(filestruct),DIMENSION(:),ALLOCATABLE :: files
TYPE(filestruct),DIMENSION(MAXFILES) :: files
END TYPE filelist_struct
WAUTELET Philippe
committed
TYPE workfield
CHARACTER(LEN=FM_FIELD_SIZE) :: name ! nom du champ
INTEGER :: TYPE ! type (entier ou reel)
CHARACTER(LEN=:), POINTER :: comment
LOGICAL :: found ! T if found in the input file
LOGICAL :: calc ! T if computed from other variables
LOGICAL :: tbw ! to be written or not
LOGICAL :: tbr ! to be read or not
INTEGER,DIMENSION(MAXRAW) :: src ! List of variables used to compute the variable (needed only if calc=.true.)
INTEGER :: tgt ! Target: id of the variable that use it (calc variable)
END TYPE workfield
#ifndef LOWMEM
TYPE lfidata
INTEGER(KIND=8), DIMENSION(:), POINTER :: iwtab
END TYPE lfidata
TYPE(lfidata), DIMENSION(:), ALLOCATABLE :: lfiart
#endif
LOGICAL(KIND=LFI_INT), PARAMETER :: ltrue = .TRUE.
LOGICAL(KIND=LFI_INT), PARAMETER :: lfalse = .FALSE.
CONTAINS
FUNCTION str_replace(hstr, hold, hnew)
CHARACTER(LEN=*) :: hstr, hold, hnew
CHARACTER(LEN=LEN_TRIM(hstr)+MAX(0,LEN(hnew)-LEN(hold))) :: str_replace
INTEGER :: pos
pos = INDEX(hstr,hold)
IF (pos /= 0) THEN
str_replace = hstr(1:pos-1)//hnew//hstr(pos+LEN(hold):)
ELSE
str_replace = hstr
END IF
END FUNCTION str_replace
SUBROUTINE FMREADLFIN1(klu,hrecfm,kval,kresp)
INTEGER(KIND=LFI_INT), INTENT(IN) :: klu ! logical fortran unit au lfi file
CHARACTER(LEN=*),INTENT(IN) :: hrecfm ! article name to be read
INTEGER, INTENT(OUT) :: kval ! integer value for hrecfm article
INTEGER(KIND=LFI_INT), INTENT(OUT):: kresp! return code null if OK
!
INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE::iwork
INTEGER :: icomlen
INTEGER(KIND=LFI_INT) :: iresp,ilenga,iposex
!
CALL LFINFO(iresp,klu,hrecfm,ilenga,iposex)
IF (iresp /=0 .OR. ilenga == 0) THEN
kresp = -1
kval = 0
ELSE
ALLOCATE(IWORK(ilenga))
CALL LFILEC(iresp,klu,hrecfm,iwork,ilenga)
icomlen = iwork(2)
kval = iwork(3+icomlen)
kresp = iresp
DEALLOCATE(IWORK)
END IF
END SUBROUTINE FMREADLFIN1
WAUTELET Philippe
committed
SUBROUTINE parse_infiles(infiles, nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw, tpreclist, kbuflen, options, icurrent_level)
TYPE(filelist_struct), INTENT(IN) :: infiles
INTEGER, INTENT(IN) :: nbvar_infile, nbvar_tbr, nbvar_calc, nbvar_tbw
TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
INTEGER, INTENT(OUT) :: kbuflen
WAUTELET Philippe
committed
TYPE(option),DIMENSION(:), INTENT(IN) :: options
INTEGER, INTENT(IN), OPTIONAL :: icurrent_level
INTEGER :: ndb, nde, ndey, idx, idx_var, maxvar
INTEGER :: idims, idimtmp, jdim, status, var_id
INTEGER :: comment_size, fsize, sizemax
CHARACTER(LEN=FM_FIELD_SIZE) :: yrecfm
CHARACTER(LEN=4) :: suffix
#ifdef LOWMEM
INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
#endif
INTEGER(KIND=LFI_INT) :: iresp,ilu,ileng,ipos
CHARACTER(LEN=FM_FIELD_SIZE) :: var_calc
CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
INTEGER, DIMENSION(10) :: idim_id
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
IF (infiles%files(1)%format == LFI_FORMAT) THEN
ilu = infiles%files(1)%lun_id
CALL FMREADLFIN1(ilu,'JPHEXT',JPHEXT,iresp)
IF (iresp /= 0) JPHEXT=1
! First check if IMAX,JMAX,KMAX exist in LFI file
! to handle 3D, 2D variables -> update IDIMX,IDIMY,IDIMZ
CALL FMREADLFIN1(ilu,'IMAX',IDIMX,iresp)
IF (iresp == 0) IDIMX = IDIMX+2*JPHEXT ! IMAX + 2*JPHEXT
!
CALL FMREADLFIN1(ilu,'JMAX',IDIMY,iresp)
IF (iresp == 0) IDIMY = IDIMY+2*JPHEXT ! JMAX + 2*JPHEXT
!
CALL FMREADLFIN1(ilu,'KMAX',IDIMZ,iresp)
IF (iresp == 0) IDIMZ = IDIMZ+2 ! KMAX + 2*JPVEXT
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
kcdf_id = infiles%files(1)%lun_id
status = NF90_INQ_DIMID(kcdf_id, "DIMX", idim_id(1))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(1),len = IDIMX)
status = NF90_INQ_DIMID(kcdf_id, "DIMY", idim_id(2))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(2),len = IDIMY)
status = NF90_INQ_DIMID(kcdf_id, "DIMZ", idim_id(3))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(3),len = IDIMZ)
END IF
GUSEDIM = (IDIMX*IDIMY > 0)
IF (GUSEDIM) THEN
PRINT *,'MESONH 3D, 2D articles DIMENSIONS used :'
PRINT *,'DIMX =',IDIMX
PRINT *,'DIMY =',IDIMY
PRINT *,'DIMZ =',IDIMZ ! IDIMZ may be equal to 0 (PGD files)
ELSE
PRINT *,'BEWARE : ALL MesoNH arrays are handled as 1D arrays !'
END IF
sizemax = 0
IF (present(icurrent_level)) THEN
write(suffix,'(I4.4)') icurrent_level
current_level = icurrent_level
ElSE
suffix=''
current_level = -1
END IF
! Phase 1 : build articles list to convert.
!
! Pour l'instant tous les articles du fichier LFI sont
! convertis. On peut modifier cette phase pour prendre en
! compte un sous-ensemble d'article (liste definie par
WAUTELET Philippe
committed
IF (options(OPTVAR)%set) THEN
#ifndef LOWMEM
IF(.NOT.ALLOCATED(lfiart) .AND. infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(nbvar_tbr+nbvar_calc))
#endif
ALLOCATE(tpreclist(nbvar_tbr+nbvar_calc))
DO ji=1,nbvar_tbr+nbvar_calc
tpreclist(ji)%found = .FALSE.
tpreclist(ji)%calc = .FALSE. !By default variables are not computed from others
tpreclist(ji)%tbw = .TRUE. !By default variables are written
tpreclist(ji)%tbr = .TRUE. !By default variables are read
tpreclist(ji)%src(:) = -1
tpreclist(ji)%tgt = -1
END DO
! A variable list is provided with -v var1,...
ndb = 1
idx_var = 1
DO ji=1,nbvar_tbw
WAUTELET Philippe
committed
!crash compiler GCC 4.2.0: nde = INDEX(TRIM(options(OPTVAR)%cvalue(ndb:)),',')
nde = INDEX(TRIM(options(OPTVAR)%cvalue(ndb:len(trim(options(OPTVAR)%cvalue)))),',')
IF (nde == 0) nde = LEN( TRIM(options(OPTVAR)%cvalue(ndb:len(trim(options(OPTVAR)%cvalue)))) ) + 1
WAUTELET Philippe
committed
yrecfm = options(OPTVAR)%cvalue(ndb:ndb+nde-2)
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
!Detect operations on variables (only + is supported now)
ndey = INDEX(TRIM(yrecfm),'=')
idx = 1
IF (ndey /= 0) THEN
var_calc = yrecfm(1:ndey-1)
DO WHILE (ndey /= 0)
IF (idx>MAXRAW) THEN
print *,'Error: MAXRAW exceeded (too many raw variables for 1 computed one)'
STOP
END IF
yrecfm = yrecfm(ndey+1:)
ndey = INDEX(TRIM(yrecfm),'+')
IF (ndey /= 0) THEN
var_raw(idx) = yrecfm(1:ndey-1)
ELSE
var_raw(idx) = TRIM(yrecfm)
END IF
idx = idx + 1
END DO
tpreclist(idx_var)%name = trim(var_calc)
tpreclist(idx_var)%calc = .TRUE.
tpreclist(idx_var)%tbw = .TRUE.
tpreclist(idx_var)%tbr = .FALSE.
idx_var=idx_var+1
DO jj = 1, idx-1
tpreclist(idx_var-jj)%src(jj) = idx_var
tpreclist(idx_var)%name = trim(var_raw(jj))
tpreclist(idx_var)%calc = .FALSE.
tpreclist(idx_var)%tbw = .FALSE.
tpreclist(idx_var)%tbr = .TRUE.
tpreclist(idx_var)%tgt = idx_var-jj
idx_var=idx_var+1
END DO
ELSE
tpreclist(idx_var)%name = trim(yrecfm)
tpreclist(idx_var)%calc = .FALSE.
tpreclist(idx_var)%tbw = .TRUE.
idx_var=idx_var+1
END IF
ndb = nde+ndb
END DO
DO ji=1,nbvar_tbr+nbvar_calc
IF (tpreclist(ji)%calc) CYCLE
yrecfm = TRIM(tpreclist(ji)%name)
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
IF (infiles%files(1)%format == LFI_FORMAT) THEN
CALL LFINFO(iresp,ilu,trim(yrecfm)//trim(suffix),ileng,ipos)
IF (iresp == 0 .AND. ileng /= 0) tpreclist(ji)%found = .true.
leng = ileng
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
status = NF90_INQ_VARID(kcdf_id,trim(yrecfm)//trim(suffix),tpreclist(ji)%id_in)
IF (status == NF90_NOERR) THEN
tpreclist(ji)%found = .true.
status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in,ndims = idims,dimids = idim_id)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
!TODO:useful?
!DUPLICATED
IF (idims == 0) THEN
! variable scalaire
leng = 1
ELSE
! infos sur dimensions
leng = 1
DO jdim=1,idims
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
leng = leng*idimtmp
END DO
END IF
END IF
END IF
IF (.NOT.tpreclist(ji)%found) THEN
PRINT *,'Article ',TRIM(yrecfm), ' not found!'
tpreclist(ji)%tbw = .FAlSE.
tpreclist(ji)%tbr = .FAlSE.
ELSE
! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
#ifndef LOWMEM
!TODO:useful for netcdf?
IF (infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(ji)%iwtab(leng))
maxvar = nbvar_tbr+nbvar_calc
DO ji=1,nbvar_tbr+nbvar_calc
print *,ji,'name=',trim(tpreclist(ji)%name),' calc=',tpreclist(ji)%calc,' tbw=',tpreclist(ji)%tbw,&
' tbr=',tpreclist(ji)%tbr,' found=',tpreclist(ji)%found
END DO
ELSE
! Entire file is converted
#ifndef LOWMEM
IF(.NOT.ALLOCATED(lfiart) .AND. infiles%files(1)%format == LFI_FORMAT) ALLOCATE(lfiart(nbvar_infile))
#endif
ALLOCATE(tpreclist(nbvar_infile))
DO ji=1,nbvar_infile
tpreclist(ji)%calc = .FALSE. !By default variables are not computed from others
tpreclist(ji)%tbw = .TRUE. !By default variables are written
tpreclist(ji)%src(:) = -1
END DO
IF (infiles%files(1)%format == LFI_FORMAT) THEN
CALL LFIPOS(iresp,ilu)
ladvan = .TRUE.
DO ji=1,nbvar_infile
CALL LFICAS(iresp,ilu,yrecfm,ileng,ipos,ladvan)
! PRINT *,'Article ',ji,' : ',TRIM(yrecfm),', longueur = ',ileng
tpreclist(ji)%name = trim(yrecfm)
tpreclist(ji)%found = .TRUE.
IF (ileng > sizemax) sizemax = ileng
#ifndef LOWMEM
ALLOCATE(lfiart(ji)%iwtab(ileng))
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
END DO
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
DO ji=1,nbvar_infile
tpreclist(ji)%id_in = ji
status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in, name = tpreclist(ji)%name, ndims = idims, &
dimids = idim_id)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
! PRINT *,'Article ',ji,' : ',TRIM(tpreclist(ji)%name),', longueur = ',ileng
tpreclist(ji)%found = .TRUE.
!TODO:useful?
!DUPLICATED
IF (idims == 0) THEN
! variable scalaire
leng = 1
ELSE
! infos sur dimensions
leng = 1
DO jdim=1,idims
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
leng = leng*idimtmp
END DO
END IF
IF (leng > sizemax) sizemax = leng
END DO
END IF
maxvar = nbvar_infile
WRITE(*,'("Taille maximale du buffer :",f10.3," Mio")') sizemax*8./1048576.
ALLOCATE(iwork(sizemax))
#endif
! Phase 2 : Extract comments and dimensions for valid articles.
! Infos are put in tpreclist.
CALL init_dimCDF()
DO ji=1,maxvar
IF (tpreclist(ji)%calc .OR. .NOT.tpreclist(ji)%found) CYCLE
IF (infiles%files(1)%format == LFI_FORMAT) THEN
yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
tpreclist(ji)%grid = iwork(1)
comment_size = iwork(2)
CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
comment_size = lfiart(ji)%iwtab(2)
tpreclist(ji)%TYPE = get_ftype(yrecfm,current_level)
ALLOCATE(character(len=comment_size) :: tpreclist(ji)%comment)
DO jj=1,comment_size
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
tpreclist(ji)%comment(jj:jj) = CHAR(ich)
END DO
fsize = ileng-(2+comment_size)
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
! GRID attribute definition
status = NF90_GET_ATT(kcdf_id,tpreclist(ji)%id_in,'GRID',tpreclist(ji)%grid)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
! COMMENT attribute definition
status = NF90_INQUIRE_ATTRIBUTE(kcdf_id,tpreclist(ji)%id_in,'COMMENT',len=comment_size)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
ALLOCATE(character(len=comment_size) :: tpreclist(ji)%comment)
status = NF90_GET_ATT(kcdf_id,tpreclist(ji)%id_in,'COMMENT',tpreclist(ji)%comment)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_INQUIRE_VARIABLE(kcdf_id,tpreclist(ji)%id_in, xtype = itype, ndims = idims, &
dimids = idim_id)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
SELECT CASE(itype)
CASE(NF90_CHAR)
tpreclist(ji)%TYPE = TEXT
CASE(NF90_INT)
tpreclist(ji)%TYPE = INT
CASE(NF90_FLOAT,NF90_DOUBLE)
tpreclist(ji)%TYPE = FLOAT
CASE default
PRINT *, 'Attention : variable ',TRIM(tpreclist(ji)%name), ' a un TYPE non reconnu par le convertisseur.'
PRINT *, '--> TYPE force a REAL(KIND 8) dans LFI !'
END SELECT
!DUPLICATED
IF (idims == 0) THEN
! variable scalaire
leng = 1
ELSE
! infos sur dimensions
leng = 1
DO jdim=1,idims
status = NF90_INQUIRE_DIMENSION(kcdf_id,idim_id(jdim),len = idimtmp)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
leng = leng*idimtmp
END DO
END IF
fsize = leng
END IF
tpreclist(ji)%dim=>get_dimCDF(fsize)
END DO
!Complete info for calculated variables
IF (nbvar_calc>0) THEN
DO ji=1,maxvar
IF (.NOT.tpreclist(ji)%calc) CYCLE
tpreclist(ji)%TYPE = tpreclist(tpreclist(ji)%src(1))%TYPE
tpreclist(ji)%grid = tpreclist(tpreclist(ji)%src(1))%grid
tpreclist(ji)%dim => tpreclist(tpreclist(ji)%src(1))%dim
!TODO: cleaner length!
ALLOCATE(character(len=256) :: tpreclist(ji)%comment)
tpreclist(ji)%comment='Constructed from'
jj = 1
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' '//trim(tpreclist(tpreclist(ji)%src(jj))%name)
IF (jj<MAXRAW .AND. tpreclist(ji)%src(jj+1)>0) THEN
tpreclist(ji)%comment = trim(tpreclist(ji)%comment)//' +'
END IF
jj=jj+1
END DO
END DO
END IF
PRINT *,'Nombre de dimensions = ', size_dimCDF()
#ifdef LOWMEM
DEALLOCATE(iwork)
#endif
WAUTELET Philippe
committed
SUBROUTINE read_data_lfi(infiles, nbvar, tpreclist, kbuflen, current_level)
TYPE(filelist_struct), INTENT(IN) :: infiles
INTEGER, INTENT(INOUT) :: nbvar
TYPE(workfield), DIMENSION(:), POINTER :: tpreclist
INTEGER, INTENT(IN) :: kbuflen
INTEGER, INTENT(IN), OPTIONAL :: current_level
INTEGER :: ji,jj
INTEGER :: ndb, nde
LOGICAL :: ladvan
INTEGER :: ich
INTEGER :: fsize,sizemax
CHARACTER(LEN=FM_FIELD_SIZE) :: yrecfm
CHARACTER(LEN=4) :: suffix
#ifdef LOWMEM
INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
#endif
INTEGER(KIND=LFI_INT) :: iresp,ilu,ileng,ipos
CHARACTER(LEN=FM_FIELD_SIZE) :: var_calc
CHARACTER(LEN=FM_FIELD_SIZE),dimension(MAXRAW) :: var_raw
IF (present(current_level)) THEN
write(suffix,'(I4.4)') current_level
ElSE
suffix=''
END IF
#ifdef LOWMEM
ALLOCATE(iwork(kbuflen))
#endif
DO ji=1,nbvar
IF (.NOT.tpreclist(ji)%tbr) CYCLE
yrecfm = trim(tpreclist(ji)%name)//trim(suffix)
CALL LFINFO(iresp,ilu,yrecfm,ileng,ipos)
#ifdef LOWMEM
CALL LFILEC(iresp,ilu,yrecfm,iwork,ileng)
tpreclist(ji)%grid = iwork(1)
#else
CALL LFILEC(iresp,ilu,yrecfm,lfiart(ji)%iwtab,ileng)
tpreclist(ji)%grid = lfiart(ji)%iwtab(1)
#endif
END DO
#ifdef LOWMEM
DEALLOCATE(iwork)
#endif
END SUBROUTINE read_data_lfi
SUBROUTINE HANDLE_ERR(status,line)
INTEGER :: status,line
IF (status /= NF90_NOERR) THEN
PRINT *, 'line ',line,': ',NF90_STRERROR(status)
END IF
END SUBROUTINE HANDLE_ERR
WAUTELET Philippe
committed
SUBROUTINE def_ncdf(outfiles,tpreclist,nbvar,options)
TYPE(filelist_struct), INTENT(IN) :: outfiles
TYPE(workfield),DIMENSION(:),INTENT(INOUT) :: tpreclist
INTEGER, INTENT(IN) :: nbvar
WAUTELET Philippe
committed
TYPE(option),DIMENSION(:), INTENT(IN) :: options
WAUTELET Philippe
committed
INTEGER :: compress_level, status
WAUTELET Philippe
committed
INTEGER :: idx, ji, nbfiles
INTEGER:: kcdf_id
TYPE(dimCDF), POINTER :: tzdim
INTEGER :: invdims
WAUTELET Philippe
committed
INTEGER :: type_float
INTEGER, DIMENSION(10) :: ivdims
CHARACTER(LEN=20) :: ycdfvar
WAUTELET Philippe
committed
WAUTELET Philippe
committed
IF (options(OPTREDUCE)%set) THEN
WAUTELET Philippe
committed
type_float = NF90_REAL
ELSE
type_float = NF90_DOUBLE
END IF
WAUTELET Philippe
committed
DO ji = 1,nbfiles
WAUTELET Philippe
committed
! global attributes
status = NF90_PUT_ATT(kcdf_id,NF90_GLOBAL,'Title',VERSION_ID)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
WAUTELET Philippe
committed
! define DIMENSIONS
tzdim=>first_DimCDF()
DO WHILE(ASSOCIATED(tzdim))
IF (tzdim%create) THEN
status = NF90_DEF_DIM(kcdf_id,tzdim%name,tzdim%len,tzdim%id)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
END IF
tzdim=>tzdim%next
END DO
END DO
PRINT *,'------------- NetCDF DEFINITION ---------------'
! define VARIABLES and ATTRIBUTES
WAUTELET Philippe
committed
idx = 1
DO ji=1,nbvar
IF (.NOT.tpreclist(ji)%tbw) CYCLE
IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
IF (tpreclist(ji)%dim%create) THEN
invdims = 1
ivdims(1) = tpreclist(ji)%dim%id
ELSE
invdims = tpreclist(ji)%dim%ndims
WAUTELET Philippe
committed
IF(options(OPTMERGE)%set) invdims=invdims+1 !when merging variables from LFI splitted files
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
SELECT CASE(invdims)
CASE(2)
ivdims(1)=ptdimx%id
ivdims(2)=ptdimy%id
CASE(3)
ivdims(1)=ptdimx%id
ivdims(2)=ptdimy%id
ivdims(3)=ptdimz%id
CASE(12)
ivdims(1)=ptdimx%id
ivdims(2)=ptdimz%id
invdims = 2 ! on retablit la bonne valeur du nbre de dimension
CASE default
PRINT *,'Fatal error in NetCDF dimension definition'
STOP
END SELECT
END IF
ELSE
! scalar variables
invdims = 0
ivdims(1) = 0 ! ignore dans ce cas
END IF
! Variables definition
!! NetCDF n'aime pas les '%' dans le nom des variables
!! "%" remplaces par '__'
ycdfvar = str_replace(tpreclist(ji)%name,'%','__')
!! ni les '.' remplaces par '--'
ycdfvar = str_replace(ycdfvar,'.','--')
WAUTELET Philippe
committed
SELECT CASE(tpreclist(ji)%TYPE)
CASE (TEXT)
! PRINT *,'TEXT : ',tpreclist(ji)%name
status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_CHAR,&
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
CASE (INT,BOOL)
! PRINT *,'INT,BOOL : ',tpreclist(ji)%name
status = NF90_DEF_VAR(kcdf_id,ycdfvar,NF90_INT,&
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
CASE(FLOAT)
! PRINT *,'FLOAT : ',tpreclist(ji)%name
WAUTELET Philippe
committed
status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
CASE default
PRINT *,'ATTENTION : ',TRIM(tpreclist(ji)%name),' est de&
& TYPE inconnu --> force a REAL'
WAUTELET Philippe
committed
status = NF90_DEF_VAR(kcdf_id,ycdfvar,type_float,&
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
! Compress data (costly operation for the CPU)
WAUTELET Philippe
committed
IF (options(OPTCOMPRESS)%set .AND. invdims>0) THEN
compress_level = options(OPTCOMPRESS)%ivalue
status = NF90_DEF_VAR_DEFLATE(kcdf_id,tpreclist(ji)%id_out,1,1,compress_level)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
END IF
status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'GRID',tpreclist(ji)%grid)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_PUT_ATT(kcdf_id,tpreclist(ji)%id_out,'COMMENT',trim(tpreclist(ji)%comment))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
WAUTELET Philippe
committed
IF (options(OPTSPLIT)%set) idx = idx + 1
WAUTELET Philippe
committed
DO ji = 1,nbfiles
WAUTELET Philippe
committed
status = NF90_ENDDEF(kcdf_id)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
END DO
WAUTELET Philippe
committed
SUBROUTINE fill_ncdf(infiles,outfiles,tpreclist,knaf,kbuflen,options,current_level)
TYPE(filelist_struct), INTENT(IN):: infiles, outfiles
TYPE(workfield), DIMENSION(:),INTENT(IN):: tpreclist
INTEGER, INTENT(IN):: knaf
WAUTELET Philippe
committed
TYPE(option),DIMENSION(:), INTENT(IN):: options
INTEGER, INTENT(IN), OPTIONAL :: current_level
#ifdef LOWMEM
INTEGER(KIND=8),DIMENSION(:),ALLOCATABLE :: iwork
#endif
WAUTELET Philippe
committed
INTEGER :: idx, ji,jj
INTEGER :: kcdf_id
INTEGER :: extent, ndims
INTEGER :: src
INTEGER :: level
INTEGER(KIND=LFI_INT) :: iresp,ilu,ileng,ipos
CHARACTER(LEN=4) :: suffix
WAUTELET Philippe
committed
INTEGER,DIMENSION(:),ALLOCATABLE :: itab
REAL(KIND=8),DIMENSION(:),ALLOCATABLE :: xtab
CHARACTER, DIMENSION(:), ALLOCATABLE :: ytab
REAL(KIND=8), DIMENSION(:,:,:), ALLOCATABLE :: xtab3d, xtab3d2
INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: itab3d, itab3d2
WAUTELET Philippe
committed
IF (infiles%files(1)%format == LFI_FORMAT) ilu = infiles%files(1)%lun_id
IF (present(current_level)) THEN
write(suffix,'(I4.4)') current_level
level = current_level
ElSE
suffix=''
level = 1
END IF
#if LOWMEM
ALLOCATE(iwork(kbuflen))
#endif
ALLOCATE(itab(kbuflen))
ALLOCATE(xtab(kbuflen))
WAUTELET Philippe
committed
idx = 1
IF (.NOT.tpreclist(ji)%tbw) CYCLE
WAUTELET Philippe
committed
IF (ASSOCIATED(tpreclist(ji)%dim)) THEN
extent = tpreclist(ji)%dim%len
ndims = tpreclist(ji)%dim%ndims
ndims = 0
idims(:) = 1
if(ndims>0) idims(1) = ptdimx%len
if(ndims>1) idims(2) = ptdimy%len
if(ndims>2) idims(3) = ptdimz%len
if(ndims>3) then
PRINT *,'Too many dimensions'
STOP
endif
SELECT CASE(tpreclist(ji)%TYPE)
CASE (INT,BOOL)
IF (infiles%files(1)%format == LFI_FORMAT) THEN
IF (.NOT.tpreclist(ji)%calc) THEN
CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
itab(1:extent) = iwork(3+iwork(2):3+iwork(2)+extent-1)
ELSE
src=tpreclist(ji)%src(1)
CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
itab(1:extent) = iwork(3+iwork(2):3+iwork(2)+extent-1)
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
itab(1:extent) = itab(1:extent) + iwork(3+iwork(2):3+iwork(2)+extent-1)
IF (.NOT.tpreclist(ji)%calc) THEN
itab(1:extent) = lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):)
ELSE
src=tpreclist(ji)%src(1)
itab(1:extent) = lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
itab(1:extent) = xtab(1:extent) + lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):)
jj=jj+1
END DO
END IF
!TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
SELECT CASE(ndims)
CASE (0)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1))
CASE (1)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab(1:extent),count=(/extent/))
CASE (2)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len/)), &
start = (/1,1,level/) )
CASE (3)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(itab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
CASE DEFAULT
print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
END SELECT
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
ALLOCATE( itab3d(idims(1),idims(2),idims(3)) )
IF (.NOT.tpreclist(ji)%calc) THEN
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,itab3d)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
ELSE
ALLOCATE( itab3d2(idims(1),idims(2),idims(3)) )
src=tpreclist(ji)%src(1)
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,itab3d2)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
itab3d(:,:,:) = itab3d(:,:,:) + itab3d2(:,:,:)
jj=jj+1
END DO
DEALLOCATE(itab3d2)
END IF
!TODO: not clean, should be done only if merging z-levels
IF (ndims == 2) THEN
start = (/1,1,level/)
ELSE
start = (/1,1,1/)
ENDIF
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,itab3d,start=start)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
DEALLOCATE(itab3d)
END IF
IF (infiles%files(1)%format == LFI_FORMAT) THEN
IF (.NOT.tpreclist(ji)%calc) THEN
CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
ELSE
src=tpreclist(ji)%src(1)
CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
jj=jj+1
END DO
ENDIF
IF (.NOT.tpreclist(ji)%calc) THEN
xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
ELSE
src=tpreclist(ji)%src(1)
xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
jj=jj+1
END DO
END IF
!TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
SELECT CASE(ndims)
CASE (0)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
CASE (1)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
CASE (2)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
start = (/1,1,level/) )
CASE (3)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
CASE DEFAULT
print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
END SELECT
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
ALLOCATE( xtab3d(idims(1),idims(2),idims(3)) )
IF (.NOT.tpreclist(ji)%calc) THEN
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,xtab3d)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
ELSE
ALLOCATE( xtab3d2(idims(1),idims(2),idims(3)) )
src=tpreclist(ji)%src(1)
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(src)%id_in,xtab3d2)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
xtab3d(:,:,:) = xtab3d(:,:,:) + xtab3d2(:,:,:)
jj=jj+1
END DO
DEALLOCATE(xtab3d2)
END IF
!TODO: not clean, should be done only if merging z-levels
IF (ndims == 2) THEN
start = (/1,1,level/)
ELSE
start = (/1,1,1/)
ENDIF
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab3d,start=start)
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
DEALLOCATE(xtab3d)
END IF
IF (infiles%files(1)%format == LFI_FORMAT) THEN
#if LOWMEM
CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
#endif
ALLOCATE(ytab(extent))
DO jj=1,extent
#if LOWMEM
ich = iwork(2+iwork(2)+jj)
#else
ich = lfiart(ji)%iwtab(2+lfiart(ji)%iwtab(2)+jj)
#endif
ytab(jj) = CHAR(ich)
END DO
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
status = NF90_GET_VAR(infiles%files(1)%lun_id,tpreclist(ji)%id_in,ytab,count=(/extent/))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,ytab,count=(/extent/))
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
END IF
IF (infiles%files(1)%format == LFI_FORMAT) THEN
IF (.NOT.tpreclist(ji)%calc) THEN
CALL LFINFO(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(ji)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
ELSE
src=tpreclist(ji)%src(1)
CALL LFINFO(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),ileng,ipos)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
CALL LFILEC(iresp,ilu,trim(tpreclist(src)%name)//trim(suffix),iwork,ileng)
xtab(1:extent) = xtab(1:extent) + TRANSFER(iwork(3+iwork(2):),(/ 0.0_8 /))
jj=jj+1
END DO
ENDIF
IF (.NOT.tpreclist(ji)%calc) THEN
xtab(1:extent) = TRANSFER(lfiart(ji)%iwtab(3+lfiart(ji)%iwtab(2):),(/ 0.0_8 /))
ELSE
src=tpreclist(ji)%src(1)
xtab(1:extent) = TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
jj = 2
DO WHILE (tpreclist(ji)%src(jj)>0 .AND. jj.LE.MAXRAW)
src=tpreclist(ji)%src(jj)
xtab(1:extent) = xtab(1:extent) + TRANSFER(lfiart(src)%iwtab(3+lfiart(src)%iwtab(2):),(/ 0.0_8 /))
jj=jj+1
END DO
END IF
!TODO: works in all cases??? (X, Y, Z dimensions assumed to be ptdimx,y or z)
SELECT CASE(ndims)
CASE (0)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1))
CASE (1)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,xtab(1:extent),count=(/extent/))
CASE (2)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len/)), &
start = (/1,1,level/) )
CASE (3)
status = NF90_PUT_VAR(kcdf_id,tpreclist(ji)%id_out,reshape(xtab,(/ptdimx%len,ptdimy%len,ptdimz%len/)))
CASE DEFAULT
print *,'Error: arrays with ',tpreclist(ji)%dim%ndims,' dimensions are not supported'
END SELECT
IF (status /= NF90_NOERR) CALL HANDLE_ERR(status,__LINE__)
ELSE IF (infiles%files(1)%format == NETCDF_FORMAT) THEN
print *,'Error: unknown datatype'
STOP
END IF
WAUTELET Philippe
committed
if (options(OPTSPLIT)%set) idx = idx + 1
END DO
DEALLOCATE(itab,xtab)
#if LOWMEM
DEALLOCATE(iwork)
#endif