directaccess.f90

      MODULE directaccess
      USE safe_open_mod
      USE stel_kinds, ONLY: dp
 
!     DIRECT ACCESS FILE HANDLING ROUTINES
      INTEGER :: rec_length, data_size, block_size, num_rows
      INTEGER :: iunit_da=0, blocks_per_row, recs_per_block
      INTEGER :: irec_pos, byte_size_rec, byte_size_dp
      CHARACTER(LEN=256) :: filename
 
      CONTAINS 
 
      SUBROUTINE opendafile(datasize, blksize, blocksperrow,            &
     &                      filename_in, iunit, iflag)
      INTEGER, INTENT(in) :: datasize, blksize, blocksperrow,           &
     &                       iflag
      INTEGER, INTENT(inout) :: iunit 
      CHARACTER*(*), INTENT(in) :: filename_in
      INTEGER, PARAMETER :: CreateNew=0, openexisting=1, scratch=2
      INTEGER :: ierr
      REAL(dp) :: dummy
      CHARACTER(LEN=10) :: Status
   
!GET EFFECTIVE "byte_size" OF UNFORMATTED dp VARIABLE (machine dependent!)
!ON Compaq, rec length size is in units of 4 byte chunks, so byte_size_rec = 4
!ON lf95, byte_size_rec = 1
      INQUIRE(iolength=byte_size_rec) dummy
      byte_size_dp = kind(dummy)
   
      data_size = datasize
      rec_length = byte_size_rec*datasize
      block_size = byte_size_dp*blksize
      blocks_per_row = blocksperrow
      filename = filename_in
!skip this distance to next block, if datasize different from blksize
      recs_per_block  = max(1,blksize/datasize)
      irec_pos = 0
 
!  create disk file for doing direct access i/o.
 
      IF (iflag .eq. createnew) THEN
         status = "replace"
      ELSE IF (iflag .eq. openexisting) THEN
         status = "old"
      ELSE
         status = "scratch"
      END IF
 
 
      CALL safe_open(iunit, ierr, filename, status, 'unformatted',   &
     &     rec_length, 'DIRECT')
 
      iunit_da = iunit
 
      IF (ierr .ne. 0) THEN
         WRITE (6, '(a7,i4)') 'Status code: ', status, ' Error stat: ', ierr
         stop 'Error creating Direct Access file!'
      END IF
 
      END SUBROUTINE opendafile
 
      
      SUBROUTINE changedafileparams(datasize, blksize,                  &
     &                              blocksperrow, new_filename, nrows)
      INTEGER, INTENT(in) :: datasize, blksize, blocksperrow, nrows
      INTEGER :: ierr, new_data_size, new_block_size, new_rec_length,   &
     &           new_blocks_per_row, new_recs_per_block, inew_da
      INTEGER :: i, j, k, boffset, recloc, nsplit, isplit, new_row_offset
      CHARACTER*(*)  :: new_filename
      REAL(dp), ALLOCATABLE :: DataItem(:)
 
    
!store new parameters
      new_data_size = datasize
      new_rec_length = byte_size_rec*datasize
      new_block_size = byte_size_dp*blksize
      new_blocks_per_row = blocksperrow
!skip this distance to next block, if rec_length different from block_size
      new_recs_per_block  = max(1,blksize/datasize)
      inew_da      = iunit_da+1
 
      IF ((rec_length.eq.new_rec_length) .and.                          &
     &    (block_size.eq.new_block_size)) RETURN
 
 
!  create disk file for doing direct access i/o.
      ierr = index(filename,'.',back=.true.)
 
      filename = new_filename
 
      CALL safe_open(inew_da, ierr, filename, 'replace', 'unformatted', &
     &     new_rec_length, 'DIRECT')
      IF (ierr .ne. 0) stop 'Error opening existing Direct Access file!'
 
      ALLOCATE (dataitem(new_data_size), stat=ierr)
      IF (ierr .ne. 0) stop 'Allocation error in ChangeDAFileParams'
 
!Load up block by reading a column at a time
      IF (irec_pos .eq. 0) THEN
         DO i = 1, nrows
            DO j = 1, blocks_per_row
               boffset = 1
               DO k = 1, recs_per_block
                  CALL readdaitem1(dataitem(boffset), i, j, k)
                  boffset = boffset + recs_per_block
               END DO
               !write full block in new da file
               recloc = 1 + new_recs_per_block*((j-1) + new_blocks_per_row*(i-1))
               WRITE (inew_da, rec=recloc, iostat=ierr) dataitem
            END DO
         END DO
 
      ELSE 
         num_rows = nrows
 
!        Data comes out as follows for each COLUMN (m,n,ntype):
!        L2,D1,U0 L5,D4,U3  .... L(3*n1+2),D(3*n1+1),U(3*n1)
!        L3,D2,U1 L6,D5,U4       L(3*n2+3),D(3*n2+2),U(3*n2+1)
!        L4,D3,U2 L7,D6,U5       L(3*n3+4),D(3*n3+3),U(3*n3+2)
!
!        Therefore we must split the i=1,nrows loop into 3
         isplit = 0
         DO nsplit = 1, 3
         DO i = nsplit, nrows, 3
            isplit = isplit+1         !sequential index of records in original file
            DO j = 1, blocks_per_row    
!        j=1=>L  2=>D  3=>U block
               boffset = 1
               new_row_offset = 2-j+i
               IF (new_row_offset.lt.1 .or. new_row_offset.gt.nrows) cycle
               DO k = 1, recs_per_block
