00001 MODULE TestDatasets
00002
00003 USE Xmdf
00004
00005 CHARACTER(LEN=*), PARAMETER :: DATASETS_LOCATION = 'Datasets'
00006 CHARACTER(LEN=*), PARAMETER :: SCALAR_A_LOCATION = 'Scalars/ScalarA'
00007 CHARACTER(LEN=*), PARAMETER :: SCALAR_B_LOCATION = 'Scalars/ScalarB'
00008 CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_LOCATION = 'Vectors/Vector2D_A'
00009 CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_LOCATION = 'Vectors/Vector2D_B'
00010
00011 CONTAINS
00012
00013 ! --------------------------------------------------------------------------
00014 ! FUNCTION tdReadDatasets
00015 ! PURPOSE Read a dataset group from an XMDF file and output information to
00016 ! to a text file
00017 ! NOTES
00018 ! --------------------------------------------------------------------------
00019 RECURSIVE SUBROUTINE TD_READ_DATASETS (a_xGroupId, a_FileUnit, error)
00020 INTEGER, INTENT(IN) :: a_xGroupId
00021 INTEGER, INTENT(IN) :: a_FileUnit
00022 INTEGER, INTENT(OUT) :: error
00023 INTEGER nPaths, nMaxPathLength, j
00024 CHARACTER, ALLOCATABLE, DIMENSION(:) :: Paths
00025 CHARACTER(LEN=500) IndividualPath
00026 INTEGER nStatus, i
00027 INTEGER xScalarId, xVectorId, xMultiId
00028 INTEGER nMultiDatasets
00029
00030 xScalarId = NONE
00031 xVectorId = NONE
00032
00033 nMultiDatasets = 0
00034 nPaths = 0
00035 nMaxPathLength = 0
00036
00037 ! Look for scalar datasets
00038 call XF_GET_SCALAR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
00039 if (nStatus >= 0 .AND. nPaths > 0) then
00040 allocate(Paths(nPaths*nMaxPathLength))
00041 call XF_GET_SCALAR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, &
00042 error)
00043 endif
00044 if (nStatus < 0) then
00045 error = -1
00046 return
00047 endif
00048
00049 ! Output number and paths to scalar datasets
00050 WRITE(a_FileUnit,*) 'Number of Scalars ', nPaths
00051 do i=2, nPaths
00052 IndividualPath = ''
00053 do j=1, nMaxPathLength-1
00054 IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
00055 enddo
00056 WRITE(a_FileUnit,*) 'Reading scalar: ', IndividualPath(1:nMaxPathLength-1)
00057 call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
00058 xScalarId, nStatus)
00059 if (nStatus < 0) then
00060 error = -1
00061 return
00062 endif
00063
00064 call TDI_READ_SCALAR(xScalarId, a_FileUnit, nStatus)
00065 call XF_CLOSE_GROUP(xScalarId, error)
00066 if (nStatus < 0) then
00067 WRITE(*,*) 'Error reading scalar dataset.'
00068 error = -1
00069 return
00070 endif
00071 enddo
00072
00073 if (allocated(Paths)) deallocate(Paths)
00074 ! Look for vector datasets
00075 call XF_GET_VECTOR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
00076 if (nStatus >= 0 .AND. nPaths > 0) then
00077 allocate(Paths(nPaths*nMaxPathLength))
00078 call XF_GET_VECTOR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, error)
00079 endif
00080 if (nStatus < 0) then
00081 error = -1
00082 return
00083 endif
00084
00085 ! Output number and paths to scalar datasets
00086 WRITE(a_FileUnit,*) 'Number of Vectors ', nPaths
00087 do i=2, nPaths
00088 do j=1, nMaxPathLength-1
00089 IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
00090 enddo
00091 WRITE(a_FileUnit,*) 'Reading Vector: ', &
00092 IndividualPath(1:nMaxPathLength-1)
00093 call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
00094 xVectorId, nStatus)
00095 if (nStatus < 0) then
00096 error = -1
00097 return
00098 endif
00099 call TDI_READ_VECTOR(xVectorId, a_FileUnit, nStatus)
00100 call XF_CLOSE_GROUP(xVectorId, error)
00101 if (nStatus < 0) then
00102 WRITE(*,*) 'Error reading vector dataset.'
00103 error = -1
00104 return
00105 endif
00106 enddo
00107
00108 if (allocated(Paths)) deallocate(Paths)
00109
00110 ! find multidataset folders
00111 call XF_GET_GRP_PTHS_SZ_MLT_DSETS(a_xGroupId, nMultiDatasets, &
00112 nMaxPathLength, nStatus)
00113 if (nStatus >= 0 .AND. nMultiDatasets > 0) then
00114 allocate(Paths(nMultiDatasets*nMaxPathLength))
00115 call XF_GET_ALL_GRP_PATHS_MLT_DSETS(a_xGroupId, nMultiDatasets, &
00116 nMaxPathLength, Paths, error)
00117 if (nStatus < 0) then
00118 error = -1
00119 return
00120 endif
00121
00122 ! Output number and paths to multidatasets
00123 WRITE(a_FileUnit,*) 'Number of Multidatasets ', nMultiDatasets
00124 do i=2, nMultiDatasets
00125 IndividualPath = ''
00126 do j=1, nMaxPathLength-1
00127 IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
00128 enddo
00129 WRITE(a_FileUnit,*) 'Reading multidataset: ', &
00130 IndividualPath(1:nMaxPathLength-1)
00131 call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
00132 xMultiId, nStatus)
00133 if (nStatus < 0) then
00134 error = -1
00135 return
00136 endif
00137
00138 call TD_READ_DATASETS(xMultiId, a_FileUnit, nStatus)
00139 call XF_CLOSE_GROUP(xMultiId, error)
00140 if (nStatus < 0) then
00141 WRITE(*,*) 'Error reading multidatasets.'
00142 error = -1
00143 return
00144 endif
00145 enddo
00146 endif
00147 if (allocated(Paths)) deallocate(Paths)
00148
00149 error = 1
00150 return
00151
00152 END SUBROUTINE
00153 !tdReadDatasets
00154 ! --------------------------------------------------------------------------
00155 ! FUNCTION tdReadActivityScalarAIndex
00156 ! PURPOSE Read all timestep values for a particular index
00157 ! NOTES
00158 ! --------------------------------------------------------------------------
00159 SUBROUTINE TD_READ_ACTIVITY_SCALAR_A_INDEX(a_Filename, a_Index, error)
00160 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00161 INTEGER, INTENT(IN) :: a_Index
00162 INTEGER, INTENT(OUT) :: error
00163 INTEGER status
00164 INTEGER xFileId, xDsetsId, xScalarAId
00165 INTEGER nTimesteps, i
00166 INTEGER, ALLOCATABLE :: bActive(:)
00167
00168 xFileId = NONE
00169 xDsetsId = NONE
00170 xScalarAId = NONE
00171
00172 ! open the file
00173 call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
00174 if (status < 0) then
00175 error = -1
00176 return
00177 endif
00178
00179 ! open the dataset group
00180 call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00181 if (status >= 0) then
00182 call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
00183 endif
00184 if (status < 0) then
00185 error = status
00186 return
00187 endif
00188
00189 ! Find out the number of timesteps in the file
00190 CALL XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
00191 if (status < 0) then
00192 error = status
00193 return
00194 endif
00195
00196 if (nTimesteps < 1) then
00197 error = -1
00198 return
00199 endif
00200
00201 ! Read the values for the index
00202 allocate(bActive(nTimesteps))
00203 call XF_READ_ACTIVE_VALS_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
00204 bActive, status)
00205 ! output the data
00206 WRITE(*,*) ''
00207 WRITE(*,*) 'Reading activity for scalar A slice at index: ', a_Index
00208 do i=1, nTimesteps
00209 WRITE(*,*) bActive(i), ' '
00210 enddo
00211
00212 deallocate(bActive)
00213
00214 error = status
00215 return
00216
00217 END SUBROUTINE
00218 ! tdReadActivityScalarAIndex
00219
00220 ! --------------------------------------------------------------------------
00221 ! FUNCTION tdReadScalarAIndex
00222 ! PURPOSE Read all timestep values for a particular index
00223 ! NOTES
00224 ! --------------------------------------------------------------------------
00225 SUBROUTINE TD_READ_SCALAR_A_INDEX (a_Filename, a_Index, error)
00226 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00227 INTEGER, INTENT(IN) :: a_Index
00228 INTEGER, INTENT(OUT) :: error
00229 INTEGER status
00230 INTEGER xFileId, xDsetsId, xScalarAId
00231 INTEGER nTimesteps, i
00232 REAL, ALLOCATABLE :: fValues(:)
00233
00234 xFileId = NONE
00235 xDsetsId = NONE
00236 xScalarAId = NONE
00237
00238 ! open the file
00239 call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
00240 if (status < 0) then
00241 error = -1
00242 return
00243 endif
00244
00245 ! open the dataset group
00246 call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00247 if (status >= 0) then
00248 call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
00249 endif
00250 if (status < 0) then
00251 error = status
00252 return
00253 endif
00254
00255 ! Find out the number of timesteps in the file
00256 call XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
00257 if (status < 0) then
00258 error = status
00259 return
00260 endif
00261
00262 if (nTimesteps < 1) then
00263 error = -1
00264 return
00265 endif
00266
00267 ! Read the values for the index
00268 allocate (fValues(nTimesteps))
00269 call XF_READ_SCALAR_VALUES_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
00270 fValues, status)
00271
00272 ! output the data
00273 WRITE(*,*) ''
00274 WRITE(*,*) 'Reading scalar A slice at index: ', a_Index
00275 do i=1, nTimesteps
00276 WRITE(*,*) fValues(i), ' '
00277 enddo
00278
00279 deallocate(fValues)
00280
00281 error = status
00282 return
00283
00284 END SUBROUTINE
00285 ! tdReadScalarAtIndex
00286
00287 ! --------------------------------------------------------------------------
00288 ! FUNCTION tdWriteScalarA
00289 ! PURPOSE Write scalar Dataset to an HDF5 File
00290 ! NOTES This tests dynamic data sets, and activity
00291 ! This dataset is dynamic concentrations (mg/L) with output times
00292 ! in minutes.
00293 ! Dataset is for a mesh and so nActive is the number of elements
00294 ! which is not the same as the nValues which would be number of nodes
00295 ! reads/writes a reference time in julian days
00296 ! --------------------------------------------------------------------------
00297 SUBROUTINE TD_WRITE_SCALAR_A (a_Filename, a_Compression, error)
00298 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00299 INTEGER, INTENT(IN) :: a_Compression
00300 INTEGER, INTENT(OUT) :: error
00301 INTEGER xFileId, xDsetsId, xScalarAId, xCoordId
00302 INTEGER nValues, nTimes, nActive
00303 REAL(DOUBLE) dTime, dJulianReftime
00304 INTEGER iTimestep, iActive, iHpgnZone
00305 REAL fValues(10) ! nValues
00306 INTEGER*1 bActivity(10) ! activity
00307 INTEGER i, status
00308
00309 ! initialize the data
00310 nValues = 10
00311 nTimes = 3
00312 nActive = 8
00313 dTime = 0.0
00314
00315 ! 5th item in data set is always inactive, others active
00316 do iActive = 1, nActive
00317 bActivity(iActive) = 1
00318 enddo
00319 bActivity(6) = 0
00320
00321
00322 ! create the file
00323 call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
00324 if (error .LT. 0) then
00325 ! close the file
00326 call XF_CLOSE_FILE(xFileId, error)
00327 return
00328 endif
00329
00330 ! create the group where we will put all the datasets
00331 call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00332 if (status < 0) then
00333 call XF_CLOSE_FILE(xFileId, error)
00334 error = -1
00335 return
00336 endif
00337
00338 ! Create the scalar A dataset group
00339 call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_A_LOCATION, 'mg/L', &
00340 TS_HOURS, a_Compression, xScalarAId, status)
00341 if (status .LT. 0) then
00342 ! close the dataset
00343 call XF_CLOSE_GROUP(xScalarAId, error)
00344 call XF_CLOSE_GROUP(xDsetsId, error)
00345 call XF_CLOSE_FILE(xFileId, error)
00346 error = status
00347 return
00348 endif
00349
00350 ! Add in a reftime. This is a julian day for:
00351 ! noon July 1, 2003
00352 dJulianReftime = 2452822.0;
00353 call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
00354 if (status < 0) then
00355 call XF_CLOSE_GROUP(xScalarAId, error)
00356 call XF_CLOSE_GROUP(xDsetsId, error)
00357 call XF_CLOSE_FILE(xFileId, error)
00358 endif
00359
00360 ! Loop through timesteps adding them to the file
00361 do iTimestep = 1, nTimes
00362 ! We will have an 0.5 hour timestep
00363 dTime = iTimestep * 0.5
00364
00365 fValues(1) = dTime
00366 do i = 2, nValues
00367 fValues(i) = fValues(i-1)*2.5
00368 end do
00369
00370 ! write the dataset array values
00371 call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, error)
00372 if (error .GE. 0) then
00373 ! write activity array
00374 call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, error)
00375 end if
00376 enddo
00377
00378 ! Write Coordinate file - for ScalarA, we will set the coordinate system
00379 ! to be Geographic HPGN, with HPGN settings written to the file.
00380 call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
00381 if (status < 0) then
00382 call XF_CLOSE_GROUP(xScalarAId, error)
00383 call XF_CLOSE_GROUP(xDsetsId, error)
00384 call XF_CLOSE_FILE(xFileId, error)
00385 error = -1
00386 return
00387 endif
00388
00389 ! set HPGN Zone for test
00390 iHpgnZone = 29 ! Utah
00391 ! Write Coordinate Information to file
00392 call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
00393 call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
00394 call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
00395 call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)
00396
00397 ! write additional information
00398 call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)
00399
00400 call XF_CLOSE_GROUP(xCoordId, error)
00401 xCoordId = 0;
00402
00403 ! close the dataset
00404 call XF_CLOSE_GROUP(xScalarAId, error)
00405 call XF_CLOSE_GROUP(xDsetsId, error)
00406 call XF_CLOSE_FILE(xFileId, error)
00407
00408 return
00409 END SUBROUTINE
00410 ! tdWriteScalarA
00411
00412 ! --------------------------------------------------------------------------
00413 ! FUNCTION TD_WRITE_SCALAR_B
00414 ! PURPOSE Write scalar Dataset to an HDF5 File
00415 ! NOTES This tests dynamic data sets, and activity
00416 ! This dataset is dynamic concentrations (mg/L) with output times
00417 ! in minutes.
00418 ! Dataset is for a mesh and so nActive is the number of elements
00419 ! which is not the same as the nValues which would be number of nodes
00420 ! reads/writes a reference time in julian days
00421 ! --------------------------------------------------------------------------
00422 SUBROUTINE TD_WRITE_SCALAR_B (a_Filename, a_Compression, a_Overwrite, error)
00423 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00424 INTEGER, INTENT(IN) :: a_Compression
00425 LOGICAL, INTENT(IN) :: a_Overwrite
00426 INTEGER, INTENT(OUT) :: error
00427 INTEGER xFileId, xDsetsId, xScalarBId, xCoordId
00428 INTEGER nValues, nTimes, nActive
00429 REAL(DOUBLE) dTime, dJulianReftime
00430 INTEGER iTimestep, iActive
00431 REAL fValues(10) ! nValues
00432 INTEGER*1 bActivity(10) ! activity
00433 INTEGER i, status
00434
00435 ! initialize the data
00436 nValues = 10
00437 nTimes = 3
00438 nActive = 8
00439 dTime = 0.0
00440 i = 0
00441
00442 ! 5th item in data set is always inactive, others active
00443 do iActive = 1, nActive
00444 bActivity(iActive) = 1
00445 enddo
00446 bActivity(6) = 0
00447
00448 if (a_Overwrite) then
00449 ! open the already-existing file
00450 call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
00451 if (status < 0) then
00452 error = -1
00453 return
00454 endif
00455 ! open the group where we have all the datasets
00456 call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00457 if (status < 0) then
00458 call XF_CLOSE_FILE(xFileId, error)
00459 error = -1
00460 return
00461 endif
00462 else
00463 ! create the file
00464 call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
00465 if (error .LT. 0) then
00466 ! close the file
00467 call XF_CLOSE_FILE(xFileId, error)
00468 return
00469 endif
00470
00471 ! create the group where we will put all the datasets
00472 call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00473 if (status < 0) then
00474 call XF_CLOSE_FILE(xFileId, error)
00475 error = -1
00476 return
00477 endif
00478 endif
00479
00480 ! Create/Overwrite the scalar B dataset group
00481 call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_B_LOCATION, 'mg/L', &
00482 TS_HOURS, a_Compression, xScalarBId, status)
00483 if (status < 0) then
00484 ! close the dataset
00485 call XF_CLOSE_GROUP(xScalarBId, error)
00486 call XF_CLOSE_GROUP(xDsetsId, error)
00487 call XF_CLOSE_FILE(xFileId, error)
00488 error = status
00489 return
00490 endif
00491
00492 ! Add in a reftime. This is a julian day for:
00493 ! noon July 1, 2003
00494 dJulianReftime = 2452822.0;
00495 call XF_WRITE_REFTIME(xScalarBId, dJulianReftime, status)
00496 if (status < 0) then
00497 call XF_CLOSE_GROUP(xScalarBId, error)
00498 call XF_CLOSE_GROUP(xDsetsId, error)
00499 call XF_CLOSE_FILE(xFileId, error)
00500 endif
00501
00502 if (.NOT. a_Overwrite) then
00503 ! Loop through timesteps adding them to the file
00504 do iTimestep = 1, nTimes
00505 ! We will have an 0.5 hour timestep
00506 dTime = iTimestep * 0.5
00507
00508 fValues(1) = dTime
00509 do i = 2, nValues
00510 fValues(i) = fValues(i-1)*2.5
00511 end do
00512
00513 ! write the dataset array values
00514 call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
00515 if (error .GE. 0) then
00516 ! write activity array
00517 call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
00518 end if
00519 if (error < 0) then
00520 call XF_CLOSE_GROUP(xScalarBId, error)
00521 call XF_CLOSE_GROUP(xDsetsId, error)
00522 call XF_CLOSE_FILE(xFileId, error)
00523 endif
00524 enddo
00525 else
00526 ! Loop through timesteps adding them to the file
00527 do iTimestep = 1, nTimes
00528 ! We will have an 1.5 hour timestep
00529 dTime = iTimestep * 1.5
00530
00531 fValues(1) = dTime
00532 do i = 2, nValues
00533 fValues(i) = fValues(i-1)*1.5
00534 end do
00535
00536 ! write the dataset array values
00537 call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
00538 if (error .GE. 0) then
00539 ! write activity array
00540 call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
00541 end if
00542 if (error < 0) then
00543 call XF_CLOSE_GROUP(xScalarBId, error)
00544 call XF_CLOSE_GROUP(xDsetsId, error)
00545 call XF_CLOSE_FILE(xFileId, error)
00546 endif
00547 enddo
00548 endif
00549
00550 if (.NOT. a_Overwrite) then
00551 ! Write Coordinate file - for ScalarB, we will set the coordinate system
00552 ! to be UTM, with UTM Zone settings written to the file.
00553 call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
00554 if (status < 0) then
00555 call XF_CLOSE_GROUP(xScalarBId, error)
00556 call XF_CLOSE_GROUP(xDsetsId, error)
00557 call XF_CLOSE_FILE(xFileId, error)
00558 error = -1
00559 return
00560 endif
00561
00562 ! Write Coord Info to file
00563 call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
00564 call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
00565
00566 call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
00567 call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)
00568
00569 ! write additional information - we'll use the max value for this test
00570 call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)
00571
00572 call XF_CLOSE_GROUP(xCoordId, error)
00573 xCoordId = 0
00574 endif
00575
00576 ! close the dataset
00577 call XF_CLOSE_GROUP(xScalarBId, error)
00578 call XF_CLOSE_GROUP(xDsetsId, error)
00579 call XF_CLOSE_FILE(xFileId, error)
00580
00581 error = 1
00582 return
00583 END SUBROUTINE
00584 ! tdWriteScalarB
00585 !------------------------------------------------------------------------------
00586 ! FUNCTION TD_WRITE_COORDS_TO_MULTI
00587 ! PURPOSE Write coordinate system to a multidataset file
00588 ! NOTES
00589 !------------------------------------------------------------------------------
00590 SUBROUTINE TD_WRITE_COORDS_TO_MULTI (a_xFileId, error)
00591 INTEGER, INTENT(IN) :: a_xFileId
00592 INTEGER, INTENT(OUT) :: error
00593 INTEGER xCoordId
00594 INTEGER status
00595
00596 ! Write Coordinate file - for Multidatasets, we will set the coordinate system
00597 ! to be UTM, with UTM Zone settings written to the file.
00598 call XF_CREATE_COORDINATE_GROUP(a_xFileId, xCoordId, status)
00599 if (status < 0) then
00600 error = status
00601 return
00602 endif
00603
00604 ! Write Coord Info to file
00605 call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
00606 call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
00607
00608 call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
00609 call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)
00610
00611 ! write additional information - we'll use the max value for this test
00612 call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)
00613
00614 call XF_CLOSE_GROUP(xCoordId, error)
00615 xCoordId = 0
00616
00617 return
00618 END SUBROUTINE
00619
00620 ! --------------------------------------------------------------------------
00621 ! FUNCTION tdWriteScalarAToMulti
00622 ! PURPOSE Write scalar Dataset to an HDF5 File
00623 ! NOTES This tests dynamic data sets, and activity
00624 ! This dataset is dynamic concentrations (mg/L) with output times
00625 ! in minutes.
00626 ! Dataset is for a mesh and so nActive is the number of elements
00627 ! which is not the same as the nValues which would be number of nodes
00628 ! reads/writes a reference time in julian days
00629 ! --------------------------------------------------------------------------
00630 SUBROUTINE TD_WRITE_SCALAR_A_TO_MULTI (a_GroupID, status)
00631 ! CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00632 ! INTEGER, INTENT(IN) :: a_Compression
00633 ! INTEGER, INTENT(OUT) :: error
00634 INTEGER xFileId, xDsetsId, xScalarAId
00635 INTEGER a_GroupID
00636 INTEGER nValues, nTimes, nActive
00637 REAL(DOUBLE) dTime, dJulianReftime
00638 INTEGER iTimestep, iActive
00639 REAL fValues(10) ! nValues
00640 INTEGER*1 bActivity(10) ! activity
00641 INTEGER i, status
00642
00643 ! initialize the data
00644 nValues = 10
00645 nTimes = 3
00646 nActive = 8
00647 dTime = 0.0
00648
00649 ! 5th item in data set is always inactive, others active
00650 do iActive = 1, nActive
00651 bActivity(iActive) = 1
00652 enddo
00653 bActivity(6) = 0
00654
00655 ! Create the scalar A dataset group
00656 call XF_CREATE_SCALAR_DATASET(a_GroupID, SCALAR_A_LOCATION, 'mg/L', &
00657 TS_HOURS, NONE, xScalarAId, status)
00658 if (status .LT. 0) then
00659 ! close the dataset
00660 call XF_CLOSE_GROUP(xScalarAId, status)
00661 call XF_CLOSE_GROUP(xDsetsId, status)
00662 call XF_CLOSE_FILE(xFileId, status)
00663 return
00664 endif
00665
00666 ! Add in a reftime. This is a julian day for:
00667 ! noon July 1, 2003
00668 dJulianReftime = 2452822.0;
00669 call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
00670 if (status < 0) then
00671 call XF_CLOSE_GROUP(xScalarAId, status)
00672 call XF_CLOSE_GROUP(xDsetsId, status)
00673 call XF_CLOSE_FILE(xFileId, status)
00674 endif
00675
00676 ! Loop through timesteps adding them to the file
00677 do iTimestep = 1, nTimes
00678 ! We will have an 0.5 hour timestep
00679 dTime = iTimestep * 0.5
00680
00681 fValues(1) = dTime
00682 do i = 2, nValues
00683 fValues(i) = fValues(i-1)*2.5
00684 end do
00685
00686 ! write the dataset array values
00687 call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, status)
00688 if (status .GE. 0) then
00689 ! write activity array
00690 call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, status)
00691 end if
00692 enddo
00693
00694 ! close the dataset
00695 call XF_CLOSE_GROUP(xScalarAId, status)
00696 !call XF_CLOSE_GROUP(a_GroupID, status)
00697 !call XF_CLOSE_FILE(a_FileID, status)
00698
00699 return
00700 END SUBROUTINE
00701 ! tdWriteScalarAToMulti
00702 ! --------------------------------------------------------------------------
00703 ! FUNCTION tdReadVector2DAIndex
00704 ! PURPOSE Read all timestep values for a particular index
00705 ! NOTES
00706 ! --------------------------------------------------------------------------
00707 SUBROUTINE TD_READ_VECTOR2D_A_INDEX (a_Filename, a_Index, error)
00708 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00709 INTEGER, INTENT(IN) :: a_Index
00710 INTEGER, INTENT(OUT) :: error
00711 INTEGER status
00712 INTEGER xFileId, xDsetsId, xVector2DA
00713 INTEGER nTimesteps, i
00714 REAL, ALLOCATABLE :: fValues(:)
00715
00716 xFileId = NONE
00717 xDsetsId = NONE
00718 xVector2DA = NONE
00719
00720 ! open the file
00721 call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
00722 if (status < 0) then
00723 error = -1
00724 return
00725 endif
00726
00727 ! open the dataset group
00728 call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00729 if (status >= 0) then
00730 call XF_OPEN_GROUP(xDsetsId, VECTOR2D_A_LOCATION, xVector2DA, status)
00731 endif
00732 if (status < 0) then
00733 error = status
00734 return
00735 endif
00736
00737 ! Find out the number of timesteps in the file
00738 call XF_GET_DATASET_NUM_TIMES(xVector2DA, nTimesteps, status)
00739 if (status < 0) then
00740 error = status
00741 return
00742 endif
00743
00744 if (nTimesteps < 1) then
00745 error = -1
00746 return
00747 endif
00748
00749 ! Read the values for the index
00750 allocate(fValues(nTimesteps*2))
00751 call XF_READ_VECTOR_VALUES_AT_INDEX(xVector2DA, a_Index, 1, nTimesteps, 2, &
00752 fValues, status)
00753
00754 ! output the data
00755 WRITE(*,*) ''
00756 WRITE(*,*) 'Reading vector 2D A slice at index: ', a_Index
00757 do i=1, nTimesteps
00758 WRITE(*,*) fValues(i*2-1), ' ', fValues(i*2)
00759 enddo
00760 WRITE(*,*) ''
00761
00762 deallocate(fValues)
00763
00764 error = status
00765 return
00766
00767 END SUBROUTINE
00768 !tdReadVector2DAIndex
00769
00770 ! --------------------------------------------------------------------------
00771 ! FUNCTION tdWriteVector2D_A
00772 ! PURPOSE Write scalar Dataset to an HDF5 File
00773 ! NOTES This tests dynamic data sets, and activity
00774 ! This dataset is dynamic concentrations (mg/L) with output times
00775 ! in minutes.
00776 ! Dataset is for a mesh and so nActive is the number of elements
00777 ! which is not the same as the nValues which would be number of nodes
00778 ! reads/writes a reference time in julian days
00779 ! --------------------------------------------------------------------------
00780 SUBROUTINE TD_WRITE_VECTOR2D_A (a_Filename, a_Compression, error)
00781 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00782 INTEGER, INTENT(IN) :: a_Compression
00783 INTEGER, INTENT(OUT) :: error
00784 INTEGER xFileId, xDsetsId, xVector2D_A, xCoordId
00785 INTEGER nValues, nTimes, nComponents, nActive
00786 REAL(DOUBLE) dTime
00787 INTEGER iTimestep, iActive
00788 REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
00789 INTEGER*1 bActivity(100) ! activity
00790 INTEGER i, j, status
00791 INTEGER iHpgnZone
00792
00793 ! initialize the data
00794 nComponents = 2
00795 nValues = 100
00796 nTimes = 6
00797 nActive = 75
00798 dTime = 0.0
00799
00800 ! 5th item in data set is always inactive, others active
00801 bActivity(1) = 0
00802 do iActive = 2, nActive
00803 if (mod(iActive-1, 3) == 0) then
00804 bActivity(iActive) = 0
00805 else
00806 bActivity(iActive) = 1
00807 endif
00808 enddo
00809
00810 ! create the file
00811 call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
00812 if (error .LT. 0) then
00813 ! close the dataset
00814 call XF_CLOSE_FILE(xFileId, error)
00815 return
00816 endif
00817
00818 ! create the group where we will put all the datasets
00819 call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00820 if (status < 0) then
00821 call XF_CLOSE_FILE(xFileId, error)
00822 error = -1
00823 return
00824 endif
00825
00826 ! Create the vector dataset group
00827 call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_A_LOCATION, 'ft/s', &
00828 TS_SECONDS, a_Compression, xVector2D_A, status)
00829 if (status .LT. 0) then
00830 ! close the dataset
00831 call XF_CLOSE_GROUP(xVector2D_A, error)
00832 call XF_CLOSE_GROUP(xDsetsId, error)
00833 call XF_CLOSE_FILE(xFileId, error)
00834 error = status
00835 return
00836 endif
00837
00838 ! Loop through timesteps adding them to the file
00839 do iTimestep = 1, nTimes
00840 ! We will have an 0.5 hour timestep
00841 dTime = iTimestep * 0.5
00842
00843 do i = 1, nValues
00844 do j = 1, nComponents
00845 fValues(j,i) = ((i-1)*nComponents + j)*dTime
00846 end do
00847 end do
00848
00849 ! write the dataset array values
00850 call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
00851 fValues, error)
00852 if (error .GE. 0) then
00853 ! write activity array
00854 call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, error)
00855 end if
00856 enddo
00857
00858 ! Write Coordinate file - for Vector2D_A, we will set the coordinate system
00859 ! to be Geographic HPGN, with HPGN settings written to the file.
00860 call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
00861 if (status < 0) then
00862 call XF_CLOSE_GROUP(xVector2D_A, error)
00863 call XF_CLOSE_GROUP(xDsetsId, error)
00864 call XF_CLOSE_FILE(xFileId, error)
00865 error = -1
00866 return
00867 endif
00868
00869 ! set HPGN info for test
00870 iHpgnZone = 29 ! Utah
00871
00872 call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
00873 call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
00874 call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
00875 call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)
00876
00877 ! write additional information
00878 call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)
00879
00880 call XF_CLOSE_GROUP(xCoordId, error)
00881 xCoordId = 0
00882
00883 ! close the dataset
00884 call XF_CLOSE_GROUP(xVector2D_A, error)
00885 call XF_CLOSE_GROUP(xDsetsId, error)
00886 call XF_CLOSE_FILE(xFileId, error)
00887
00888 return
00889 END SUBROUTINE
00890 ! tdWriteVector2D_A
00891
00892 ! --------------------------------------------------------------------------
00893 ! FUNCTION TD_WRITE_VECTOR2D_B
00894 ! PURPOSE Write scalar Dataset to an HDF5 File
00895 ! NOTES This tests dynamic data sets, and activity
00896 ! This dataset is dynamic concentrations (mg/L) with output times
00897 ! in minutes.
00898 ! Dataset is for a mesh and so nActive is the number of elements
00899 ! which is not the same as the nValues which would be number of nodes
00900 ! reads/writes a reference time in julian days
00901 ! --------------------------------------------------------------------------
00902 SUBROUTINE TD_WRITE_VECTOR2D_B (a_Filename, a_Compression, a_Overwrite, error)
00903 CHARACTER(LEN=*), INTENT(IN) :: a_Filename
00904 INTEGER, INTENT(IN) :: a_Compression
00905 LOGICAL, INTENT(IN) :: a_Overwrite
00906 INTEGER, INTENT(OUT) :: error
00907 INTEGER xFileId, xDsetsId, xVector2D_B, xCoordId
00908 INTEGER nValues, nTimes, nComponents, nActive
00909 REAL(DOUBLE) dTime
00910 INTEGER iTimestep, iActive
00911 REAL, DIMENSION(2, 100) :: fValues
00912 INTEGER*1 bActivity(100)
00913 INTEGER i, j, status
00914
00915 ! initialize the data
00916 nComponents = 2
00917 nValues = 100
00918 nTimes = 6
00919 nActive = 75
00920 dTime = 0.0
00921
00922 ! 5th item in data set is always inactive, others active
00923 bActivity(1) = 0
00924 do iActive = 2, nActive
00925 if (mod(iActive-1, 3) == 0) then
00926 bActivity(iActive) = 0
00927 else
00928 bActivity(iActive) = 1
00929 endif
00930 enddo
00931
00932 if (a_Overwrite) then
00933 ! open the already-existing file
00934 call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
00935 if (status < 0) then
00936 error = -1
00937 return
00938 endif
00939 ! open the group where we have all the datasets
00940 call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00941 if (status < 0) then
00942 call XF_CLOSE_FILE(xFileId, error)
00943 error = -1
00944 return
00945 endif
00946 else
00947 ! create the file
00948 call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
00949 if (error .LT. 0) then
00950 ! close the dataset
00951 call XF_CLOSE_FILE(xFileId, error)
00952 return
00953 endif
00954
00955 ! create the group where we will put all the datasets
00956 call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
00957 if (status < 0) then
00958 call XF_CLOSE_FILE(xFileId, error)
00959 error = -1
00960 return
00961 endif
00962 endif
00963
00964 ! Create/Overwrite the vector dataset group
00965 call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_B_LOCATION, 'ft/s', &
00966 TS_SECONDS, a_Compression, xVector2D_B, status)
00967 if (status .LT. 0) then
00968 ! close the dataset
00969 call XF_CLOSE_GROUP(xVector2D_B, error)
00970 call XF_CLOSE_GROUP(xDsetsId, error)
00971 call XF_CLOSE_FILE(xFileId, error)
00972 error = status
00973 return
00974 endif
00975
00976 if (.NOT. a_Overwrite) then
00977 ! Loop through timesteps adding them to the file
00978 do iTimestep = 1, nTimes
00979 ! We will have an 0.5 hour timestep
00980 dTime = iTimestep * 0.5
00981 do i = 1, nValues
00982 do j = 1, nComponents
00983 fValues(j,i) = ((i-1)*nComponents + j)*dTime
00984 end do
00985 end do
00986 ! write the dataset array values
00987 call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
00988 fValues, error)
00989 if (error .GE. 0) then
00990 ! write activity array
00991 call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
00992 end if
00993 if (error < 0) then
00994 call XF_CLOSE_GROUP(xVector2D_B, error)
00995 call XF_CLOSE_GROUP(xDsetsId, error)
00996 call XF_CLOSE_FILE(xFileId, error)
00997 endif
00998 enddo
00999 else
01000 ! Loop through timesteps adding them to the file
01001 do iTimestep = 1, nTimes
01002 ! We will have an 1.5 hour timestep
01003 dTime = iTimestep * 1.5
01004 do i = 1, nValues
01005 do j = 1, nComponents
01006 fValues(j,i) = ((i-1)*nComponents + j)*dTime
01007 end do
01008 end do
01009 ! write the dataset array values
01010 call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
01011 fValues, error)
01012 if (error .GE. 0) then
01013 ! write activity array
01014 call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
01015 end if
01016 if (error < 0) then
01017 call XF_CLOSE_GROUP(xVector2D_B, error)
01018 call XF_CLOSE_GROUP(xDsetsId, error)
01019 call XF_CLOSE_FILE(xFileId, error)
01020 endif
01021 enddo
01022 endif
01023
01024 if (.NOT. a_Overwrite) then
01025 ! Write Coordinate file - for ScalarB, we will set the coordinate system
01026 ! to be UTM, with UTM Zone settings written to the file.
01027 call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
01028 if (status < 0) then
01029 call XF_CLOSE_GROUP(xVector2D_B, error)
01030 call XF_CLOSE_GROUP(xDsetsId, error)
01031 call XF_CLOSE_FILE(xFileId, error)
01032 error = -1
01033 return
01034 endif
01035
01036 ! write the coordinate data to the file
01037 call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
01038 call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
01039 call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
01040 call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)
01041
01042 ! write additional information - we'll use the max UTM zone for the test
01043 call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)
01044
01045 call XF_CLOSE_GROUP(xCoordId, error)
01046 xCoordId = 0
01047 endif
01048
01049 ! close the dataset
01050 call XF_CLOSE_GROUP(xVector2D_B, error)
01051 call XF_CLOSE_GROUP(xDsetsId, error)
01052 call XF_CLOSE_FILE(xFileId, error)
01053
01054 return
01055 END SUBROUTINE
01056 ! tdWriteVector2D_B
01057
01058 ! --------------------------------------------------------------------------
01059 ! FUNCTION tdWriteVector2D_AToMulti
01060 ! PURPOSE Write scalar Dataset to an HDF5 File
01061 ! NOTES This tests dynamic data sets, and activity
01062 ! This dataset is dynamic concentrations (mg/L) with output times
01063 ! in minutes.
01064 ! Dataset is for a mesh and so nActive is the number of elements
01065 ! which is not the same as the nValues which would be number of nodes
01066 ! reads/writes a reference time in julian days
01067 ! --------------------------------------------------------------------------
01068 SUBROUTINE TD_WRITE_VECTOR2D_A_TO_MULTI (a_FileID, a_GroupID, status)
01069 INTEGER xVector2D_A
01070 INTEGER a_FileID, a_GroupID
01071 INTEGER nValues, nTimes, nComponents, nActive
01072 REAL(DOUBLE) dTime
01073 INTEGER iTimestep, iActive
01074 REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
01075 INTEGER*1 bActivity(100) ! activity
01076 INTEGER i, j, status
01077
01078 ! initialize the data
01079 nComponents = 2
01080 nValues = 100
01081 nTimes = 6
01082 nActive = 75
01083 dTime = 0.0
01084
01085 ! 5th item in data set is always inactive, others active
01086 bActivity(1) = 0
01087 do iActive = 2, nActive
01088 if (mod(iActive-1, 3) == 0) then
01089 bActivity(iActive) = 0
01090 else
01091 bActivity(iActive) = 1
01092 endif
01093 enddo
01094
01095 ! Create the vector dataset group
01096 call XF_CREATE_VECTOR_DATASET(a_GroupID, VECTOR2D_A_LOCATION, 'ft/s', &
01097 TS_SECONDS, NONE, xVector2D_A, status)
01098 if (status .LT. 0) then
01099 ! close the dataset
01100 call XF_CLOSE_GROUP(xVector2D_A, status)
01101 call XF_CLOSE_GROUP(a_GroupID, status)
01102 call XF_CLOSE_FILE(a_FileID, status)
01103 return
01104 endif
01105
01106 ! Loop through timesteps adding them to the file
01107 do iTimestep = 1, nTimes
01108 ! We will have an 0.5 hour timestep
01109 dTime = iTimestep * 0.5
01110
01111 do i = 1, nValues
01112 do j = 1, nComponents
01113 fValues(j,i) = ((i-1)*nComponents + j)*dTime
01114 end do
01115 end do
01116
01117 ! write the dataset array values
01118 call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
01119 fValues, status)
01120 if (status .GE. 0) then
01121 ! write activity array
01122 call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, status)
01123 end if
01124 enddo
01125
01126 ! close the dataset
01127 call XF_CLOSE_GROUP(xVector2D_A, status)
01128 return
01129 END SUBROUTINE
01130 ! tdWriteVector2D_AToMulti
01131 ! --------------------------------------------------------------------------
01132 ! FUNCTION tdiReadScalar
01133 ! PURPOSE Read a scalar from an XMDF file and output information to
01134 ! to a text file
01135 ! NOTES
01136 ! --------------------------------------------------------------------------
01137 SUBROUTINE TDI_READ_SCALAR (a_xScalarId, FileUnit, error)
01138 INTEGER, INTENT(IN) :: a_xScalarId
01139 INTEGER, INTENT(IN) :: FileUnit
01140 INTEGER, INTENT(OUT) :: error
01141 INTEGER nTimes, nValues, nActive
01142 LOGICAL*2 bUseReftime
01143 INTEGER iTime
01144 CHARACTER(LEN=100) TimeUnits
01145 REAL(DOUBLE), ALLOCATABLE :: Times(:)
01146 REAL, ALLOCATABLE :: Values(:), Minimums(:), Maximums(:)
01147 INTEGER, ALLOCATABLE :: Active(:)
01148 REAL(DOUBLE) Reftime
01149 nTimes = NONE
01150 nValues = NONE
01151 nActive = None
01152
01153 ! read the time units
01154 call XF_GET_DATASET_TIME_UNITS(a_xScalarId, TimeUnits, error)
01155 if (error < 0) return
01156
01157 WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))
01158
01159 ! see if we are using a reftime
01160 call XF_USE_REFTIME (a_xScalarId, bUseReftime, error)
01161 if (error < 0) then
01162 return
01163 endif
01164 if (bUseReftime) then
01165 call XF_READ_REFTIME (a_xScalarId, Reftime, error)
01166 if (error < 0) then
01167 return
01168 endif
01169 WRITE(FileUnit,*) 'Reftime: ', Reftime
01170 endif
01171
01172 ! read in the number of values and number of active values
01173 call XF_GET_DATASET_NUMVALS(a_xScalarId, nValues, error)
01174 if (error .GE. 0) then
01175 call XF_GET_DATASET_NUMACTIVE(a_xScalarId, nActive, error)
01176 endif
01177 if (error .LT. 0) return
01178
01179 if (nValues <= 0) then
01180 WRITE(FileUnit, *) 'No data to read in.'
01181 error = -1
01182 return
01183 endif
01184
01185 ! read in the number of times
01186 call XF_GET_DATASET_NUM_TIMES(a_xScalarId, nTimes, error)
01187 if (error < 0) then
01188 return
01189 endif
01190
01191 ! Read in the individual time values
01192 allocate(Times(nTimes))
01193
01194 call XF_GET_DATASET_TIMES(a_xScalarId, nTimes, Times, error)
01195 if (error < 0) return
01196
01197 ! Read in the minimum and maximum values
01198 allocate(Minimums(nTimes))
01199 allocate(Maximums(nTimes))
01200
01201 call XF_GET_DATASET_MINS(a_xScalarId, nTimes, Minimums, error)
01202 if (error >= 0) then
01203 call XF_GET_DATASET_MAXS(a_xScalarId, nTimes, Maximums, error)
01204 endif
01205 if (error < 0) then
01206 deallocate(Times)
01207 deallocate(Minimums)
01208 deallocate(Maximums)
01209 return
01210 endif
01211
01212 allocate(Values(nValues))
01213 if (nActive .GT. 0) then
01214 allocate(Active(nActive))
01215 endif
01216
01217 WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
01218 WRITE(FileUnit,*) 'Number Values: ', nValues
01219 WRITE(FileUnit,*) 'Number Active: ', nActive
01220 WRITE(FileUnit,*) ''
01221
01222 ! loop through the timesteps, read the values and active values and write
01223 ! them to the text file
01224 do iTime = 1, nTimes
01225 call XF_READ_SCALAR_VALUES_TIMESTEP(a_xScalarId, iTime, nValues, Values, error)
01226 if (error >= 0 .AND. nActive > 0) then
01227 call XF_READ_ACTIVITY_TIMESTEP(a_xScalarId, iTime, nActive, Active, error)
01228 endif
01229
01230 ! Write the time, min, max, values and active values to the text output
01231 ! file.
01232 WRITE(FileUnit,*) 'Timestep at ', Times(iTime)
01233 WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
01234 WRITE(FileUnit,*) 'Max: ', Maximums(iTime)
01235
01236 WRITE(FileUnit,*) 'Values:'
01237 WRITE(FileUnit,*) Values(1:nValues)
01238 WRITE(FileUnit,*) ''
01239
01240 WRITE(FileUnit,*) 'Activity:'
01241 WRITE(FileUnit,*) Active(1:nActive)
01242 WRITE(FileUnit,*) ''
01243 end do
01244
01245 if (allocated(Times)) then
01246 deallocate(Times)
01247 endif
01248
01249 if (allocated(Minimums)) then
01250 deallocate(Minimums)
01251 endif
01252
01253 if (allocated(Maximums)) then
01254 deallocate(Maximums)
01255 endif
01256
01257 if (allocated(Values)) then
01258 deallocate(Values)
01259 endif
01260
01261 if (allocated(Active)) then
01262 deallocate(Active)
01263 endif
01264
01265 return
01266 END SUBROUTINE
01267 ! tdiReadScalar
01268
01269 ! --------------------------------------------------------------------------
01270 ! FUNCTION TDI_READ_VECTOR
01271 ! PURPOSE Read a vector from an XMDF file and output information to
01272 ! to a text file
01273 ! NOTES
01274 ! --------------------------------------------------------------------------
01275 SUBROUTINE TDI_READ_VECTOR (a_xVectorId, FileUnit, error)
01276 INTEGER, INTENT(IN) :: a_xVectorId
01277 INTEGER, INTENT(IN) :: FileUnit
01278 INTEGER, INTENT(OUT) :: error
01279 INTEGER nTimes, nValues, nActive, nComponents
01280 INTEGER iTime, i
01281 LOGICAL*2 bUseReftime
01282 CHARACTER(LEN=100) TimeUnits
01283 REAL(DOUBLE), ALLOCATABLE :: Times(:)
01284 REAL, ALLOCATABLE, DIMENSION (:, :) :: Values
01285 REAL, ALLOCATABLE :: Minimums(:), Maximums(:)
01286 INTEGER, ALLOCATABLE :: Active(:)
01287 REAL(DOUBLE) Reftime
01288
01289 nTimes = NONE
01290 nValues = NONE
01291 nActive = NONE
01292 nComponents = NONE
01293
01294 ! read the time units
01295 call XF_GET_DATASET_TIME_UNITS(a_xVectorId, TimeUnits, error)
01296 if (error < 0) return
01297
01298 WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))
01299
01300 ! see if we are using a reftime
01301 call XF_USE_REFTIME (a_xVectorId, bUseReftime, error)
01302 if (error < 0) then
01303 return
01304 endif
01305 if (bUseReftime) then
01306 call XF_READ_REFTIME (a_xVectorId, Reftime, error)
01307 if (error < 0) then
01308 return
01309 endif
01310 WRITE(FileUnit,*) 'Reftime: ', Reftime
01311 endif
01312
01313 ! read in the number of values and number of active values
01314 call XF_GET_DATASET_NUMVALS(a_xVectorId, nValues, error)
01315 if (error .GE. 0) then
01316 call XF_GET_DATASET_NUMCOMPONENTS(a_xVectorId, nComponents, error)
01317 if (error .GE. 0) then
01318 call XF_GET_DATASET_NUMACTIVE(a_xVectorId, nActive, error)
01319 endif
01320 endif
01321 if (error .LT. 0) return
01322
01323 if (nValues <= 0) then
01324 WRITE(FileUnit, *) 'No data to read in.'
01325 error = -1
01326 return
01327 endif
01328
01329 ! read in the number of times
01330 call XF_GET_DATASET_NUM_TIMES(a_xVectorId, nTimes, error)
01331 if (error < 0) then
01332 return
01333 endif
01334
01335 ! Read in the individual time values
01336 allocate(Times(nTimes))
01337
01338 call XF_GET_DATASET_TIMES(a_xVectorId, nTimes, Times, error)
01339 if (error < 0) return
01340
01341 ! Read in the minimum and maximum values
01342 allocate(Minimums(nTimes))
01343 allocate(Maximums(nTimes))
01344
01345 call XF_GET_DATASET_MINS(a_xVectorId, nTimes, Minimums, error)
01346 if (error >= 0) then
01347 call XF_GET_DATASET_MAXS(a_xVectorId, nTimes, Maximums, error)
01348 endif
01349 if (error < 0) then
01350 deallocate(Times)
01351 deallocate(Minimums)
01352 deallocate(Maximums)
01353 return
01354 endif
01355
01356 allocate(Values(nComponents, nValues))
01357 if (nActive .GT. 0) then
01358 allocate(Active(nActive))
01359 endif
01360
01361 WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
01362 WRITE(FileUnit,*) 'Number Values: ', nValues
01363 WRITE(FileUnit,*) 'Number Components: ', nComponents
01364 WRITE(FileUnit,*) 'Number Active: ', nActive
01365
01366 ! loop through the timesteps, read the values and active values and write
01367 ! them to the text file
01368 do iTime = 1, nTimes
01369 call XF_READ_VECTOR_VALUES_TIMESTEP(a_xVectorId, iTime, nValues, &
01370 nComponents, Values, error)
01371 if (error >= 0 .AND. nActive > 0) then
01372 call XF_READ_ACTIVITY_TIMESTEP(a_xVectorId, iTime, nActive, Active, error)
01373 endif
01374
01375 ! Write the time, min, max, values and active values to the text output
01376 ! file.
01377 WRITE(FileUnit,*) ''
01378 WRITE(FileUnit,*) 'Timestep at ', Times(iTime)
01379 WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
01380 WRITE(FileUnit,*) 'Max: ', Maximums(iTime)
01381
01382 WRITE(FileUnit,*) 'Values:'
01383 do i=1, nValues
01384 WRITE(FileUnit,*) Values(1:nComponents,i:i)
01385 enddo
01386 WRITE(FileUnit,*) ''
01387
01388 WRITE(FileUnit,*) 'Activity:'
01389 WRITE(FileUnit,*) Active(1:nActive)
01390 WRITE(FileUnit,*) ''
01391 WRITE(FileUnit,*) ''
01392
01393 end do
01394
01395 if (allocated(Times)) then
01396 deallocate(Times)
01397 endif
01398
01399 if (allocated(Minimums)) then
01400 deallocate(Minimums)
01401 endif
01402
01403 if (allocated(Maximums)) then
01404 deallocate(Maximums)
01405 endif
01406
01407 if (allocated(Values)) then
01408 deallocate(Values)
01409 endif
01410
01411 if (allocated(Active)) then
01412 deallocate(Active)
01413 endif
01414
01415 return
01416 END SUBROUTINE
01417 ! tdiReadVector
01418
01419 END MODULE TestDatasets
00001 MODULE TestGeomPaths
00002
00003 USE TestDatasets
00004 USE Xmdf
00005 USE ErrorDefinitions
00006 USE XmdfDefs
00007
00008 CONTAINS
00009
00010 SUBROUTINE TM_READ_TEST_PATHS(Filename, OutFilename, error)
00011 CHARACTER(LEN=*), INTENT(IN) :: Filename, OutFilename
00012 INTEGER, INTENT(OUT) :: error
00013 INTEGER OutUnit
00014 INTEGER :: xFileId, xGroupId
00015 INTEGER nGroups, nMaxPathLength, nDims, nPaths, nTimes
00016 INTEGER i, j, nStatus
00017 CHARACTER,ALLOCATABLE :: Paths(:)
00018 REAL(DOUBLE), ALLOCATABLE :: times(:), locs(:)
00019 REAL(DOUBLE) Nullval
00020 CHARACTER(LEN=BIG_STRING_SIZE) :: IndividualPath
00021 INTEGER StartLoc
00022 INTEGER, DIMENSION(2) :: PathIndices
00023 INTEGER iTime, iPath
00024
00025 ! open the XMDF file
00026 CALL XF_OPEN_FILE(Filename, .FALSE., xFileId, nStatus)
00027 if (nStatus < 0) then
00028 WRITE(*) 'Unable to open XMDF file TM_READ_TEST_PATHS.'
00029 error = nStatus;
00030 return
00031 endif
00032
00033 ! open the Output file
00034 OutUnit = 79
00035 OPEN(UNIT=OutUnit, FILE=OutFilename, STATUS='REPLACE', ACTION='WRITE', &
00036 IOSTAT = error)
00037 if (OutUnit == 0) then
00038 call XF_CLOSE_FILE(xFileId, error)
00039 error = -1
00040 return
00041 endif
00042
00043 ! find the geomotric path groups
00044 ! Get the number and paths of datasets in the file.
00045 CALL XF_GRP_PTHS_SZ_FOR_GEOM_PTHS(xFileId, nGroups, &
00046 nMaxPathLength, nStatus)
00047 if (nStatus >= 0 .AND. nGroups > 0) then
00048 allocate (Paths(nMaxPathLength*nGroups))
00049 CALL XF_GRP_PTHS_FOR_GEOM_PTHS(xFileId, nGroups, &
00050 nMaxPathLength, Paths, nStatus)
00051 endif
00052 if (nStatus < 0) then
00053 CALL XF_CLOSE_FILE(xFileId, nStatus)
00054 error = -1
00055 return
00056 endif
00057 ! Report the number and paths to individual geom paths groups in the file.
00058 WRITE(OutUnit,*) 'Number of geometric paths in file: ', nGroups
00059
00060 do i = 1, nGroups
00061 StartLoc = (i-1)*nMaxPathLength + 1
00062 IndividualPath = ''
00063 do j = 1, nMaxPathLength - 1
00064 IndividualPath(j:j) = Paths(StartLoc+j-1)
00065 enddo
00066 WRITE(OutUnit,*) 'Reading particles in group: ', &
00067 IndividualPath(1:LEN_TRIM(IndividualPath))
00068
00069 CALL XF_OPEN_GROUP(xFileId, IndividualPath, xGroupId, nStatus)
00070 if (nStatus >= 0) then
00071 ! read the dimensionality of the paths
00072 CALL XF_GET_PATH_DIMENSIONALITY(xGroupId, nDims, nStatus)
00073 if (nStatus >= 0) then
00074 WRITE(OutUnit,*) 'Group dimensionality:', nDims
00075 ! read the number of paths
00076 CALL XF_GET_NUMBER_OF_PATHS(xGroupId, nPaths, nStatus)
00077 if (nStatus >= 0) then
00078 WRITE(OutUnit,*) 'Number of paths in group:', nPaths
00079 ! read the number of timesteps
00080 CALL XF_GET_NUMBER_OF_TIMES(xGroupId, nTimes, nStatus)
00081 if (nStatus >= 0) then
00082 WRITE(OutUnit,*) 'Number of timesteps in group:', nTimes
00083 ! allocate times array
00084 allocate(times(nTimes))
00085 CALL XF_GET_PATH_TIMES_ARRAY(xGroupId, nTimes, times, nStatus)
00086 if (nStatus >= 0) then
00087 ! get space for the data location values
00088 allocate(locs(nPaths*nDims))
00089 ! read null value
00090 CALL XF_GET_PATH_NULL_VAL(xGroupId, NullVal, nStatus)
00091 if (nStatus >= 0) then
00092 WRITE(OutUnit,*) 'Null value:', NullVal
00093 do iTime=1, nTimes
00094 WRITE(OutUnit,*) 'Timestep: ', times(iTime)
00095 ! read the data for this timestep
00096 CALL XF_READ_PATH_LOCATIONS_AT_TIME(xGroupId, iTime, &
00097 1, nPaths, locs, nStatus)
00098 if (nStatus >= 0) then
00099 ! write the data for this timestep
00100 WRITE(OutUnit,*) ' X Y'
00101 if (nDims == 3) then
00102 WRITE(OutUnit,*) ' Z'
00103 endif
00104 WRITE(OutUnit,*) ''
00105 do j=1, nPaths
00106 if (locs(j*nDims) == NullVal) then
00107 WRITE(OutUnit,*) 'Particle not active yet'
00108 else
00109 WRITE(OutUnit,*) locs((j-1)*nDims+1), ' ', locs((j-1)*nDims+2)
00110 if (nDims == 3) then
00111 WRITE(OutUnit,*) ' ', locs((j-1)*nDims+3)
00112 endif
00113 WRITE(OutUnit,*) ''
00114 endif
00115 enddo ! Loop through the paths
00116 endif ! if timestep read
00117 enddo ! Loop through timesteps
00118 endif ! if null value read
00119 if (allocated(locs)) deallocate (locs)
00120 endif ! if we get get the times array
00121 ! get space for the data location values - 1 particle
00122 ! all times
00123 allocate(locs(nTimes*nDims))
00124 do iPath=1, nPaths
00125 ! read the values for this particle for all timesteps
00126 CALL XF_READ_PATH_LOCS_FOR_PART(xGroupId, iPath, &
00127 1, nTimes, locs, nStatus)
00128 if (nStatus >= 0) then
00129 ! write the data for this path
00130 WRITE(OutUnit,*) 'Time X Y'
00131 if (nDims == 3) then
00132 WRITE(OutUnit,*) ' Z'
00133 endif
00134 WRITE(OutUnit, *) ''
00135 do j=1, nTimes
00136 if (locs((j-1)*nDims+1) .NE. NullVal) then
00137 WRITE(OutUnit,*) times(j), ' ', locs((j-1)*nDims + 1), ' ', &
00138 locs((j-1)*nDims+2)
00139 if (nDims == 3) then
00140 WRITE(OutUnit,*) ' ', locs((j-1)*nDims+3)
00141 endif
00142 WRITE(OutUnit,*) ''
00143 endif
00144 enddo ! Loop through the times
00145 endif ! if read path locations for particle
00146 enddo ! Loop through the paths
00147 if (allocated(locs)) deallocate(locs)
00148 endif ! if allocation of locs suceeded
00149 ! get space for the data location values - 2 particle
00150 ! all times
00151 allocate(locs(nTimes*nDims*2))
00152 PathIndices(1) = 1
00153 PathIndices(2) = nPaths
00154 CALL XF_READ_PATH_LOCS_FOR_PARTS(xGroupId, 2, &
00155 PathIndices, 1, nTimes, locs, nStatus)
00156 if (nStatus >= 0) then
00157 ! write the data for these 2 paths
00158 if (nDims == 3) then
00159 WRITE(OutUnit,*) 'Timestep X1 Y1 Z1 Xn Yn Zn'
00160 else
00161 WRITE(OutUnit,*) 'Timestep X1 Y1 Xn Yn'
00162 endif
00163 do j=1, nTimes
00164 if (nDims == 3) then
00165 WRITE(OutUnit,*) times(j), ' ', locs((j-1)*2*nDims+1), ' ', &
00166 locs((j-1)*2*nDims+2),' ', locs((j-1)*2*nDims+3),' ', &
00167 locs((j-1)*2*nDims+4), locs((j-1)*2*nDims+5), ' ',locs((j-1)*2*nDims+6)
00168 else
00169 WRITE(OutUnit,*) times(j),' ', locs((j-1)*2*nDims + 1), ' ', &
00170 locs((j-1)*2*nDims+2), ' ',locs((j-1)*2*nDims+3)
00171 endif
00172 enddo
00173 else
00174 WRITE(*,*) 'Error reading path locations for multiple particles'
00175 endif
00176 endif
00177 if (allocated(locs)) deallocate(locs)
00178 if (allocated(times)) deallocate(times)
00179 CALL XF_CLOSE_GROUP(xGroupId, nStatus)
00180 if (nStatus < 0) then
00181 WRITE(*)'Error reading geometric paths..'
00182 endif
00183 endif
00184 endif
00185 enddo ! loop through groups
00186 ! free the paths
00187 if (allocated(Paths)) deallocate(Paths)
00188 ! close the files
00189 close(OutUnit)
00190 CALL XF_CLOSE_FILE(xFileId, nStatus)
00191 error = nStatus;
00192 return
00193
00194 END SUBROUTINE
00195 ! tmReadTestPaths
00196
00197 SUBROUTINE TM_WRITE_TEST_PATHS(Filename, Compression, error)
00198 CHARACTER(LEN=*), INTENT(IN) :: Filename
00199 INTEGER, INTENT(IN) :: Compression
00200 INTEGER, INTENT(OUT) :: error
00201 INTEGER nPaths
00202 REAL(DOUBLE), DIMENSION(6) :: pLocs
00203 REAL, DIMENSION(2) :: pSpeeds
00204 REAL(DOUBLE) NullVal
00205
00206 REAL NullValSinglePrec
00207 INTEGER xFileId, xPathGroupId, xSpeedId, xPropGroupId
00208 INTEGER status
00209 REAL(DOUBLE) Timeval
00210
00211 nPaths = 0
00212 NullVal = -99999.9d0
00213 NullValSinglePrec = -99999.9
00214
00215 ! create the file
00216 call XF_CREATE_FILE(Filename, .TRUE., xFileId, status)
00217 if (status < 0) then
00218 error = -1
00219 return
00220 endif
00221
00222 ! create the group to store the particle paths
00223 call XF_CREATE_GEOMETRIC_PATH_GROUP(xFileId, "particles", &
00224 'abcdefglasdfjoaieur', Compression, xPathGroupId, NullVal, status)
00225 if (status < 0) then
00226 error = -1
00227 return
00228 endif
00229
00230 ! create the data set to store the speed
00231 call XF_CREATE_SCALAR_DSET_EXTNDBL(xPathGroupId, 'Vmag', 'm/s', &
00232 TS_SECONDS, NullValSinglePrec, Compression, xSpeedId, status)
00233 if (status < 0) then
00234 error = -1
00235 return
00236 endif
00237
00238 call XF_CREATE_PROPERTY_GROUP(xSpeedId, xPropGroupId, status)
00239 if (status < 0) then
00240 call XF_CLOSE_GROUP(xSpeedId, status)
00241 error = -1
00242 return
00243 endif
00244 call XF_WRITE_PROPERTY_FLOAT(xPropGroupId, PROP_NULL_VALUE, 1, &
00245 NullValSinglePrec, -1, status)
00246 call XF_CLOSE_GROUP(xPropGroupId, status)
00247
00248 ! Setup the arrays for the path group at timestep 0
00249 ! particle location at the first timestep
00250 nPaths = 1
00251 pLocs(1) = 1.0
00252 pLocs(2) = 2.0
00253 pLocs(3) = 3.0
00254
00255 ! store the particles for the first timestep
00256 Timeval = 0.0
00257 call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
00258 if (status < 0) then
00259 error = -1
00260 return
00261 endif
00262 ! set up and store the speed at timestep 0
00263 pSpeeds(1) = 1.1
00264 Timeval = 0.0
00265 call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 1, pSpeeds, status)
00266 if (status < 0) then
00267 error = -1
00268 return
00269 endif
00270
00271 ! Setup the arrays for the path group at timestep 1
00272 ! particle location at the first timestep
00273 pLocs(1) = 4.0
00274 pLocs(2) = 5.0
00275 pLocs(3) = 6.0
00276
00277 ! store the particles for the second timestep
00278 Timeval = 1.0
00279 call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
00280 if (status < 0) then
00281 error = -1
00282 return
00283 endif
00284
00285 ! set up and store the speed at timestep 2
00286 pSpeeds(1) = 1.2
00287 call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 1, pSpeeds, status)
00288 if (status < 0) then
00289 error = -1
00290 return
00291 endif
00292
00293 ! Setup the arrays for the path group at timestep 3-add a particle
00294 ! particle location at the first timestep
00295 nPaths = 2
00296 pLocs(1) = 7.0
00297 pLocs(2) = 8.0
00298 pLocs(3) = 9.0
00299 pLocs(4) = -1.0
00300 pLocs(5) = -2.0
00301 pLocs(6) = -3.0
00302
00303 ! store the particles for the timestep 3
00304 Timeval = 2.0
00305 call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
00306 if (status < 0) then
00307 error = -1
00308 return
00309 endif
00310
00311 ! extend the data set for speed
00312 call XF_EXTEND_SCALAR_DATASET(xSpeedId, 2, status)
00313 if (status < 0) then
00314 error = -1
00315 return
00316 endif
00317
00318 ! set up and store the speed at timestep 4
00319 pSpeeds(1) = 1.3
00320 pSpeeds(2) = 2.1
00321 call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 2, pSpeeds, status)
00322 if (status < 0) then
00323 error = -1
00324 return
00325 endif
00326
00327 ! Setup the arrays for the path group at timestep 3-inactive particle(static)
00328 ! particle location at the first timestep
00329 pLocs(1) = 7.0
00330 pLocs(2) = 8.0
00331 pLocs(3) = 9.0
00332 pLocs(4) = -4.0
00333 pLocs(5) = -5.0
00334 pLocs(6) = -6.0
00335
00336 ! store the particles for timestep 4
00337 Timeval = 3.0
00338 call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
00339 if (status < 0) then
00340 error = -1
00341 return
00342 endif
00343
00344 ! set up and store the speed at timestep 4
00345 pSpeeds(1) = NullVal
00346 pSpeeds(2) = 2.2
00347 call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 2, pSpeeds, status)
00348 if (status < 0) then
00349 error = -1
00350 return
00351 endif
00352
00353 ! close the resources
00354 call XF_CLOSE_GROUP(xSpeedId, status)
00355 call XF_CLOSE_GROUP(xPathGroupId, status)
00356 call XF_CLOSE_FILE(xFileId, status)
00357
00358 error = 1
00359 return
00360 END SUBROUTINE
00361
00362 END MODULE TestGeomPaths
00001 MODULE TestGrid
00002
00003 USE ErrorDefinitions
00004 USE XmdfDefs
00005 USE Xmdf
00006
00007 CHARACTER(LEN=*),PARAMETER :: GRID_CART2D_GROUP_NAME = 'Grid Cart2D Group'
00008 CHARACTER(LEN=*),PARAMETER :: GRID_CURV2D_GROUP_NAME = 'Grid Curv2D Group'
00009 CHARACTER(LEN=*),PARAMETER :: GRID_CART3D_GROUP_NAME = 'Grid Cart3D Group'
00010
00011 CONTAINS
00012
00013 !****************************
00014 ! ---------------------------------------------------------------------------
00015 ! FUNCTION TG_READ_GRID
00016 ! PURPOSE Read a grid and write data to a text file
00017 ! NOTES
00018 ! ---------------------------------------------------------------------------
00019 SUBROUTINE TG_READ_GRID(a_Id, a_Outfile, error)
00020 INTEGER, INTENT(IN) :: a_Id
00021 INTEGER, INTENT(IN) :: a_Outfile
00022 INTEGER, INTENT(OUT) :: error
00023 INTEGER nGridType, nExtrudeType, nDims, nCellsI, nCellsJ
00024 INTEGER nCellsK, nLayers, nOrientation, nValsI, nValsJ
00025 INTEGER nValsK, nExtrudeVals, nCompOrigin, nUDir
00026 CHARACTER(265) strGridType, strExtrudeType
00027 LOGICAL*2 bDefined
00028 REAL(DOUBLE) dOrigin(3), dBearing, dDip, dRoll
00029 REAL(DOUBLE), ALLOCATABLE :: dExtrudeVals(:), dCoordI(:), dCoordJ(:)
00030 REAL(DOUBLE), ALLOCATABLE :: dCoordK(:)
00031 INTEGER i
00032
00033 nGridType = 0
00034 nExtrudeType = 0
00035 nDims = 0
00036 nCellsI = 0
00037 nCellsJ = 0
00038 nCellsK = 0
00039 nLayers = 0
00040 nOrientation = 0
00041 nValsI = 0
00042 nValsJ = 0
00043 nValsK = 0
00044 nExtrudeVals = 0
00045 nCompOrigin = 1
00046 nUDir = 1
00047 bDefined = .FALSE.
00048 dBearing = 0.0
00049 dDip = 0.0
00050 dRoll =0.0
00051 i = 0
00052 error = 1
00053
00054 ! Grid type
00055 call XF_GET_GRID_TYPE(a_Id, nGridType, error)
00056 if (error < 0) then
00057 return
00058 endif
00059 SELECT CASE (nGridType)
00060 CASE (GRID_TYPE_CARTESIAN)
00061 strGridType = 'Cartesian'
00062 CASE (GRID_TYPE_CURVILINEAR)
00063 strGridType = 'Curvilinear'
00064 CASE (GRID_TYPE_CARTESIAN_EXTRUDED)
00065 strGridType = 'Cartesian extruded'
00066 CASE (GRID_TYPE_CURVILINEAR_EXTRUDED)
00067 strGridType = 'Curvilinear extruded'
00068 CASE DEFAULT
00069 WRITE(*,*) 'Invalid grid type'
00070 error = -1
00071 END SELECT
00072 WRITE(a_Outfile,*) 'The grid type is: ', strGridType(1:LEN_TRIM(strGridType))
00073
00074 ! Number of dimensions
00075 call XF_GET_NUMBER_OF_DIMENSIONS(a_Id, nDims, error)
00076 if (error .LT. 0) then
00077 return
00078 endif
00079 if (nDims == 2) then
00080 WRITE(a_Outfile,*) 'The grid is two-dimensional'
00081 elseif (nDims == 3) then
00082 WRITE(a_Outfile,*) 'The grid is three-dimensional'
00083 else
00084 WRITE(*,*) 'The grid dimensions are invalid'
00085 error = -1
00086 endif
00087
00088 ! Extrusion type if applicable
00089 if (nGridType .EQ. GRID_TYPE_CARTESIAN_EXTRUDED .OR. &
00090 nGridType .EQ. GRID_TYPE_CURVILINEAR_EXTRUDED) then
00091 call XF_GET_EXTRUSION_TYPE(a_Id, nExtrudeType, error)
00092 if (error < 0) then
00093 return
00094 endif
00095 SELECT CASE (nExtrudeType)
00096 case (EXTRUDE_SIGMA)
00097 strExtrudeType = 'Sigma stretch'
00098 case (EXTRUDE_CARTESIAN)
00099 strExtrudeType = 'Cartesian'
00100 case (EXTRUDE_CURV_AT_CORNERS)
00101 strExtrudeType = 'Curvilinear at Corners'
00102 case (EXTRUDE_CURV_AT_CELLS)
00103 strExtrudeType = 'Curvilinear at Cells'
00104 END SELECT
00105 WRITE(a_Outfile,*) 'The grid is extruding using: ', &
00106 strExtrudeType(1:LEN_TRIM(strExtrudeType))
00107 endif
00108
00109 ! Origin
00110 call XF_ORIGIN_DEFINED(a_Id, bDefined, error)
00111 if (error < 0) then
00112 return
00113 endif
00114 if (bDefined) then
00115 call XF_GET_ORIGIN(a_Id, dOrigin(1), dOrigin(2), dOrigin(3), error)
00116 if (error < 0) then
00117 return
00118 endif
00119 WRITE(a_Outfile,*) 'The grid origin is ', dOrigin(1), ' ',&
00120 dOrigin(2), ' ', dOrigin(3)
00121 endif
00122
00123 ! Orientation
00124 call XF_GET_ORIENTATION(a_Id, nOrientation, error)
00125 if (error < 0) then
00126 return
00127 endif
00128 if (nOrientation == ORIENTATION_RIGHT_HAND) then
00129 WRITE(a_Outfile,*) 'The grid has a right hand orientation'
00130 elseif (nOrientation == ORIENTATION_LEFT_HAND) then
00131 WRITE(a_Outfile,*) 'The grid has a left hand orientation'
00132 else
00133 WRITE(*,*) 'Invalid grid orientation';
00134 error = -1
00135 return
00136 endif
00137
00138 ! Bearing
00139 call XF_BEARING_DEFINED(a_Id, bDefined, error)
00140 if (error < 0) then
00141 return
00142 endif
00143 if (bDefined) then
00144 call XF_GET_BEARING(a_Id, dBearing, error)
00145 if (error < 0) then
00146 return
00147 endif
00148 WRITE(a_Outfile,*) 'The grid bearing is ', dBearing
00149 endif
00150
00151 ! Dip
00152 call XF_DIP_DEFINED(a_Id, bDefined, error)
00153 if (error < 0) then
00154 return
00155 endif
00156 if (bDefined) then
00157 call XF_GET_DIP(a_Id, dDip, error)
00158 if (error < 0) then
00159 return
00160 endif
00161 WRITE(a_Outfile,*) 'The grid Dip is ', dDip
00162 endif
00163
00164 if (nDims == 3) then
00165 ! Roll
00166 call XF_ROLL_DEFINED(a_Id, bDefined, error)
00167 if (error < 0) then
00168 return
00169 endif
00170 if (bDefined) then
00171 call XF_GET_ROLL(a_Id, dRoll, error)
00172 if (error < 0) then
00173 return
00174 endif
00175 WRITE(a_Outfile,*) 'The grid Roll is ', dRoll
00176 endif
00177 endif
00178
00179 ! Computational origin
00180 call XF_COMPUTATIONAL_ORIGIN_DEFINED(a_Id, bDefined, error)
00181 if (error < 0) then
00182 return
00183 endif
00184 if (bDefined) then
00185 call XF_GET_COMPUTATIONAL_ORIGIN(a_Id, nCompOrigin, error)
00186 if (error < 0) then
00187 return
00188 endif
00189 WRITE(a_Outfile,*) 'The grid Computational Origin is ', nCompOrigin
00190 else
00191 WRITE(a_Outfile,*) 'The grid Computational Origin is not defined';
00192 endif
00193
00194
00195 ! U Direction
00196 call XF_GET_U_DIRECTION_DEFINED(a_Id, bDefined, error)
00197 if (error < 0) then
00198 return
00199 endif
00200 if (bDefined) then
00201 call XF_GET_U_DIRECTION(a_Id, nUDir, error)
00202 if (error < 0) then
00203 return
00204 endif
00205 WRITE(a_Outfile,*) 'The grid U Direction is ', nUDir
00206 else
00207 WRITE(a_Outfile,*) 'The grid U Direction is not defined'
00208 endif
00209
00210 ! number of cells in each direction
00211 call XF_GET_NUMBER_CELLS_IN_I(a_Id, nCellsI, error)
00212 if (error >= 0) then
00213 call XF_GET_NUMBER_CELLS_IN_J(a_Id, nCellsJ, error)
00214 if ((error >= 0) .AND. (nDims == 3)) then
00215 call XF_GET_NUMBER_CELLS_IN_K(a_Id, nCellsK, error)
00216 endif
00217 endif
00218 if (error < 0) then
00219 return
00220 endif
00221 WRITE(a_Outfile,*) 'Number of cells in I ', nCellsI
00222 WRITE(a_Outfile,*) 'Number of cells in J ', nCellsJ
00223 if (nDims == 3) then
00224 WRITE(a_Outfile,*) 'Number of cells in K ', nCellsK
00225 endif
00226
00227 ! Grid coordinates
00228 if (nGridType == GRID_TYPE_CARTESIAN .OR. &
00229 nGridType == GRID_TYPE_CARTESIAN_EXTRUDED) then
00230 nValsI = nCellsI
00231 nValsJ = nCellsJ
00232 if (nDims == 3) then
00233 nValsK = nCellsK
00234 endif
00235 elseif (nGridType == GRID_TYPE_CURVILINEAR .OR. &
00236 nGridType == GRID_TYPE_CURVILINEAR_EXTRUDED) then
00237 if (nDims == 3) then
00238 ! three dimensions
00239 nValsK = (nCellsI + 1) * (nCellsJ + 1) * (nCellsK + 1)
00240 nValsJ = nValsK
00241 nValsI = nValsJ
00242 else
00243 ! two dimensions
00244 nValsJ = (nCellsI + 1) * (nCellsJ + 1)
00245 nValsI = nValsJ
00246 endif
00247 else
00248 WRITE(*,*) 'Invalid grid type'
00249 error = -1
00250 return
00251 endif
00252
00253 ALLOCATE(dCoordI(nValsI))
00254 ALLOCATE(dCoordJ(nValsJ))
00255 if (nDims == 3) then
00256 ALLOCATE(dCoordK(nValsK))
00257 endif
00258
00259 call XF_GET_GRID_COORDS_I(a_Id, nValsI, dCoordI, error)
00260 if (error >= 0) then
00261 call XF_GET_GRID_COORDS_J(a_Id, nValsJ, dCoordJ, error)
00262 if ((error >= 0) .AND. (nDims == 3)) then
00263 call XF_GET_GRID_COORDS_K(a_Id, nValsK, dCoordK, error)
00264 endif
00265 endif
00266 if (error < 0) then
00267 WRITE(*,*) 'Error reading coordinates'
00268 error = -1
00269 return
00270 endif
00271
00272 WRITE(a_Outfile,*) 'The Coordinates in direction I:'
00273 do i = 1, nValsI
00274 if (mod(i,5) == 0) then
00275 WRITE(a_Outfile,*) ''
00276 endif
00277 WRITE(a_Outfile,*) dCoordI(i)
00278 enddo
00279 WRITE(a_Outfile,*) ''
00280
00281 WRITE(a_Outfile,*) 'The Coordinates in direction J:'
00282 do i = 1, nValsJ
00283 if (mod(i,5) == 0) then
00284 WRITE(a_Outfile,*) ''
00285 endif
00286 WRITE(a_Outfile,*) dCoordJ(i)
00287 enddo
00288 WRITE(a_Outfile,*) ''
00289
00290 if (nDims == 3) then
00291 WRITE(a_Outfile,*) 'The Coordinates in direction K:'
00292 do i = 1, nValsK
00293 if (mod(i,5) == 0) then
00294 WRITE(a_Outfile,*) ''
00295 endif
00296 WRITE(a_Outfile,*) dCoordK(i)
00297 enddo
00298 endif
00299 WRITE(a_Outfile,*) ''
00300
00301 if (ALLOCATED(dCoordI)) DEALLOCATE(dCoordI)
00302 if (ALLOCATED(dCoordJ)) DEALLOCATE(dCoordJ)
00303 if (ALLOCATED(dCoordK)) DEALLOCATE(dCoordK)
00304
00305 !
00306 if (nGridType .EQ. GRID_TYPE_CARTESIAN_EXTRUDED .OR. &
00307 nGridType .EQ. GRID_TYPE_CURVILINEAR_EXTRUDED) then
00308 call XF_GET_EXTRUDE_NUM_LAYERS(a_Id, nLayers, error)
00309 if (error < 0) then
00310 return
00311 endif
00312
00313 SELECT CASE(nExtrudeType)
00314 case (EXTRUDE_SIGMA)
00315 nExtrudeVals = nLayers
00316 case (EXTRUDE_CURV_AT_CORNERS)
00317 nExtrudeVals = (nCellsI + 1) * (nCellsJ + 1) * nLayers
00318 case (EXTRUDE_CURV_AT_CELLS)
00319 nExtrudeVals = nCellsI * nCellsJ * nLayers
00320 END SELECT
00321
00322 ALLOCATE(dExtrudeVals(nExtrudeVals))
00323
00324 call XF_GET_EXTRUDE_VALUES(a_Id, nExtrudeVals, dExtrudeVals, error)
00325 if (error < 0) then
00326 return
00327 endif
00328
00329 WRITE(*,*) 'The extrude values are:'
00330 do i = 1, nExtrudeVals
00331 if (mod(i,5) == 0) then
00332 WRITE(a_Outfile,*) ''
00333 endif
00334 WRITE(a_Outfile,*) dExtrudeVals(i)
00335 enddo
00336 if (ALLOCATED(dExtrudeVals)) DEALLOCATE(dExtrudeVals)
00337 endif
00338
00339 return
00340
00341 END SUBROUTINE TG_READ_GRID
00342
00343 !****************************
00344 !----------------------------------------------------------------------------
00345 ! FUNCTION TG_WRITE_TEST_GRID_CART_2D
00346 ! PURPOSE Write a file that contains data for a 2D Cartesian Grid
00347 ! NOTES A picture of the grid is in the file (TestGridCart2D.gif)
00348 ! returns TRUE on success and FALSE on failure
00349 !----------------------------------------------------------------------------
00350 SUBROUTINE TG_WRITE_TEST_GRID_CART_2D(Filename, error)
00351 CHARACTER(LEN=*), INTENT(IN) :: Filename
00352 INTEGER, INTENT(OUT) :: error
00353 INTEGER nDimensions
00354 INTEGER nCellsI, nCellsJ
00355 INTEGER nGridType
00356 INTEGER nCompOrigin, nUDir
00357 REAL(DOUBLE) dOriginX, dOriginY, dOriginZ
00358 INTEGER nOrientation
00359 REAL(DOUBLE) dBearing
00360 REAL(DOUBLE) PlanesI(5), PlanesJ(5)
00361 INTEGER i, j, iSpcZone
00362 INTEGER xFileId, xGridId, xCoordId
00363 INTEGER status
00364 INTEGER tmpOut1, tmpOut2
00365
00366 nDimensions = 2
00367 nCellsI = 5
00368 nCellsJ = 5
00369 nGridType = GRID_TYPE_CARTESIAN
00370 nCompOrigin = 4
00371 nUDir = -2
00372 dOriginX = 10.0
00373 dOriginY = 10.0
00374 dOriginZ = 0.0
00375 nOrientation = ORIENTATION_RIGHT_HAND
00376 dBearing = 45.0
00377 xFileId = NONE
00378 xGridId = NONE
00379
00380 ! Fill in the grid plane data with a constant size of 30
00381 do i = 1, nCellsI
00382 PlanesI(i) = i*30.0
00383 enddo
00384 do j = 1, nCellsJ
00385 PlanesJ(j) = j*30.0
00386 enddo
00387
00388 ! create the file
00389 call XF_CREATE_FILE(Filename, .TRUE., xFileId, status)
00390 if (status < 0) then
00391 error = -1
00392 return
00393 endif
00394
00395 ! create the group to store the grid
00396 call XF_CREATE_GROUP_FOR_GRID(xFileId, GRID_CART2D_GROUP_NAME, xGridId, status)
00397 if (status < 0) then
00398 call XF_CLOSE_FILE(xFileId, error)
00399 error = -1
00400 return
00401 endif
00402
00403 ! Write the grid information to the file
00404 call XF_SET_GRID_TYPE(xGridId, nGridType, tmpOut1)
00405 call XF_SET_NUMBER_OF_DIMENSIONS(xGridId, nDimensions, tmpOut2)
00406
00407 if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
00408 call XF_CLOSE_GROUP(xGridId, error)
00409 call XF_CLOSE_FILE(xFileId, error)
00410 error = -1
00411 return
00412 endif
00413
00414 ! set origin and orientation
00415 call XF_SET_ORIGIN(xGridId, dOriginX, dOriginY, dOriginZ, tmpOut1)
00416 call XF_SET_ORIENTATION(xGridId, nOrientation, tmpOut2)
00417
00418 if ((tmpOut1 < 0) .OR. (tmpOut2 .LT. 0)) then
00419 call XF_CLOSE_GROUP(xGridId, error)
00420 call XF_CLOSE_FILE(xFileId, error)
00421 error = -1
00422 return
00423 endif
00424
00425 ! Set bearing
00426 call XF_SET_BEARING(xGridId, dBearing, tmpOut1)
00427 if (tmpOut1 < 0) then
00428 call XF_CLOSE_GROUP(xGridId, error)
00429 call XF_CLOSE_FILE(xFileId, error)
00430 error = -1
00431 return
<