V3FIT
directaccess.f90
1  MODULE directaccess
2  USE safe_open_mod
3  USE stel_kinds, ONLY: dp
4 
5 ! DIRECT ACCESS FILE HANDLING ROUTINES
6  INTEGER :: rec_length, data_size, block_size, num_rows
7  INTEGER :: iunit_da=0, blocks_per_row, recs_per_block
8  INTEGER :: irec_pos, byte_size_rec, byte_size_dp
9  CHARACTER(LEN=256) :: filename
10 
11  CONTAINS
12 
13  SUBROUTINE opendafile(datasize, blksize, blocksperrow, &
14  & filename_in, iunit, iflag)
15  INTEGER, INTENT(in) :: datasize, blksize, blocksperrow, &
16  & iflag
17  INTEGER, INTENT(inout) :: iunit
18  CHARACTER*(*), INTENT(in) :: filename_in
19  INTEGER, PARAMETER :: CreateNew=0, openexisting=1, scratch=2
20  INTEGER :: ierr
21  REAL(dp) :: dummy
22  CHARACTER(LEN=10) :: Status
23 
24 !GET EFFECTIVE "byte_size" OF UNFORMATTED dp VARIABLE (machine dependent!)
25 !ON Compaq, rec length size is in units of 4 byte chunks, so byte_size_rec = 4
26 !ON lf95, byte_size_rec = 1
27  INQUIRE(iolength=byte_size_rec) dummy
28  byte_size_dp = kind(dummy)
29 
30  data_size = datasize
31  rec_length = byte_size_rec*datasize
32  block_size = byte_size_dp*blksize
33  blocks_per_row = blocksperrow
34  filename = filename_in
35 !skip this distance to next block, if datasize different from blksize
36  recs_per_block = max(1,blksize/datasize)
37  irec_pos = 0
38 
39 ! create disk file for doing direct access i/o.
40 
41  IF (iflag .eq. createnew) THEN
42  status = "replace"
43  ELSE IF (iflag .eq. openexisting) THEN
44  status = "old"
45  ELSE
46  status = "scratch"
47  END IF
48 
49 
50  CALL safe_open(iunit, ierr, filename, status, 'unformatted', &
51  & rec_length, 'DIRECT')
52 
53  iunit_da = iunit
54 
55  IF (ierr .ne. 0) THEN
56  WRITE (6, '(a7,i4)') 'Status code: ', status, ' Error stat: ', ierr
57  stop 'Error creating Direct Access file!'
58  END IF
59 
60  END SUBROUTINE opendafile
61 
62 
63  SUBROUTINE changedafileparams(datasize, blksize, &
64  & blocksperrow, new_filename, nrows)
65  INTEGER, INTENT(in) :: datasize, blksize, blocksperrow, nrows
66  INTEGER :: ierr, new_data_size, new_block_size, new_rec_length, &
67  & new_blocks_per_row, new_recs_per_block, inew_da
68  INTEGER :: i, j, k, boffset, recloc, nsplit, isplit, new_row_offset
69  CHARACTER*(*) :: new_filename
70  REAL(dp), ALLOCATABLE :: DataItem(:)
71 
72 
73 !store new parameters
74  new_data_size = datasize
75  new_rec_length = byte_size_rec*datasize
76  new_block_size = byte_size_dp*blksize
77  new_blocks_per_row = blocksperrow
78 !skip this distance to next block, if rec_length different from block_size
79  new_recs_per_block = max(1,blksize/datasize)
80  inew_da = iunit_da+1
81 
82  IF ((rec_length.eq.new_rec_length) .and. &
83  & (block_size.eq.new_block_size)) RETURN
84 
85 
86 ! create disk file for doing direct access i/o.
87  ierr = index(filename,'.',back=.true.)
88 
89  filename = new_filename
90 
91  CALL safe_open(inew_da, ierr, filename, 'replace', 'unformatted', &
92  & new_rec_length, 'DIRECT')
93  IF (ierr .ne. 0) stop 'Error opening existing Direct Access file!'
94 
95  ALLOCATE (dataitem(new_data_size), stat=ierr)
96  IF (ierr .ne. 0) stop 'Allocation error in ChangeDAFileParams'
97 
98 !Load up block by reading a column at a time
99  IF (irec_pos .eq. 0) THEN
100  DO i = 1, nrows
101  DO j = 1, blocks_per_row
102  boffset = 1
103  DO k = 1, recs_per_block
104  CALL readdaitem1(dataitem(boffset), i, j, k)
105  boffset = boffset + recs_per_block
106  END DO
107  !write full block in new da file
108  recloc = 1 + new_recs_per_block*((j-1) + new_blocks_per_row*(i-1))
109  WRITE (inew_da, rec=recloc, iostat=ierr) dataitem
110  END DO
111  END DO
112 
113  ELSE
114  num_rows = nrows
115 
116 ! Data comes out as follows for each COLUMN (m,n,ntype):
117 ! L2,D1,U0 L5,D4,U3 .... L(3*n1+2),D(3*n1+1),U(3*n1)
118 ! L3,D2,U1 L6,D5,U4 L(3*n2+3),D(3*n2+2),U(3*n2+1)
119 ! L4,D3,U2 L7,D6,U5 L(3*n3+4),D(3*n3+3),U(3*n3+2)
120 !
121 ! Therefore we must split the i=1,nrows loop into 3
122  isplit = 0
123  DO nsplit = 1, 3
124  DO i = nsplit, nrows, 3
125  isplit = isplit+1 !sequential index of records in original file
126  DO j = 1, blocks_per_row
127 ! j=1=>L 2=>D 3=>U block
128  boffset = 1
129  new_row_offset = 2-j+i
130  IF (new_row_offset.lt.1 .or. new_row_offset.gt.nrows) cycle
131  DO k = 1, recs_per_block
132 ! Read one column (k) at a time
133  CALL readdaitem_seq(dataitem(boffset), isplit, j, k)
134  boffset = boffset + recs_per_block
135  END DO
136  !write full block to new da file
137  recloc = 1 + new_recs_per_block*((j-1) + new_blocks_per_row*(new_row_offset-1))
138  WRITE (inew_da, rec=recloc, iostat=ierr) dataitem
139  END DO
140  END DO
141  END DO
142 
143  IF (isplit .ne. nrows) stop 'isplit != nrows'
144 
145  END IF
146 
147 ! Close old scratch file when finished writing out new file
148  CLOSE (iunit_da, status='DELETE')
149 
150  iunit_da = inew_da
151  data_size = new_data_size
152  rec_length = new_rec_length
153  block_size = new_block_size
154  blocks_per_row = new_blocks_per_row
155  recs_per_block = new_recs_per_block
156  irec_pos = 0
157 
158  READ(inew_da, rec=2) dataitem
159 ! WRITE (133, '(1p,4e14.4)') DataItem
160 
161  DEALLOCATE (dataitem)
162 
163  END SUBROUTINE changedafileparams
164 
165 
166  SUBROUTINE closedafile
167 
168  IF (iunit_da .gt. 0) THEN
169  CLOSE (iunit_da)
170  iunit_da = 0
171  END IF
172 
173  END SUBROUTINE closedafile
174 
175 
176  SUBROUTINE deletedafile (filename)
177  CHARACTER*(*) :: filename
178  INTEGER :: ierr, rec_length=1
179 
180  IF (iunit_da .eq. 0) THEN
181  iunit_da=100
182  CALL safe_open(iunit_da, ierr, filename, 'replace', 'unformatted', &
183  & rec_length, 'DIRECT')
184  IF (ierr .ne. 0) THEN
185  print *,'Unable to open existing ScratchFile'
186  RETURN
187  END IF
188  END IF
189 
190  CLOSE (iunit_da, status='DELETE')
191  iunit_da = 0
192 
193  END SUBROUTINE deletedafile
194 
195 
196  SUBROUTINE writedaitem_ra(DataItem, BlockRowIndex, ColIndex, IndexInBlock)
197  REAL(dp), INTENT(in) :: DataItem(data_size)
198  INTEGER, INTENT(in) :: BlockRowIndex, ColIndex, IndexInBlock
199  INTEGER :: recloc, ierr
200  INTEGER :: StartIndex
201 
202  IF (colindex > blocks_per_row) stop 'ColIndex > Block_Per_Row in WriteDAItem'
203  IF (indexinblock > recs_per_block) stop 'IndexInBloc > skip_size in WriteDAItem'
204 
205 
206  startindex = indexinblock
207  IF (recs_per_block .eq. 1) startindex = 1
208  recloc = startindex + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
209 
210  WRITE (iunit_da, rec=recloc, iostat=ierr) dataitem
211  IF (ierr .ne. 0) THEN
212  WRITE (6,*) 'Ierr = ', ierr, ' in WriteDAItem'
213  stop
214  END IF
215 
216  END SUBROUTINE writedaitem_ra
217 
218 
219  SUBROUTINE writedaitem_seq(DataItem)
220  REAL(dp), INTENT(in) :: DataItem(data_size)
221  INTEGER :: ierr
222 
223  !Perform "sequential" write of direct access elements (buffering helps)
224 
225  irec_pos = irec_pos+1
226 
227  WRITE (iunit_da, rec=irec_pos, iostat=ierr) dataitem
228  IF (ierr .ne. 0) THEN
229  WRITE (6,*) 'Ierr = ', ierr, ' in WriteDAItem'
230  stop
231  END IF
232 
233  END SUBROUTINE writedaitem_seq
234 
235 
236  SUBROUTINE readdaitem1(DataItem, BlockRowIndex, ColIndex, StartIndex)
237  REAL(dp), INTENT(out) :: DataItem(data_size)
238  INTEGER, INTENT(in) :: BlockRowIndex, ColIndex, StartIndex
239  INTEGER :: recloc, ierr
240 ! INTEGER :: StartIndex=1
241 
242  recloc = startindex + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
243  READ (iunit_da, rec=recloc, iostat=ierr) dataitem
244  IF (ierr .ne. 0) THEN
245  WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
246  stop
247  END IF
248 
249  END SUBROUTINE readdaitem1
250 
251 
252  SUBROUTINE readdaitem2(DataItem, BlockRowIndex, ColIndex)
253  REAL(dp), INTENT(out) :: DataItem(data_size)
254  INTEGER, INTENT(in) :: BlockRowIndex, ColIndex
255  INTEGER :: recloc, ierr
256 ! INTEGER :: StartIndex=1
257 
258  recloc = 1 + recs_per_block*((colindex-1) + blocks_per_row*(blockrowindex-1))
259  READ (iunit_da, rec=recloc, iostat=ierr) dataitem
260  IF (ierr .ne. 0) THEN
261  WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
262  stop
263  END IF
264 
265  END SUBROUTINE readdaitem2
266 
267 
268  SUBROUTINE readdaitem_seq(DataItem, BlockRowIndex, BlockTypeIndex, ColIndex)
269  REAL(dp), INTENT(out) :: DataItem(data_size)
270  INTEGER, INTENT(in) :: BlockRowIndex, ColIndex, BlockTypeIndex
271  INTEGER :: recloc, ierr
272 
273  recloc = blocktypeindex + blocks_per_row*(blockrowindex-1+num_rows*(colindex-1))
274  READ (iunit_da, rec=recloc, iostat=ierr) dataitem
275  IF (ierr .ne. 0) THEN
276  WRITE (6,*) 'Ierr = ', ierr, ' in ReadDAItem'
277  stop
278  END IF
279 
280  END SUBROUTINE readdaitem_seq
281 
282 
283  END MODULE directaccess