!        Read one column (k) at a time
                  CALL readdaitem_seq(dataitem(boffset), isplit, j, k)
                  boffset = boffset + recs_per_block
               END DO
               !write full block to new da file
               recloc = 1 + new_recs_per_block*((j-1) + new_blocks_per_row*(new_row_offset-1))
               WRITE (inew_da, rec=recloc, iostat=ierr) dataitem
            END DO
         END DO
         END DO
 
         IF (isplit .ne. nrows) stop 'isplit != nrows'
 
      END IF
 
!     Close old scratch file when finished writing out new file
      CLOSE (iunit_da, status='DELETE')
 
      iunit_da = inew_da
      data_size = new_data_size
      rec_length = new_rec_length
      block_size = new_block_size
      blocks_per_row = new_blocks_per_row
      recs_per_block = new_recs_per_block
      irec_pos = 0
 
      READ(inew_da, rec=2) dataitem
!     WRITE (133, '(1p,4e14.4)') DataItem
 
      DEALLOCATE (dataitem)
 
      END SUBROUTINE changedafileparams
 
      
      SUBROUTINE closedafile
 
      IF (iunit_da .gt. 0) THEN
         CLOSE (iunit_da)
         iunit_da = 0
      END IF
 
      END SUBROUTINE closedafile
 
 
      SUBROUTINE deletedafile (filename)
      CHARACTER*(*) :: filename
      INTEGER :: ierr, rec_length=1
      
      IF (iunit_da .eq. 0) THEN
         iunit_da=100
         CALL safe_open(iunit_da, ierr, filename, 'replace', 'unformatted', &
     &                  rec_length, 'DIRECT')
         IF (ierr .ne. 0) THEN
            print *,'Unable to open existing ScratchFile'
            RETURN
         END IF
      END IF
 
      CLOSE (iunit_da, status='DELETE')
      iunit_da = 0
 
      END SUBROUTINE deletedafile
 
     
      SUBROUTINE writedaitem_ra(DataItem, BlockRowIndex, ColIndex, IndexInBlock)
      REAL(dp), INTENT(in) :: DataItem(data_size)
      INTEGER, INTENT(in)  :: BlockRowIndex, ColIndex, IndexInBlock
      INTEGER :: recloc, ierr
      INTEGER :: StartIndex
 
      IF (colindex > blocks_per_row) stop 'ColIndex > Block_Per_Row in WriteDAItem'
      IF (indexinblock > recs_per_block)  stop 'IndexInBloc > skip_size in WriteDAItem'
         
 
      startindex = indexinblock
      IF (recs_per_block .eq. 1) startindex = 1
      recloc = startindex + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
 
      WRITE (iunit_da, rec=recloc, iostat=ierr) dataitem
      IF (ierr .ne. 0) THEN
         WRITE (6,*) 'Ierr = ', ierr, ' in WriteDAItem'
         stop
      END IF
 
      END SUBROUTINE writedaitem_ra
 
      
      SUBROUTINE writedaitem_seq(DataItem)
      REAL(dp), INTENT(in) :: DataItem(data_size)
      INTEGER :: ierr
 
      !Perform "sequential" write of direct access elements (buffering helps)
 
      irec_pos = irec_pos+1
 
      WRITE (iunit_da, rec=irec_pos, iostat=ierr) dataitem
      IF (ierr .ne. 0) THEN
         WRITE (6,*) 'Ierr = ', ierr, ' in WriteDAItem'
         stop
      END IF
 
      END SUBROUTINE writedaitem_seq
 
 
      SUBROUTINE readdaitem1(DataItem, BlockRowIndex, ColIndex, StartIndex)
      REAL(dp), INTENT(out) :: DataItem(data_size)
      INTEGER, INTENT(in)  :: BlockRowIndex, ColIndex, StartIndex
      INTEGER :: recloc, ierr
!      INTEGER :: StartIndex=1
 
      recloc = startindex + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
      READ (iunit_da, rec=recloc, iostat=ierr) dataitem
      IF (ierr .ne. 0) THEN
         WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
         stop
      END IF
 
      END SUBROUTINE readdaitem1
 
 
      SUBROUTINE readdaitem2(DataItem, BlockRowIndex, ColIndex)
      REAL(dp), INTENT(out) :: DataItem(data_size)
      INTEGER, INTENT(in)  :: BlockRowIndex, ColIndex
      INTEGER :: recloc, ierr
!      INTEGER :: StartIndex=1
 
      recloc = 1 + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
      READ (iunit_da, rec=recloc, iostat=ierr) dataitem
      IF (ierr .ne. 0) THEN
         WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
         stop
      END IF
 
      END SUBROUTINE readdaitem2
 
 
      SUBROUTINE readdaitem_seq(DataItem, BlockRowIndex, BlockTypeIndex, ColIndex)
      REAL(dp), INTENT(out) :: DataItem(data_size)
      INTEGER, INTENT(in)  :: BlockRowIndex, ColIndex, BlockTypeIndex
      INTEGER :: recloc, ierr
 
      recloc = blocktypeindex + blocks_per_row*(blockrowindex-1+num_rows*(colindex-1))
      READ (iunit_da, rec=recloc, iostat=ierr) dataitem
      IF (ierr .ne. 0) THEN
         WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
         stop
      END IF
 
      END SUBROUTINE readdaitem_seq
 
 
      END MODULE directaccess