!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
SUBROUTINE join_hdf (fileheader,nxsm,nysm,nz,nxlg,nylg,buf_r,buf_rsm, & 3,27
buf_i,buf_ism,buf_i16,buf_i16sm,buf_r1,buf_r2,sstat)
IMPLICIT NONE
INCLUDE 'mp.inc'
INCLUDE 'hdf.f90' ! HDF4 library include file
CHARACTER (LEN=*) :: fileheader
INTEGER :: nxsm,nysm,nz,nxlg,nylg
INTEGER :: sstat ! split status: sstat=1 if error encountered
!-----------------------------------------------------------------------
!
! Misc. local variables
!
!-----------------------------------------------------------------------
CHARACTER (LEN=128) :: filename
INTEGER :: fi, fj, i, j, k
INTEGER :: nxin, nyin, nzin
REAL :: buf_r(nxlg,nylg,nz), buf_rsm(nxsm,nysm,nz)
REAL :: buf_r1(nxsm+nysm+nz), buf_r2(nxlg+nylg+nz), amax, amin, scale
INTEGER :: buf_i(nxlg,nylg,nz), buf_ism(nxsm,nysm,nz)
INTEGER (KIND=selected_int_kind(4)) :: &
buf_i16(nxlg,nylg,nz), buf_i16sm(nxsm,nysm,nz)
INTEGER :: ierr
INTEGER :: sd_id,sd_id1,sd_id2,ndata,nattr,istat,aindex,dindex
INTEGER :: sds_id,sds_id1,sds_id2,rank,dims(6),dtype,ndattr,adtype
INTEGER :: sds_index
CHARACTER (LEN=128) :: name, aname
CHARACTER (LEN=1024) :: char_attr
CHARACTER (LEN=1), ALLOCATABLE :: lchar_attr(:)
INTEGER :: nvalues
INTEGER :: size(3),start(3),stride(3),strideout(3),sizeout(3)
INTEGER :: x_off,y_off,i0,j0
INTEGER :: comp_code, comp_prm(1)
INTEGER :: itmp
!-----------------------------------------------------------------------
!
! Functions
!
!-----------------------------------------------------------------------
INTEGER :: sfcreate, sfrdata, sfrnatt, sfscompress, sfselect, &
sfsnatt, sfwdata, sffattr, sfendacc, sfn2index
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
! nxlg = (nxsm-3)*nproc_x+3
! nylg = (nysm-3)*nproc_y+3
stride = 1
start = 0
sstat = 0
!-----------------------------------------------------------------------
!
! Open file 0101, read in its attributes and check the
! dimensions in the file against the dimensions passed in.
! Open the joined file as well.
!
!-----------------------------------------------------------------------
CALL hdfopen
(trim(fileheader)//'_0101',1,sd_id1)
IF (sd_id1 < 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR opening ", &
trim(fileheader)//'_0101'," for reading."
sstat = 1
RETURN
END IF
CALL hdfopen(trim(fileheader),2,sd_id2)
IF (sd_id2 < 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR creating HDF4 file: ", &
trim(fileheader)
sstat = 1
GO TO 600
END IF
CALL hdfinfo(sd_id1,ndata,nattr,istat)
CALL hdfrdi(sd_id1,"nx",nxin,istat)
CALL hdfrdi(sd_id1,"ny",nyin,istat)
IF (nz > 1) THEN
CALL hdfrdi(sd_id1,"nz",nzin,istat)
ELSE
nzin = nz ! the file is 2-D so it won't have nz
ENDIF
IF ((nxin /= nxsm).OR.(nyin /= nysm).OR.(nzin /= nz)) THEN
WRITE (*,*) "ERROR: mismatch in sizes."
WRITE (*,*) "nxin,nyin,nzin: ",nxin,nyin,nzin
WRITE (*,*) "nxsm,nysm,nz: ",nxsm,nysm,nz
sstat = 1
RETURN
END IF
!-----------------------------------------------------------------------
!
! Read/write header info.
!
!-----------------------------------------------------------------------
DO aindex = 0,nattr-1
CALL hdfainfo(sd_id1,aindex,name,dtype,nvalues,istat)
IF (dtype == dfnt_char8) THEN
IF (nvalues > 1024) THEN
ALLOCATE (lchar_attr(nvalues))
CALL hdfrdc(sd_id1,nvalues,name,lchar_attr,istat)
CALL hdfwrtc(sd_id2,nvalues,name,lchar_attr,istat)
DEALLOCATE(lchar_attr)
ELSE
CALL hdfrdc(sd_id1,nvalues,name,char_attr,istat)
CALL hdfwrtc(sd_id2,nvalues,name,char_attr,istat)
ENDIF
ELSE IF (dtype == dfnt_float32) THEN
istat = sfrnatt(sd_id1, aindex, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading attribute ",trim(name)
sstat = 1
GOTO 600
ENDIF
istat = sfsnatt(sd_id2, trim(name), dfnt_float32, nvalues, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR, writing attribute ",trim(name)
sstat = 1
GOTO 600
ENDIF
ELSE IF (dtype == dfnt_int32) THEN
IF (trim(name) == 'nx') THEN
CALL hdfwrti(sd_id2, 'nx', nxlg, istat)
ELSE IF (trim(name) == 'ny') THEN
CALL hdfwrti(sd_id2, 'ny', nylg, istat)
ELSE
istat = sfrnatt(sd_id1, aindex, buf_i)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading attribute ",trim(name)
sstat = 1
GOTO 600
ENDIF
istat = sfsnatt(sd_id2, trim(name), dfnt_int32, nvalues, buf_i)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR, writing attribute ",trim(name)
sstat = 1
GOTO 600
ENDIF
ENDIF
ELSE
WRITE (6,*) "JOIN_HDF: ERROR, unknown data type for ", &
"attribute ", trim(name)
sstat = 1
GOTO 600
ENDIF
END DO
!-----------------------------------------------------------------------
!
! Read/write each data set.
!
!-----------------------------------------------------------------------
DO dindex = 0,ndata-1
sds_id1 = sfselect(sd_id1,dindex)
! data set
CALL hdfdinfo(sds_id1,name,rank,dims,dtype,ndattr,istat)
CALL hdfrdi(sds_id1,"hdf_comp_prm",comp_prm,istat)
CALL hdfrdi(sds_id1,"hdf_comp_code",comp_code,istat)
DO fj=1,nproc_y
DO fi=1,nproc_x
WRITE (filename, '(a,a,2i2.2)') &
trim(fileheader),'_',fi,fj
CALL hdfopen(filename,1,sd_id)
IF (sd_id < 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR opening HDF4 file: ", &
trim(filename)
sstat = 1
GO TO 600
END IF
sds_index = sfn2index(sd_id, trim(name))
IF (sds_index == -1) THEN
WRITE (6,*) "JOIN_HDF: ERROR, variable ", &
trim(name)," not found in file",trim(filename),"."
sstat = 1
GO TO 600
END IF
sds_id = sfselect(sd_id, sds_index)
x_off = (fi-1) * (nxsm-3)
y_off = (fj-1) * (nysm-3)
i0 = min(2,fi)
j0 = min(2,fj)
size(1) = nxsm
size(2) = nysm
size(3) = dims(3)
sizeout(1) = nxlg
sizeout(2) = nylg
sizeout(3) = dims(3)
IF (rank == 1) THEN
IF (dtype /= dfnt_float32) THEN
WRITE (6,*) "JOIN_HDF: ERROR, unsuppored data type for 1-d ", &
" variable ",trim(name)
sstat = 1
GOTO 600
ENDIF
IF (trim(name) == 'x') THEN ! x
IF (fj == 1) THEN
size(1) = nxsm
istat = sfrdata(sds_id, start, stride, size, buf_r1)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
DO i=i0,nxsm
buf_r2(i+x_off) = buf_r1(i)
END DO
ENDIF
ELSE IF (trim(name) == 'y') THEN ! y
IF (fi == 1) THEN
size(1) = nysm
istat = sfrdata(sds_id, start, stride, size, buf_r1)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
DO j=j0,nysm
buf_r2(j+y_off) = buf_r1(j)
END DO
ENDIF
ELSE IF (trim(name) == 'z') THEN ! z
IF (fi == 1 .and. fj == 1) THEN
size(1) = nz
istat = sfrdata(sds_id, start, stride, dims, buf_r1)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
ENDIF
ELSE
IF (fi == 1 .and. fj == 1) THEN
istat = sfrdata(sds_id, start, stride, dims, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
ENDIF
ENDIF
ELSE
IF (dtype == dfnt_float32) THEN
istat = sfrdata(sds_id, start, stride, size, buf_rsm)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
DO k=1,nz
DO j=j0,nysm
DO i=i0,nxsm
buf_r(i+x_off,j+y_off,k) = buf_rsm(i,j,k)
END DO
END DO
END DO
ELSE IF (dtype == dfnt_int32) THEN
istat = sfrdata(sds_id, start, stride, size, buf_ism)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
DO k=1,nz
DO j=j0,nysm
DO i=i0,nxsm
buf_i(i+x_off,j+y_off,k) = buf_ism(i,j,k)
END DO
END DO
END DO
ELSE IF (dtype == dfnt_int16) THEN
IF (rank == 2) THEN
aindex = sffattr(sds_id, "max")
istat = sfrnatt(sds_id, aindex, amax)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading max for ", &
trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
aindex = sffattr(sds_id, "min")
istat = sfrnatt(sds_id, aindex, amin)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading min for ", &
trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
istat = sfrdata(sds_id, start, stride, size, buf_i16sm)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
scale = (amax - amin) / 65534.0
DO j=j0,nysm
DO i=i0,nxsm
buf_r(i+x_off,j+y_off,1) = &
scale * (buf_i16sm(i,j,1) + 32767) + amin
END DO
END DO
ELSE IF (rank == 3) THEN
aindex = sffattr(sds_id, "max")
istat = sfrnatt(sds_id, aindex, buf_r2)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading max for ", &
trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
aindex = sffattr(sds_id, "min")
istat = sfrnatt(sds_id, aindex, buf_r1)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading min for ", &
trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
istat = sfrdata(sds_id, start, stride, size, buf_i16sm)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading ",trim(name),", exiting"
sstat = 1
GOTO 600
ENDIF
DO k=1,nz
scale = (buf_r2(k) - buf_r1(k)) / 65534.0
DO j=j0,nysm
DO i=i0,nxsm
buf_r(i+x_off,j+y_off,k) = &
scale * (buf_i16sm(i,j,k) + 32767) + buf_r1(k)
END DO
END DO
END DO
ELSE
WRITE (6,*) "JOIN_HDF: ERROR, unsupported rank,",rank, &
" for 16 bit remapped data"
GOTO 600
ENDIF
ELSE
WRITE (6,*) "JOIN_HDF: ERROR, unknown data type for ", &
"attribute ", trim(aname), " of sds ",trim(name)
sstat = 1
GOTO 600
ENDIF
ENDIF
CALL hdfclose(sd_id,istat)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR on close of file ",trim(filename)
sstat = 1
GOTO 600
ENDIF
END DO ! fi
END DO ! fj
IF (rank == 1) THEN
IF (trim(name) == 'x') THEN ! x
sizeout(1) = nxlg
sds_id2 = sfcreate(sd_id2, trim(name), dfnt_float32, 1, sizeout)
istat = sfwdata(sds_id2, start, stride, sizeout, buf_r2)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ELSE IF (trim(name) == 'y') THEN ! y
sizeout(1) = nylg
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_float32,1,sizeout)
istat = sfwdata(sds_id2, start, stride, sizeout, buf_r2)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ELSE IF (trim(name) == 'z') THEN ! z
sizeout(1) = nz
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_float32,1,sizeout)
istat = sfwdata(sds_id2, start, stride, sizeout, buf_r1)
! note, buf_r1 not buf_r2
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ELSE
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_float32,1,dims)
istat = sfwdata(sds_id2, start, stride, dims, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ENDIF
ELSE
IF (dtype == dfnt_float32) THEN
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_float32,rank,sizeout)
IF (comp_prm(1) > 0) THEN
istat = sfscompress(sds_id2, comp_code, comp_prm)
ENDIF
istat = sfwdata(sds_id2, start, stride, sizeout, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ELSE IF (dtype == dfnt_int32) THEN
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_int32,rank,sizeout)
IF (comp_prm(1) > 0) THEN
istat = sfscompress(sds_id2, comp_code, comp_prm)
ENDIF
istat = sfwdata(sds_id2, start, stride, sizeout, buf_i)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ELSE IF (dtype == dfnt_int16) THEN
sds_id2 = sfcreate(sd_id2,trim(name),dfnt_int16,rank,sizeout)
IF (comp_prm(1) > 0) THEN
istat = sfscompress(sds_id2, comp_code, comp_prm)
ENDIF
IF (rank == 2) THEN
CALL a3dmax0lcl(buf_r,1,nxlg,1,nxlg,1,nylg,1,nylg,1,1,1,1,amax,amin)
IF (ABS(amax-amin) < 1.0E-10) amax = 1.1 * amin + 1.
scale = 65534.0 / (amax - amin)
DO j=1,nylg
DO i=1,nxlg
itmp = nint(scale * (buf_r(i,j,1) - amin)) - 32767
buf_i16(i,j,1) = itmp
END DO
END DO
istat = sfsnatt(sds_id2, 'max', dfnt_float32, 1, amax)
istat = sfsnatt(sds_id2, 'min', dfnt_float32, 1, amin)
ELSE IF (rank == 3) THEN
DO k=1,nz
CALL a3dmax0lcl(buf_r(1,1,k),1,nxlg,1,nxlg,1,nylg,1,nylg, &
1,nz,1,1,buf_r2(k),buf_r1(k))
IF (ABS(buf_r2(k)-buf_r1(k)) < 1.0E-10) &
buf_r2(k) = 1.1 * buf_r1(k) + 1.
scale = 65534.0 / (buf_r2(k) - buf_r1(k))
DO j=1,nylg
DO i=1,nxlg
itmp = nint(scale * (buf_r(i,j,k) - buf_r1(k))) - 32767
buf_i16(i,j,k) = itmp
END DO
END DO
END DO
istat = sfsnatt(sds_id2, 'max', dfnt_float32, nz, buf_r2)
istat = sfsnatt(sds_id2, 'min', dfnt_float32, nz, buf_r1)
ENDIF
istat = sfwdata(sds_id2, start, stride, sizeout, buf_i16)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing ",trim(name), &
" to file ",trim(fileheader)," , exiting"
sstat = 1
GOTO 600
ENDIF
ENDIF
ENDIF
! data set attributes
DO aindex = 0,ndattr-1
CALL hdfainfo(sds_id1,aindex,aname,adtype,nvalues,istat)
IF (dtype /= dfnt_int16 .or. &
(trim(aname) /= "max" .and. trim(aname) /= "min")) THEN
! Don't write out max & min for 16 bit integer remapped reals
IF (adtype == dfnt_char8) THEN
IF (nvalues > 1024) THEN
ALLOCATE (lchar_attr(nvalues))
CALL hdfrdc(sds_id1,nvalues,aname,lchar_attr,istat)
CALL hdfwrtc(sds_id2,nvalues,aname,lchar_attr,istat)
DEALLOCATE(lchar_attr)
ELSE
CALL hdfrdc(sds_id1,nvalues,aname,char_attr,istat)
CALL hdfwrtc(sds_id2,nvalues,aname,char_attr,istat)
ENDIF
ELSE IF (adtype == dfnt_float32) THEN
istat = sfrnatt(sds_id1, aindex, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading attribute ",trim(aname)
sstat = 1
GOTO 600
ENDIF
istat = sfsnatt(sds_id2, trim(aname), dfnt_float32, nvalues, buf_r)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR, writing attribute ",trim(aname)
sstat = 1
GOTO 600
ENDIF
ELSE IF (adtype == dfnt_int32) THEN
istat = sfrnatt(sds_id1, aindex, buf_i)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR reading attribute ",trim(aname)
sstat = 1
GOTO 600
ENDIF
istat = sfsnatt(sds_id2, trim(aname), dfnt_int32, nvalues, buf_i)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR, writing attribute ",trim(aname)
sstat = 1
GOTO 600
ENDIF
ELSE
WRITE (6,*) "JOIN_HDF: ERROR, unknown data type for ", &
"attribute ", trim(aname)
sstat = 1
GOTO 600
ENDIF
ENDIF
END DO
istat = sfendacc(sds_id2)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR writing variable ",trim(name)
END IF
END DO ! dindex
!-----------------------------------------------------------------------
!
! Close I/O and deallocate.
!
!----------------------------------------------------------------------
600 CONTINUE
CALL hdfclose(sd_id1,istat)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR on close of file ",trim(fileheader)//'_0101'
sstat = 1
ENDIF
CALL hdfclose(sd_id2,istat)
IF (istat /= 0) THEN
WRITE (6,*) "JOIN_HDF: ERROR on close of file ",trim(fileheader)
sstat = 1
ENDIF
RETURN
END SUBROUTINE join_hdf