1 !/***************************************************************************
2 !* The HDF5 WRF IO module was written by the the HDF Group at NCSA, the *
3 !* National Center for Supercomputing Applications. *
5 !* National Center for Supercomputing Applications *
6 !* University of Illinois at Urbana-Champaign *
7 !* 605 E. Springfield, Champaign IL 61820 *
8 !* http://hdf.ncsa.uiuc.edu/ *
10 !* Copyright 2004 by the Board of Trustees, University of Illinois, *
12 !* Redistribution or use of this IO module, with or without modification, *
13 !* is permitted for any purpose, including commercial purposes. *
15 !* This software is an unsupported prototype. Use at your own risk. *
16 !* http://hdf.ncsa.uiuc.edu/apps/WRF-ROMS *
18 !* This work was funded by the MEAD expedition at the National Center *
19 !* for Supercomputing Applications, NCSA. For more information see: *
20 !* http://www.ncsa.uiuc.edu/expeditions/MEAD *
23 !****************************************************************************/
28 integer , parameter :: FATAL = 1
29 integer , parameter :: WARN = 1
30 integer , parameter :: WrfDataHandleMax = 99
31 integer , parameter :: MaxDims = 2000 ! = NF_MAX_VARS
32 integer , parameter :: MaxTabDims = 100 ! temporary,changable
33 integer , parameter :: MaxVars = 2000
34 integer , parameter :: MaxTimes = 9999 ! temporary, changable
35 integer , parameter :: MaxTimeSLen = 6 ! not exceed 1,000,000 timestamp
36 integer , parameter :: DateStrLen = 19
37 integer , parameter :: VarNameLen = 31
38 integer , parameter :: NO_DIM = 0
39 integer , parameter :: NVarDims = 4
40 integer , parameter :: NMDVarDims = 2
41 integer , parameter :: CompDsetSize = 64256 ! set to 63K
42 character (8) , parameter :: NO_NAME = 'NULL'
43 character(4) , parameter :: hdf5_true ='TRUE'
44 character(5) , parameter :: hdf5_false ='FALSE'
45 integer , parameter :: MemOrdLen = 3
46 character (DateStrLen) , parameter :: ZeroDate = '0000-00-00-00:00:00'
48 #include "wrf_io_flags.h"
49 ! This is a hack. WRF IOAPI no longer supports WRF_CHARACTER. Rip this out!
50 integer, parameter :: WRF_CHARACTER = 1080
52 character (120) :: msg
54 ! derived data type for dimensional table
56 character (len = 256) :: dim_name
61 type :: wrf_phdf5_data_handle
62 character (256) :: FileName
65 integer(hid_t) :: FileID
66 integer(hid_t) :: GroupID
67 integer(hid_t) :: DimGroupID
68 integer(hid_t) :: EnumID
69 character (256) :: GroupName
70 character (256) :: DimGroupName
73 character (5) :: TimesName
75 integer :: MaxTimeCount
76 integer :: CurrentTime !Only used for read
77 integer :: NumberTimes !Only used for read
78 character (DateStrLen), pointer :: Times(:)
79 integer(hid_t) :: TimesID
80 integer(hid_t) :: str_id
81 integer , pointer :: DimLengths(:)
82 integer , pointer :: DimIDs(:)
83 character (31) , pointer :: DimNames(:)
85 character (9) :: DimUnlimName
86 type (dim_scale) , pointer :: DIMTABLE(:)
87 integer , dimension(NVarDims) :: DimID
88 integer , dimension(NVarDims) :: Dimension
89 ! integer , pointer :: MDDsetIDs(:)
90 integer , pointer :: MDVarDimLens(:)
91 character (256) , pointer :: MDVarNames(:)
92 integer(hid_t) , pointer :: TgroupIDs(:)
93 integer(hid_t) , pointer :: DsetIDs(:)
94 integer(hid_t) , pointer :: MDDsetIDs(:)
95 ! integer(hid_t) :: DimTableID
96 integer , pointer :: VarDimLens(:,:)
97 character (VarNameLen), pointer :: VarNames(:)
98 integer :: CurrentVariable !Only used for read
100 ! first_operation is set to .TRUE. when a new handle is allocated
101 ! or when open-for-write or open-for-read are committed. It is set
102 ! to .FALSE. when the first field is read or written.
103 logical :: first_operation
104 end type wrf_phdf5_data_handle
105 type(wrf_phdf5_data_handle),target :: WrfDataHandles(WrfDataHandleMax)
107 end module wrf_phdf5_data
110 module ext_phdf5_support_routines
116 subroutine allocHandle(DataHandle,DH,Comm,Status)
120 include 'wrf_status_codes.h'
122 integer ,intent(out) :: DataHandle
123 type(wrf_phdf5_data_handle),pointer:: DH
124 integer ,intent(IN) :: Comm
125 integer ,intent(out) :: Status
129 integer(hid_t) :: enum_type
130 ! character (256) :: NullName
134 do i=1,WrfDataHandleMax
135 if(WrfDataHandles(i)%Free) then
136 DH => WrfDataHandles(i)
144 call SetUp_EnumID(enum_type,Status)
146 Status = WRF_HDF5_ERR_ALLOCATION
147 write(msg,*) 'Fatal enum ALLOCATION ERROR in ',__FILE__,', line',__LINE__
148 call wrf_debug ( FATAL , msg)
151 DH%EnumID = enum_type
153 allocate(DH%Times(MaxTimes), STAT=stat)
155 Status = WRF_HDF5_ERR_ALLOCATION
156 write(msg,*) 'Fatal ALLOCATION ERROR in ',__FILE__,', line',__LINE__
157 call wrf_debug ( FATAL , msg)
161 ! DH%Times(1:MaxTimes) = NullName
163 allocate(DH%DimLengths(MaxDims), STAT=stat)
165 Status = WRF_HDF5_ERR_ALLOCATION
166 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line',__LINE__
168 call wrf_debug ( FATAL , msg)
172 allocate(DH%DimIDs(MaxDims), STAT=stat)
174 Status = WRF_HDF5_ERR_ALLOCATION
175 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
176 call wrf_debug ( FATAL , msg)
180 allocate(DH%DimNames(MaxDims), STAT=stat)
182 Status = WRF_HDF5_ERR_ALLOCATION
183 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
184 call wrf_debug ( FATAL , msg)
188 allocate(DH%DIMTABLE(MaxTabDims), STAT = stat)
190 Status = WRF_HDF5_ERR_ALLOCATION
191 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
192 call wrf_debug ( FATAL , msg)
197 DH%DIMTABLE(j)%dim_name = NO_NAME
198 DH%DIMTABLE(j)%unlimited = -1
201 allocate(DH%MDDsetIDs(MaxVars), STAT=stat)
203 Status = WRF_HDF5_ERR_ALLOCATION
204 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
205 call wrf_debug ( FATAL , msg)
209 allocate(DH%MDVarDimLens(MaxVars), STAT=stat)
211 Status = WRF_HDF5_ERR_ALLOCATION
212 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
213 call wrf_debug ( FATAL , msg)
217 allocate(DH%MDVarNames(MaxVars), STAT=stat)
219 Status = WRF_HDF5_ERR_ALLOCATION
220 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
221 call wrf_debug ( FATAL , msg)
225 allocate(DH%DsetIDs(MaxVars), STAT=stat)
227 Status = WRF_HDF5_ERR_ALLOCATION
228 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
229 call wrf_debug ( FATAL , msg)
234 allocate(DH%TgroupIDs(MaxTimes), STAT=stat)
236 Status = WRF_HDF5_ERR_ALLOCATION
237 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
238 call wrf_debug ( FATAL , msg)
243 allocate(DH%VarDimLens(NVarDims-1,MaxVars), STAT=stat)
245 Status = WRF_HDF5_ERR_ALLOCATION
246 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
247 call wrf_debug ( FATAL , msg)
251 allocate(DH%VarNames(MaxVars), STAT=stat)
253 Status = WRF_HDF5_ERR_ALLOCATION
254 write(msg,*) 'Fatal ALLOCATION ERROR in ',"__FILE__",', line', __LINE__
255 call wrf_debug ( FATAL , msg)
261 if(i==WrfDataHandleMax) then
262 Status = WRF_HDF5_ERR_TOO_MANY_FILES
263 write(msg,*) 'Warning TOO MANY FILES in ',"__FILE__",', line', __LINE__
264 call wrf_debug ( WARN , msg)
273 DH%first_operation = .TRUE.
275 end subroutine allocHandle
277 ! Obtain data handler
278 subroutine GetDH(DataHandle,DH,Status)
281 include 'wrf_status_codes.h'
282 integer ,intent(in) :: DataHandle
283 type(wrf_phdf5_data_handle) ,pointer :: DH
284 integer ,intent(out) :: Status
286 if(DataHandle < 1 .or. DataHandle > WrfDataHandleMax) then
287 Status = WRF_HDF5_ERR_BAD_DATA_HANDLE
290 DH => WrfDataHandles(DataHandle)
292 Status = WRF_HDF5_ERR_BAD_DATA_HANDLE
299 ! Set up eumerate datatype for possible logical type
300 subroutine SetUp_EnumID(enum_type,Status)
305 include 'wrf_status_codes.h'
306 integer(hid_t) ,intent(out) :: enum_type
307 integer ,intent(out) :: Status
309 integer, dimension(2) :: data
314 call h5tenum_create_f(H5T_NATIVE_INTEGER,enum_type,hdf5err)
315 if(hdf5err.lt.0) then
316 Status = WRF_HDF5_ERR_DATATYPE
317 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
318 call wrf_debug ( WARN , msg)
322 call h5tenum_insert_f(enum_type,hdf5_true,data(1),hdf5err)
323 if(hdf5err.lt.0) then
324 Status = WRF_HDF5_ERR_DATATYPE
325 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
326 call wrf_debug ( WARN , msg)
330 call h5tenum_insert_f(enum_type,hdf5_false,data(2),Status)
331 if(hdf5err.lt.0) then
332 Status = WRF_HDF5_ERR_DATATYPE
333 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
334 call wrf_debug ( WARN , msg)
340 end subroutine SetUp_EnumID
342 ! Returns .TRUE. iff it is OK to write time-independent domain metadata to the
343 ! file referenced by DataHandle. If DataHandle is invalid, .FALSE. is
345 LOGICAL FUNCTION phdf5_ok_to_put_dom_ti( DataHandle )
347 include 'wrf_status_codes.h'
348 INTEGER, INTENT(IN) :: DataHandle
349 CHARACTER*80 :: fname
352 LOGICAL :: dryrun, first_output, retval
353 call ext_phdf5_inquire_filename( DataHandle, fname, filestate, Status )
354 IF ( Status /= WRF_NO_ERR ) THEN
355 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__, &
357 call wrf_debug ( WARN , TRIM(msg) )
360 dryrun = ( filestate .EQ. WRF_FILE_OPENED_NOT_COMMITTED )
361 first_output = phdf5_is_first_operation( DataHandle )
362 retval = .NOT. dryrun .AND. first_output
364 phdf5_ok_to_put_dom_ti = retval
366 END FUNCTION phdf5_ok_to_put_dom_ti
368 ! Returns .TRUE. iff it is OK to read time-independent domain metadata from the
369 ! file referenced by DataHandle. If DataHandle is invalid, .FALSE. is
371 LOGICAL FUNCTION phdf5_ok_to_get_dom_ti( DataHandle )
373 include 'wrf_status_codes.h'
374 INTEGER, INTENT(IN) :: DataHandle
375 CHARACTER*80 :: fname
378 LOGICAL :: dryrun, retval
379 call ext_phdf5_inquire_filename( DataHandle, fname, filestate, Status )
380 IF ( Status /= WRF_NO_ERR ) THEN
381 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__, &
383 call wrf_debug ( WARN , TRIM(msg) )
386 dryrun = ( filestate .EQ. WRF_FILE_OPENED_NOT_COMMITTED )
387 retval = .NOT. dryrun
389 phdf5_ok_to_get_dom_ti = retval
391 END FUNCTION phdf5_ok_to_get_dom_ti
393 ! Returns .TRUE. iff nothing has been read from or written to the file
394 ! referenced by DataHandle. If DataHandle is invalid, .FALSE. is returned.
395 LOGICAL FUNCTION phdf5_is_first_operation( DataHandle )
397 INCLUDE 'wrf_status_codes.h'
398 INTEGER, INTENT(IN) :: DataHandle
399 TYPE(wrf_phdf5_data_handle) ,POINTER :: DH
402 CALL GetDH( DataHandle, DH, Status )
403 IF ( Status /= WRF_NO_ERR ) THEN
404 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__, &
406 call wrf_debug ( WARN , TRIM(msg) )
409 retval = DH%first_operation
411 phdf5_is_first_operation = retval
413 END FUNCTION phdf5_is_first_operation
415 end module ext_phdf5_support_routines
417 !module wrf_phdf5_opt_data
418 ! integer ,parameter :: MaxOptVars = 100
419 !end module wrf_phdf5_opt_data
421 !module opt_data_module
423 !use wrf_phdf5_opt_data
427 ! integer,pointer :: darrays(:)
430 ! type(field),target :: fieldhandle(MaxOptVars)
431 !end module opt_data_module
433 !module opt_support_module
436 ! subroutine alloc_opt_handle(ODH)
437 ! use opt_data_module
438 ! type(field),pointer::DH
445 !end module opt_support_module
447 ! check the date, only use the length
448 subroutine DateCheck(Date,Status)
450 include 'wrf_status_codes.h'
451 character*(*) ,intent(in) :: Date
452 integer ,intent(out) :: Status
454 if(len(Date) /= DateStrLen) then
455 Status = WRF_HDF5_ERR_DATESTR_BAD_LENGTH
460 end subroutine DateCheck
462 ! This routine is for meta-data time dependent varible attribute
463 subroutine GetName(Element,Var,Name,Status)
466 include 'wrf_status_codes.h'
467 character*(*) ,intent(in) :: Element
468 character*(*) ,intent(in) :: Var
469 character*(*) ,intent(out) :: Name
470 integer ,intent(out) :: Status
471 character (VarNameLen) :: VarName
474 integer, parameter :: upper_to_lower =IACHAR('a')-IACHAR('A')
477 Name = 'MD___'//trim(Element)//VarName
480 if('A'<=c .and. c <='Z') Name(i:i)=achar(iachar(c)+upper_to_lower)
481 if(c=='-'.or.c==':') Name(i:i)='_'
485 end subroutine GetName
488 subroutine GetDataTimeIndex(IO,DataHandle,DateStr,TimeIndex,Status)
492 use ext_phdf5_support_routines
495 include 'wrf_status_codes.h'
497 character (*) ,intent(in) :: IO
498 integer ,intent(in) :: DataHandle
499 character*(*) ,intent(in) :: DateStr
500 character (DateStrLen), pointer :: TempTimes(:)
501 integer ,intent(out) :: TimeIndex
502 integer ,intent(out) :: Status
504 type(wrf_phdf5_data_handle) ,pointer :: DH
509 integer :: PreTimeCount
512 integer(hsize_t), dimension(1) :: chunk_dims =(/1/)
513 integer(hsize_t), dimension(1) :: dims
514 integer(hsize_t), dimension(1) :: hdf5_maxdims
515 integer(hsize_t), dimension(1) :: offset
516 integer(hsize_t), dimension(1) :: count
517 integer(hsize_t), dimension(1) :: sizes
519 INTEGER(HID_T) :: dset_id ! Dataset ID
520 INTEGER(HID_T) :: dspace_id ! Dataspace ID
521 INTEGER(HID_T) :: fspace_id ! Dataspace ID
522 INTEGER(HID_T) :: crp_list ! chunk ID
523 integer(hid_t) :: str_id ! string ID
526 integer(hid_t) :: group_id
527 character(Len = 512) :: groupname
531 character(len=100) :: buf
532 integer(size_t) :: name_size
533 integer(size_t) :: datelen_size
534 ! suppose the output will not exceed 100,0000 timesteps.
535 character(Len = MaxTimeSLen) :: tname
538 ! DH => WrfDataHandles(DataHandle), don't know why NetCDF doesn't use GetDH
539 call GetDH(DataHandle,DH,Status)
540 if(Status /= WRF_NO_ERR) then
541 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
542 call wrf_debug ( WARN , msg)
546 call DateCheck(DateStr,Status)
547 if(Status /= WRF_NO_ERR) then
548 Status = WRF_HDF5_ERR_DATESTR_ERROR
549 write(msg,*) 'Warning DATE STRING ERROR in ',"__FILE__",', line', __LINE__
550 call wrf_debug ( WARN , msg)
554 if(IO == 'write') then
555 TimeIndex = DH%TimeIndex
556 if(TimeIndex <= 0) then
558 elseif(DateStr < DH%Times(TimeIndex)) then
559 Status = WRF_HDF5_ERR_DATE_LT_LAST_DATE
560 write(msg,*) 'Warning DATE < LAST DATE in ',"__FILE__",', line', __LINE__
561 call wrf_debug ( WARN , msg)
563 elseif(DateStr == DH%Times(TimeIndex)) then
567 TimeIndex = TimeIndex + 1
568 ! If exceeding the maximum timestep, updating the maximum timestep
569 if(TimeIndex > MaxTimes*(DH%MaxTimeCount)) then
570 PreTimeCount = DH%MaxTimeCount
571 allocate(TempTimes(PreTimeCount*MaxTimes))
572 TempTimes(1:MaxTimes*PreTimeCount)=DH%Times(1:MaxTimes &
574 DH%MaxTimeCount = DH%MaxTimeCount +1
576 allocate(DH%Times(DH%MaxTimeCount*MaxTimes))
577 DH%Times(1:MaxTimes*PreTimeCount)=TempTimes(1:MaxTimes &
579 deallocate(TempTimes)
582 DH%TimeIndex = TimeIndex
583 DH%Times(TimeIndex) = DateStr
584 ! From NetCDF implementation, keep it in case it can be used.
586 ! VStart(2) = TimeIndex
587 ! VCount(1) = DateStrLen
590 ! create memory dataspace id and file dataspace id
593 offset(1) = TimeIndex -1
596 ! create group id for different time stamp
597 call numtochar(TimeIndex,tname)
598 groupname = 'TIME_STAMP_'//tname
599 ! call h5gn_members_f(DH%GroupID,DH%GroupName,nmembers,hdf5err)
600 ! do i = 0, nmembers - 1
601 ! call h5gget_obj_info_idx_f(DH%GroupID,DH%GroupName,i,ObjName, ObjType, &
604 ! if(ObjName(1:17) == groupname) then
605 ! call h5gopen_f(DH%GroupID,groupname,tgroupid,hdf5err)
610 if(DH%Tgroupids(TimeIndex) == -1) then
611 call h5gcreate_f(DH%groupid,groupname,group_id,hdf5err)
612 if(hdf5err .lt. 0) then
613 Status = WRF_HDF5_ERR_GROUP
614 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
615 call wrf_debug ( WARN , msg)
618 DH%Tgroupids(TimeIndex) = group_id
620 ! call h5gopen_f(DH%groupid,groupname,group_id,
621 group_id = DH%Tgroupids(TimeIndex)
624 call h5screate_simple_f(1,dims,dspace_id,hdf5err,dims)
625 if(hdf5err.lt.0) then
626 Status = WRF_HDF5_ERR_DATASPACE
627 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
628 call wrf_debug ( WARN , msg)
633 ! create HDF5 string handler for time
634 if(TimeIndex == 1) then
635 call h5tcopy_f(H5T_NATIVE_CHARACTER, str_id, hdf5err)
636 if(hdf5err.lt.0) then
637 Status = WRF_HDF5_ERR_DATATYPE
638 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
639 call wrf_debug ( WARN , msg)
643 datelen_size = DateStrLen
644 call h5tset_size_f(str_id,datelen_size,hdf5err)
645 if(hdf5err.lt.0) then
646 Status = WRF_HDF5_ERR_DATATYPE
647 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
648 call wrf_debug ( WARN , msg)
655 call h5dcreate_f(group_id,DH%TimesName,str_id,dspace_id,&
657 if(hdf5err.lt.0) then
658 Status = WRF_HDF5_ERR_DATASET_CREATE
659 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
660 call wrf_debug ( WARN , msg)
665 ! write the data in memory space to file space
666 CALL h5dwrite_f(DH%TimesID,str_id,DateStr,dims,hdf5err,dspace_id,dspace_id)
667 if(hdf5err.lt.0) then
668 Status = WRF_HDF5_ERR_DATASET_WRITE
669 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
670 call wrf_debug ( WARN , msg)
674 if(TimeIndex == 1) then
679 call h5sclose_f(dspace_id,hdf5err)
680 if(hdf5err.lt.0) then
681 Status = WRF_HDF5_ERR_DATASPACE
682 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
683 call wrf_debug ( WARN , msg)
687 call h5dclose_f(DH%TimesID,hdf5err)
688 if(hdf5err.lt.0) then
689 Status = WRF_HDF5_ERR_DATASET_GENERAL
690 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
691 call wrf_debug ( WARN , msg)
696 ! This is for IO read
697 ! Find the timeIndex(very expensive for large
698 ! time stamp, should use hashing table)
700 do i=1,MaxTimes*DH%MaxTimeCount
702 ! For handling reading maximum timestamp greater than 9000 in the future
703 ! if(DH%Times(i) == NullName) then
704 ! Status = WRF_HDF5_ERR_TIME
705 ! write(msg,*) 'Warning TIME ',DateStr,' NOT FOUND in ',"__FILE__",&
707 ! call wrf_debug ( WARN , msg)
711 if(DH%Times(i) == DateStr) then
717 ! Need a recursive function to handle this
718 ! This is a potential bug
719 if(i == MaxTimes*DH%MaxTimeCount) then
720 ! PreTimeCount = DH%MaxTimeCount
721 ! allocate(TempTimes(PreTimeCount*MaxTimes))
722 ! TempTimes(1:MaxTimes*PreTimeCount)=DH%Times(1:MaxTimes &
724 ! DH%MaxTimeCount = DH%MaxTimeCount +1
725 ! deallocate(DH%Times)
726 ! allocate(DH%Times(DH%MaxTimeCount*MaxTimes))
727 ! DH%Times(1:MaxTimes*PreTimeCount)=TempTimes(1:MaxTimes &
729 ! deallocate(TempTimes)
730 Status = WRF_HDF5_ERR_TIME
731 write(msg,*) 'Warning TIME ',DateStr,' NOT FOUND in ',"__FILE__",&
733 call wrf_debug ( WARN , msg)
738 ! do the hyperslab selection
742 end subroutine GetDataTimeIndex
744 subroutine GetAttrTimeIndex(IO,DataHandle,DateStr,TimeIndex,Status)
748 use ext_phdf5_support_routines
751 include 'wrf_status_codes.h'
753 character (*) ,intent(in) :: IO
754 integer ,intent(in) :: DataHandle
755 character*(*) ,intent(in) :: DateStr
756 character (DateStrLen), pointer :: TempTimes(:)
757 integer ,intent(out) :: TimeIndex
758 integer ,intent(out) :: Status
760 type(wrf_phdf5_data_handle) ,pointer :: DH
765 integer :: PreTimeCount
768 integer(hsize_t), dimension(1) :: chunk_dims =(/1/)
769 integer(hsize_t), dimension(1) :: dims
770 integer(hsize_t), dimension(1) :: hdf5_maxdims
771 integer(hsize_t), dimension(1) :: offset
772 integer(hsize_t), dimension(1) :: count
773 integer(hsize_t), dimension(1) :: sizes
775 INTEGER(HID_T) :: dset_id ! Dataset ID
776 INTEGER(HID_T) :: dspace_id ! Dataspace ID
777 INTEGER(HID_T) :: fspace_id ! Dataspace ID
778 INTEGER(HID_T) :: crp_list ! chunk ID
779 integer(hid_t) :: str_id ! string ID
782 integer(size_t) :: datelen_size
783 integer(hid_t) :: group_id
784 character(Len = 512) :: groupname
786 ! suppose the output will not exceed 100,0000 timesteps.
787 character(Len = MaxTimeSLen) :: tname
789 ! DH => WrfDataHandles(DataHandle), don't know why NetCDF doesn't use GetDH
790 call GetDH(DataHandle,DH,Status)
791 if(Status /= WRF_NO_ERR) then
792 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
793 call wrf_debug ( WARN , msg)
797 call DateCheck(DateStr,Status)
798 if(Status /= WRF_NO_ERR) then
799 Status = WRF_HDF5_ERR_DATESTR_ERROR
800 write(msg,*) 'Warning DATE STRING ERROR in ',"__FILE__",', line', __LINE__
801 call wrf_debug ( WARN , msg)
805 if(IO == 'write') then
806 TimeIndex = DH%TimeIndex
807 if(TimeIndex <= 0) then
809 elseif(DateStr < DH%Times(TimeIndex)) then
810 Status = WRF_HDF5_ERR_DATE_LT_LAST_DATE
811 write(msg,*) 'Warning DATE < LAST DATE in ',"__FILE__",', line', __LINE__
812 call wrf_debug ( WARN , msg)
814 elseif(DateStr == DH%Times(TimeIndex)) then
818 TimeIndex = TimeIndex + 1
819 ! If exceeding the maximum timestep, updating the maximum timestep
820 if(TimeIndex > MaxTimes*(DH%MaxTimeCount)) then
821 PreTimeCount = DH%MaxTimeCount
822 allocate(TempTimes(PreTimeCount*MaxTimes))
823 TempTimes(1:MaxTimes*PreTimeCount)=DH%Times(1:MaxTimes &
825 DH%MaxTimeCount = DH%MaxTimeCount +1
827 allocate(DH%Times(DH%MaxTimeCount*MaxTimes))
828 DH%Times(1:MaxTimes*PreTimeCount)=TempTimes(1:MaxTimes &
830 deallocate(TempTimes)
833 DH%TimeIndex = TimeIndex
834 DH%Times(TimeIndex) = DateStr
836 ! From NetCDF implementation, keep it in case it can be used.
838 ! VStart(2) = TimeIndex
839 ! VCount(1) = DateStrLen
842 ! create memory dataspace id and file dataspace id
845 offset(1) = TimeIndex -1
848 ! create group id for different time stamp
849 call numtochar(TimeIndex,tname)
850 groupname = 'TIME_STAMP_'//tname
852 if(DH%Tgroupids(TimeIndex) == -1) then
853 call h5gcreate_f(DH%groupid,groupname,group_id,hdf5err)
854 if(hdf5err .lt. 0) then
855 Status = WRF_HDF5_ERR_GROUP
856 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
857 call wrf_debug ( WARN , msg)
860 DH%Tgroupids(TimeIndex) = group_id
862 ! call h5gopen_f(DH%groupid,groupname,group_id,
863 group_id = DH%Tgroupids(TimeIndex)
866 call h5screate_simple_f(1,dims,dspace_id,hdf5err,dims)
867 if(hdf5err.lt.0) then
868 Status = WRF_HDF5_ERR_DATASPACE
869 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
870 call wrf_debug ( WARN , msg)
875 ! create HDF5 string handler for time
876 if(TimeIndex == 1) then
877 call h5tcopy_f(H5T_NATIVE_CHARACTER, str_id, hdf5err)
878 if(hdf5err.lt.0) then
879 Status = WRF_HDF5_ERR_DATATYPE
880 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
881 call wrf_debug ( WARN , msg)
885 datelen_size = DateStrLen
886 call h5tset_size_f(str_id,datelen_size,hdf5err)
887 if(hdf5err.lt.0) then
888 Status = WRF_HDF5_ERR_DATATYPE
889 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
890 call wrf_debug ( WARN , msg)
897 call h5dcreate_f(group_id,DH%TimesName,str_id,dspace_id,&
899 if(hdf5err.lt.0) then
900 Status = WRF_HDF5_ERR_DATASET_CREATE
901 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
902 call wrf_debug ( WARN , msg)
907 ! write the data in memory space to file space
908 CALL h5dwrite_f(DH%TimesID,str_id,DateStr,dims,hdf5err,dspace_id,dspace_id)
909 if(hdf5err.lt.0) then
910 Status = WRF_HDF5_ERR_DATASET_WRITE
911 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
912 call wrf_debug ( WARN , msg)
916 if(TimeIndex == 1) then
921 call h5sclose_f(dspace_id,hdf5err)
922 if(hdf5err.lt.0) then
923 Status = WRF_HDF5_ERR_DATASPACE
924 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
925 call wrf_debug ( WARN , msg)
929 call h5dclose_f(DH%TimesID,hdf5err)
930 if(hdf5err.lt.0) then
931 Status = WRF_HDF5_ERR_DATASET_GENERAL
932 write(msg,*) 'Warning Status = ',Status,' in ',__FILE__,', line', __LINE__
933 call wrf_debug ( WARN , msg)
938 ! This is for IO read
939 ! Find the timeIndex(very expensive for large
940 ! time stamp, should use hashing table)
942 do i=1,MaxTimes*DH%MaxTimeCount
945 if(DH%Times(i) == DateStr) then
951 ! Need a recursive function to handle this
952 ! This is a potential bug
953 if(i == MaxTimes*DH%MaxTimeCount) then
954 Status = WRF_HDF5_ERR_TIME
955 write(msg,*) 'Warning TIME ',DateStr,' NOT FOUND in ',"__FILE__",&
957 call wrf_debug ( WARN , msg)
962 ! do the hyperslab selection
966 end subroutine GetAttrTimeIndex
969 ! Obtain the rank of the dimension
970 subroutine GetDim(MemoryOrder,NDim,Status)
972 include 'wrf_status_codes.h'
973 character*(*) ,intent(in) :: MemoryOrder
974 integer ,intent(out) :: NDim
975 integer ,intent(out) :: Status
976 character*3 :: MemOrd
978 call LowerCase(MemoryOrder,MemOrd)
980 case ('xyz','xzy','yxz','yzx','zxy','zyx','xsz','xez','ysz','yez')
982 case ('xy','yx','xs','xe','ys','ye')
987 Status = WRF_HDF5_ERR_BAD_MEMORYORDER
992 end subroutine GetDim
994 ! Obtain the index for transposing
995 subroutine GetIndices(NDim,Start,End,i1,i2,j1,j2,k1,k2)
996 integer ,intent(in) :: NDim
997 integer ,dimension(*),intent(in) :: Start,End
998 integer ,intent(out) :: i1,i2,j1,j2,k1,k2
1008 if(NDim == 1) return
1011 if(NDim == 2) return
1015 end subroutine GetIndices
1017 ! shuffling the memory order to XYZ order
1018 subroutine ExtOrder(MemoryOrder,Vector,Status)
1020 include 'wrf_status_codes.h'
1021 character*(*) ,intent(in) :: MemoryOrder
1022 integer,dimension(*) ,intent(inout) :: Vector
1023 integer ,intent(out) :: Status
1025 integer,dimension(NVarDims) :: temp
1026 character*3 :: MemOrd
1028 call GetDim(MemoryOrder,NDim,Status)
1029 temp(1:NDim) = Vector(1:NDim)
1030 call LowerCase(MemoryOrder,MemOrd)
1031 select case (MemOrd)
1033 case ('xyz','xsz','xez','ysz','yez','xy','xs','xe','ys','ye','z','c')
1058 Status = WRF_HDF5_ERR_BAD_MEMORYORDER
1063 end subroutine ExtOrder
1065 ! shuffling the dimensional name order
1066 subroutine ExtOrderStr(MemoryOrder,Vector,ROVector,Status)
1068 include 'wrf_status_codes.h'
1069 character*(*) ,intent(in) :: MemoryOrder
1070 character*(*),dimension(*) ,intent(in) :: Vector
1071 character(256),dimension(NVarDims),intent(out) :: ROVector
1072 integer ,intent(out) :: Status
1074 character*3 :: MemOrd
1076 call GetDim(MemoryOrder,NDim,Status)
1077 ROVector(1:NDim) = Vector(1:NDim)
1078 call LowerCase(MemoryOrder,MemOrd)
1079 select case (MemOrd)
1081 case ('xyz','xsz','xez','ysz','yez','xy','xs','xe','ys','ye','z','c')
1084 ROVector(1) = 'ext_scalar'
1086 ROVector(2) = Vector(3)
1087 ROVector(3) = Vector(2)
1089 ROVector(1) = Vector(2)
1090 ROVector(2) = Vector(1)
1092 ROVector(1) = Vector(3)
1093 ROVector(2) = Vector(1)
1094 ROVector(3) = Vector(2)
1096 ROVector(1) = Vector(2)
1097 ROVector(2) = Vector(3)
1098 ROVector(3) = Vector(1)
1100 ROVector(1) = Vector(3)
1101 ROVector(3) = Vector(1)
1103 ROVector(1) = Vector(2)
1104 ROVector(2) = Vector(1)
1106 Status = WRF_HDF5_ERR_BAD_MEMORYORDER
1111 end subroutine ExtOrderStr
1114 subroutine LowerCase(MemoryOrder,MemOrd)
1115 character*(*) ,intent(in) :: MemoryOrder
1116 character*(*) ,intent(out) :: MemOrd
1118 integer ,parameter :: upper_to_lower =IACHAR('a')-IACHAR('A')
1122 N = len(MemoryOrder)
1123 MemOrd(1:N) = MemoryOrder(1:N)
1125 c = MemoryOrder(i:i)
1126 if('A'<=c .and. c <='Z') MemOrd(i:i)=achar(iachar(c)+upper_to_lower)
1129 end subroutine LowerCase
1131 subroutine UpperCase(MemoryOrder,MemOrd)
1132 character*(*) ,intent(in) :: MemoryOrder
1133 character*(*) ,intent(out) :: MemOrd
1135 integer ,parameter :: lower_to_upper =IACHAR('A')-IACHAR('a')
1139 N = len(MemoryOrder)
1140 MemOrd(1:N) = MemoryOrder(1:N)
1142 c = MemoryOrder(i:i)
1143 if('a'<=c .and. c <='z') MemOrd(i:i)=achar(iachar(c)+lower_to_upper)
1146 end subroutine UpperCase
1148 ! subroutine used in transpose routine
1149 subroutine reorder (MemoryOrder,MemO)
1150 character*(*) ,intent(in) :: MemoryOrder
1151 character*3 ,intent(out) :: MemO
1152 character*3 :: MemOrd
1153 integer :: N,i,i1,i2,i3
1156 N = len_trim(MemoryOrder)
1158 call lowercase(MemoryOrder,MemOrd)
1162 if(ichar(MemOrd(i:i)) < ichar(MemOrd(i1:i1))) I1 = i
1163 if(ichar(MemOrd(i:i)) > ichar(MemOrd(i3:i3))) I3 = i
1170 MemO(1:1) = MemoryOrder(i1:i1)
1171 MemO(2:2) = MemoryOrder(i2:i2)
1172 if(N == 3) MemO(3:3) = MemoryOrder(i3:i3)
1173 if(MemOrd(i1:i1) == 's' .or. MemOrd(i1:i1) == 'e') then
1174 MemO(1:N-1) = MemO(2:N)
1175 MemO(N:N ) = MemoryOrder(i1:i1)
1178 end subroutine reorder
1180 subroutine Transpose_hdf5(IO,MemoryOrder,di, Field,l1,l2,m1,m2,n1,n2 &
1181 ,XField,x1,x2,y1,y2,z1,z2 &
1182 ,i1,i2,j1,j2,k1,k2 )
1183 character*(*) ,intent(in) :: IO
1184 character*(*) ,intent(in) :: MemoryOrder
1185 integer ,intent(in) :: l1,l2,m1,m2,n1,n2
1186 integer ,intent(in) :: di
1187 integer ,intent(in) :: x1,x2,y1,y2,z1,z2
1188 integer ,intent(in) :: i1,i2,j1,j2,k1,k2
1189 integer ,intent(inout) :: Field(di,l1:l2,m1:m2,n1:n2)
1190 !jm 010827 integer ,intent(inout) :: XField(di,x1:x2,y1:y2,z1:z2)
1191 integer ,intent(inout) :: XField(di,(i2-i1+1)*(j2-j1+1)*(k2-k1+1))
1192 character*3 :: MemOrd
1194 integer ,parameter :: MaxUpperCase=IACHAR('Z')
1195 integer :: i,j,k,ix,jx,kx
1197 call LowerCase(MemoryOrder,MemOrd)
1198 select case (MemOrd)
1200 !#define XDEX(A,B,C) A-A ## 1+1+(A ## 2-A ## 1+1)*((B-B ## 1)+(C-C ## 1)*(B ## 2-B ## 1+1))
1208 call reorder(MemoryOrder,MemO)
1209 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1210 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1211 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1215 if(IO == 'write') then
1216 XField(1:di,(i-i1+1+(i2-i1+1)*((k-k1)+(j-j1)*(k2-k1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1218 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(i-i1+1+(i2-i1+1)*((k-k1)+(j-j1)*(k2-k1+1))))
1225 case ('xyz','xsz','xez','ysz','yez','xy','xs','xe','ys','ye','z','c','0')
1231 call reorder(MemoryOrder,MemO)
1232 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1233 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1234 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1238 if(IO == 'write') then
1239 XField(1:di,(i-i1+1+(i2-i1+1)*((j-j1)+(k-k1)*(j2-j1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1241 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(i-i1+1+(i2-i1+1)*((j-j1)+(k-k1)*(j2-j1+1))))
1254 call reorder(MemoryOrder,MemO)
1255 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1256 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1257 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1261 if(IO == 'write') then
1262 XField(1:di,(j-j1+1+(j2-j1+1)*((i-i1)+(k-k1)*(i2-i1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1264 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(j-j1+1+(j2-j1+1)*((i-i1)+(k-k1)*(i2-i1+1))))
1277 call reorder(MemoryOrder,MemO)
1278 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1279 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1280 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1284 if(IO == 'write') then
1285 XField(1:di,(k-k1+1+(k2-k1+1)*((i-i1)+(j-j1)*(i2-i1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1287 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(k-k1+1+(k2-k1+1)*((i-i1)+(j-j1)*(i2-i1+1))))
1300 call reorder(MemoryOrder,MemO)
1301 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1302 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1303 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1307 if(IO == 'write') then
1308 XField(1:di,(j-j1+1+(j2-j1+1)*((k-k1)+(i-i1)*(k2-k1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1310 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(j-j1+1+(j2-j1+1)*((k-k1)+(i-i1)*(k2-k1+1))))
1323 call reorder(MemoryOrder,MemO)
1324 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1325 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1326 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1330 if(IO == 'write') then
1331 XField(1:di,(k-k1+1+(k2-k1+1)*((j-j1)+(i-i1)*(j2-j1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1333 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(k-k1+1+(k2-k1+1)*((j-j1)+(i-i1)*(j2-j1+1))))
1346 call reorder(MemoryOrder,MemO)
1347 if(IACHAR(MemO(1:1)) > MaxUpperCase) ix=i2+i1
1348 if(IACHAR(MemO(2:2)) > MaxUpperCase) jx=j2+j1
1349 if(IACHAR(MemO(3:3)) > MaxUpperCase) kx=k2+k1
1353 if(IO == 'write') then
1354 XField(1:di,(j-j1+1+(j2-j1+1)*((i-i1)+(k-k1)*(i2-i1+1)))) = Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k))
1356 Field(1:di,abs(ix-i),abs(jx-j),abs(kx-k)) = XField(1:di,(j-j1+1+(j2-j1+1)*((i-i1)+(k-k1)*(i2-i1+1))))
1365 end subroutine Transpose_hdf5
1367 subroutine numtochar(TimeIndex,tname,Status)
1370 integer, intent(in) :: TimeIndex
1371 character(len=MaxTimeSLen),intent(out)::tname
1372 integer ,intent(out)::Status
1373 integer :: i,ten_pow,temp
1374 integer :: maxtimestep
1378 maxtimestep = maxtimestep * 10
1380 if(TimeIndex >= maxtimestep) then
1381 Status = WRF_HDF5_ERR_OTHERS
1382 write(msg,*) 'Cannot exceed the maximum timestep',maxtimestep,'in',__FILE__,' line',__LINE__
1383 call wrf_debug(FATAL,msg)
1390 tname(MaxTimeSLen+1-i:MaxTimeSLen+1-i) = achar(modulo(TimeIndex/ten_pow,temp)+iachar('0'))
1391 ten_pow = 10* ten_pow
1395 end subroutine numtochar