!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