diff --git a/CHANGELOG.md b/CHANGELOG.md index a6b309c9ffa0..ec3166efe70d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added subroutine to read nc4 tile file - Added optional start_date and start_time to control the output window for each History collection. No output will be written before then. If not specified, these default to the beginning of the experiment. - Added utility to prepare inputs for ExtDatDriver.x so that ExtData can simulate a real GEOS run - Added loggers when writing or reading weight files diff --git a/base/MAPL_LocStreamMod.F90 b/base/MAPL_LocStreamMod.F90 index d0f654f73900..37a8963cc2d2 100644 --- a/base/MAPL_LocStreamMod.F90 +++ b/base/MAPL_LocStreamMod.F90 @@ -69,8 +69,8 @@ module MAPL_LocStreamMod !EOP -integer, parameter :: NumGlobalVars=4 -integer, parameter :: NumLocalVars =4 +!integer, parameter :: NumGlobalVars=4 +!integer, parameter :: NumLocalVars =4 type MAPL_GeoLocation @@ -343,18 +343,19 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, integer :: UNIT integer :: N, I, K, L, NT type(MAPL_LocStreamType), pointer :: STREAM - real, pointer :: AVR(:,:), AVR_transpose(:,:) + real, pointer :: AVR(:,:) logical, pointer :: MSK(:) real :: X, Y, X0, Y0, XE, DX, DY integer :: II, JJ logical :: DoCoeffs character(len=MAPL_TileNameLength):: gname - + character(len=MAPL_TileNameLength):: GridNames(2) + integer :: IMs(2), JMs(2) integer :: irec integer(kind=1) :: byte(4) - integer :: I1, IN, J1, JN - integer :: iostat - logical :: isascii + integer :: I1, IN, J1, JN, N_Grids, N_PfafCat + integer :: iostat, filetype + logical :: isascii, isnc4, isbinary logical :: read_always logical, pointer :: ISMINE(:) type(MAPL_Tiling ), pointer :: TILING @@ -404,94 +405,63 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, ! Use some heuristics to determine filetype (choices are BINARY and ASCII) !------------------------------------------------------------------------ - ! just get the unit - UNIT = GETFILE(FILENAME, DO_OPEN=0, ALL_PES=.true., RC=STATUS) - _VERIFY(STATUS) + call MAPL_NCIOGetFileType(FILENAME, filetype, _RC) - INQUIRE(IOLENGTH=IREC) BYTE - open (UNIT=UNIT, FILE=FILENAME, FORM='unformatted', ACCESS='DIRECT', RECL=IREC, IOSTAT=status) - _VERIFY(STATUS) - read (UNIT, REC=1, ERR=100) BYTE - call FREE_FILE(UNIT) - - isascii = .true. - do i=1,size(byte) - if (BYTE(I) < 7) then - isascii = .false. - exit - end if - end do + isnc4 = (filetype == MAPL_FILETYPE_NC4) + isascii = (filetype == MAPL_FILETYPE_TXT) + isbinary= (filetype == MAPL_FILETYPE_BIN) - if (isascii) then + if ( .not. isbinary) then ! Open file and read header info !------------------------------- - UNIT = GETFILE(FILENAME, form='FORMATTED', RC=status) - _VERIFY(STATUS) - -! Total number of tiles in exchange grid -!--------------------------------------- + if (isnc4) then + if (MAPL_AM_I_root()) then + call MAPL_ReadTilingNC4(FILENAME, GridName=GridNames, IM=IMs, JM=JMs, N_Grids=N_Grids, N_PfafCat=N_PfafCat, AVR=AVR, _RC) + NT = size(AVR,1) + endif + call MAPL_CommsBcast(layout, NT, 1, MAPL_Root, status) + call MAPL_CommsBcast(layout, N_Grids, 1, MAPL_Root, status) + call MAPL_CommsBcast(layout, IMs, 2, MAPL_Root, status) + call MAPL_CommsBcast(layout, JMs, 2, MAPL_Root, status) - if (use_pfaf_) then - call READ_PARALLEL(layout, hdr, UNIT=UNIT, rc=status) - _VERIFY(STATUS) - nt=hdr(1) - stream%pfafstetter_catchments=hdr(2) - else - call READ_PARALLEL(layout, nt, UNIT=UNIT, rc=status) - _VERIFY(STATUS) - end if + if (use_pfaf_) then + call MAPL_CommsBcast(layout, N_PfafCat, 1, MAPL_Root, status) + endif -! Number of grids that can be attached -!------------------------------------- + do N = 1, N_Grids + call MAPL_CommsBcast(layout, GridNames(N), MAPL_TileNameLength, MAPL_Root, status) + enddo - call READ_PARALLEL(layout, STREAM%N_GRIDS, unit=UNIT, rc=status) - _VERIFY(STATUS) + if (.not. MAPL_AM_I_root()) then + allocate(AVR(NT, NumGlobalVars+NumLocalVars*N_GRIDS)) + endif + call MAPL_CommsBcast(layout, AVR, NT*(NumGlobalVars+NumLocalVars*N_GRIDS), MAPL_Root, status) + endif -! The exchange grid is used to tile each attached grid -!----------------------------------------------------- + if (isascii) then + call MAPL_ReadTilingASCII(layout, FILENAME, GridNames, NT, IMs, JMs, N_Grids, N_PfafCat, AVR, _RC) + endif + STREAM%N_GRIDS = N_Grids allocate(STREAM%TILING(STREAM%N_GRIDS), STAT=STATUS) _VERIFY(STATUS) -! The names and sizes of the grids to be tiled -!--------------------------------------------- - - do N=1,STREAM%N_GRIDS - call READ_PARALLEL(layout, STREAM%TILING(N)%NAME, unit=UNIT, rc=status) - _VERIFY(STATUS) + do N = 1, N_Grids + STREAM%TILING(N)%NAME = GridNames(N) if (NewGridNames_) then - call GenOldGridName_(STREAM%TILING(N)%NAME) + call GenOldGridName_(STREAM%TILING(N)%NAME) end if - call READ_PARALLEL(layout, STREAM%TILING(N)%IM, unit=UNIT, rc=status) - _VERIFY(STATUS) - call READ_PARALLEL(layout, STREAM%TILING(N)%JM, unit=UNIT, rc=status) - _VERIFY(STATUS) + STREAM%TILING(N)%IM = IMs(N) + STREAM%TILING(N)%JM = JMs(N) enddo if (use_pfaf_) then + stream%pfafstetter_catchments = N_PfafCat STREAM%TILING(2)%IM = stream%pfafstetter_catchments STREAM%TILING(2)%JM = 1 STREAM%TILING(2)%name = "CATCHMENT_GRID" end if - -! Read location stream file into AVR -!--------------------------------------- - if(index(STREAM%TILING(1)%NAME,'EASE') /=0 ) then - allocate(AVR(NT,9), STAT=STATUS) ! 9 columns for EASE grid - _VERIFY(STATUS) - allocate(AVR_transpose(9,NT)) - else - allocate(AVR(NT,NumGlobalVars+NumLocalVars*STREAM%N_GRIDS), STAT=STATUS) - _VERIFY(STATUS) - allocate(AVR_transpose(NumGlobalVars+NumLocalVars*STREAM%N_GRIDS,NT), STAT=STATUS) - _VERIFY(STATUS) - endif - - call READ_PARALLEL(layout, AVR_transpose(:,:), unit=UNIT, rc=status) - AVR= transpose(AVR_transpose) - deallocate(AVR_transpose) - ! adjust EASE grid starting index. Internally, the starting index is 1 instead of 0. do N=1,STREAM%N_GRIDS if(index(STREAM%TILING(N)%NAME,'EASE') /=0 ) then @@ -499,12 +469,11 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, AVR(:,NumGlobalVars+2+NumLocalVars*(N-1)) = AVR(:,NumGlobalVars+2+NumLocalVars*(N-1))+1 endif enddo - + !if (use_pfaf_) then !AVR(:,NumGlobalVars+2+NumLocalVars) = 1 !end if - call FREE_FILE(UNIT) ! Allocate msk for which tiles to include in the stream being created. !-------------------------------------------------------------------- @@ -880,7 +849,7 @@ subroutine MAPL_LocStreamCreateFromFile(LocStream, LAYOUT, FILENAME, NAME, MASK, ! and are expensive to reread. !------------------------------------------------------------------ - if(.not.isascii) then + if( isbinary ) then if ( MAPL_am_I_root() ) then rewind(UNIT) diff --git a/base/NCIO.F90 b/base/NCIO.F90 index e638d662a4b1..f6657a582ff1 100644 --- a/base/NCIO.F90 +++ b/base/NCIO.F90 @@ -22,6 +22,8 @@ module NCIOMod use MAPL_ExceptionHandling use netcdf use pFIO + use BinIOMod, only: GETFILE, READ_PARALLEL, FREE_FILE + use MAPL_Constants !use pFIO_ClientManagerMod use gFTL_StringIntegerMap use gFTL_StringVector @@ -43,6 +45,8 @@ module NCIOMod public MAPL_NCIOGetFileType public MAPL_VarReadNCPar public MAPL_VarWriteNCPar + public MAPL_ReadTilingNC4 + public MAPL_ReadTilingASCII interface MAPL_VarReadNCPar module procedure MAPL_StateVarReadNCPar @@ -78,7 +82,10 @@ module NCIOMod module procedure MAPL_VarWriteNCpar_R8_4d end interface - contains + integer, parameter, public :: NumGlobalVars=4 + integer, parameter, public :: NumLocalVars =4 + +contains subroutine MAPL_FieldReadNCPar(formatter,name,FIELD, ARRDES, HomePE, RC) @@ -4509,8 +4516,12 @@ subroutine MAPL_StateVarWriteNCPar(filename, STATE, ARRDES, CLOCK, NAME, forceWr end subroutine MAPL_StateVarWriteNCPar subroutine MAPL_NCIOGetFileType(filename,filetype,rc) - implicit none + ! filetype = 0, hdf5 + ! filetype = 1, ascii + ! filetype = 2, binary or unknown + ! filetype = -1, unknown + implicit none ! Arguments !---------- character(len=*), intent(IN ) :: filename @@ -4528,7 +4539,7 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) integer :: irec integer :: unit integer :: i, cwrd - logical :: typehdf5 + logical :: typehdf5, isascii character(len=12) :: fmt INQUIRE(IOLENGTH=IREC) WORD @@ -4540,9 +4551,9 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) read (UNIT, REC=2, ERR=100) TwoWords(5:8) close(UNIT) - typehdf5 = .true. - filetype = -1 ! Unknown + filetype = MAPL_FILETYPE_UNK ! Unknown + typehdf5 = .true. do i = 1, 8 if (iachar(TwoWords(i)) /= hdf5(i)) then typehdf5 = .false. @@ -4550,20 +4561,23 @@ subroutine MAPL_NCIOGetFileType(filename,filetype,rc) end if end do if (typehdf5) then - filetype = 0 ! HDF5 + filetype = MAPL_FILETYPE_NC4 _RETURN(ESMF_SUCCESS) - - end if - ! Attempt to identify as fortran binary - cwrd = transfer(TwoWords(1:4), irec) - ! check if divisible by 4 - irec = cwrd/4 - filetype = irec - if (cwrd /= 4*irec) then - _RETURN(ESMF_FAILURE) end if - filetype = -1 + isascii = .true. + do i = 1, 4 + if (iachar(TwoWords(i)) < 7) then + isascii = .false. + exit + end if + end do + if (isascii) then + filetype = MAPL_FILETYPE_TXT + _RETURN(ESMF_SUCCESS) + endif + + filetype = MAPL_FILETYPE_BIN _RETURN(ESMF_SUCCESS) 100 continue @@ -5064,5 +5078,241 @@ function create_flipped_field(field,rc) result(flipped_field) _RETURN(_SUCCESS) end function create_flipped_field + subroutine MAPL_ReadTilingNC4(File, GridName, im, jm, nx, ny, n_Grids, iTable, rTable, N_PfafCat, AVR,rc) + character(*), intent(IN) :: File + character(*), optional, intent(out) :: GridName(:) + integer, optional, intent(out) :: IM(:), JM(:) + integer, optional, intent(out) :: nx, ny, n_Grids + integer, optional, allocatable, intent(out) :: iTable(:,:) + real(kind=8), optional, allocatable, intent(out) :: rTable(:,:) + integer, optional, intent(out) :: N_PfafCat + real, optional, pointer, intent(out) :: AVR(:,:) ! used by GEOSgcm + integer, optional, intent(out) :: rc + + type (Attribute), pointer :: ref + character(len=:), allocatable :: attr + type (NetCDF4_FileFormatter) :: formatter + type (FileMetadata) :: meta + character(len=1) :: str_num + integer :: ng, ntile, status, ll + class(*), pointer :: attr_val(:) + class(*), pointer :: char_val + integer, allocatable :: tmp_int(:) + real(kind=8), allocatable :: fr(:) + + integer :: NumCol + integer, allocatable :: iTable_(:,:) + real(kind=8), allocatable :: rTable_(:,:) + + call formatter%open(File, pFIO_READ, rc=status) + meta = formatter%read(rc=status) + + ng = 1 + if (meta%has_attribute("Grid1_Name")) ng = 2 ! for ng=1 (i.e., EASE), attribute is just "Grid_Name" + ntile = meta%get_dimension('tile') + + if (present(n_Grids)) then + n_Grids = ng + endif + + if (present(nx)) then + ref => meta%get_attribute('raster_nx') + attr_val => ref%get_values() + select type(attr_val) + type is (integer(INT32)) + nx = attr_val(1) + endselect + endif + if (present(ny)) then + ref => meta%get_attribute('raster_ny') + attr_val => ref%get_values() + select type (attr_val) + type is (integer(INT32)) + ny = attr_val(1) + endselect + endif + + if (present(N_PfafCat)) then + ref => meta%get_attribute('N_PfafCat') + attr_val => ref%get_values() + select type (attr_val) + type is (integer(INT32)) + N_PfafCat = attr_val(1) + endselect + endif + do ll = 1, ng + if (ng == 1) then + str_num = '' + else + write(str_num, '(i0)') ll + endif + + if (present(GridName)) then + attr = 'Grid'//trim(str_num)//'_Name' + ref =>meta%get_attribute(attr) + char_val => ref%get_value() + select type(char_val) + type is(character(*)) + GridName(ll) = char_val + class default + print('unsupported subclass (not string) of attribute named '//attr) + end select + endif + if (present(IM)) then + attr = 'IM'//trim(str_num) + ref =>meta%get_attribute(attr) + attr_val => ref%get_values() + select type(attr_val) + type is( integer(INT32)) + IM(ll) = attr_val(1) + end select + endif + if (present(JM)) then + attr = 'JM'//trim(str_num) + ref =>meta%get_attribute(attr) + attr_val => ref%get_values() + select type(attr_val) + type is(integer(INT32)) + JM(ll) = attr_val(1) + end select + endif + enddo + + if (present(iTable) .or. present(AVR) ) then + allocate(iTable_(ntile,0:7)) + allocate(tmp_int(ntile)) + call formatter%get_var('typ', iTable_(:,0)) + do ll = 1, ng + if (ng == 1) then + str_num='' + else + write(str_num, '(i0)') ll + endif + + call formatter%get_var('i_indg' //trim(str_num), tmp_int, rc=status) + iTable_(:,ll*2) = tmp_int + call formatter%get_var('j_indg' //trim(str_num), tmp_int, rc=status) + iTable_(:,ll*2+1) = tmp_int + call formatter%get_var('pfaf_index'//trim(str_num), tmp_int, rc=status) + if ( ng == 1) then + iTable_(:,4) = tmp_int + else + iTable_(:,5+ll) = tmp_int + endif + enddo + endif + if (present(rTable) .or. present(AVR) ) then + allocate(rTable_(ntile,10)) + call formatter%get_var('com_lon', rTable_(:,1), rc=status) + call formatter%get_var('com_lat', rTable_(:,2), rc=status) + call formatter%get_var('area', rTable_(:,3), rc=status) + do ll = 1, ng + if (ng == 1) then + str_num='' + else + write(str_num, '(i0)') ll + endif + call formatter%get_var('frac_cell' //trim(str_num), rTable_(:,3+ll), rc=status) + enddo + call formatter%get_var('min_lon', rTable_(:, 6), rc=status) + call formatter%get_var('max_lon', rTable_(:, 7), rc=status) + call formatter%get_var('min_lat', rTable_(:, 8), rc=status) + call formatter%get_var('max_lat', rTable_(:, 9), rc=status) + call formatter%get_var('elev', rTable_(:,10), rc=status) + endif + if (present(AVR)) then + ! In GEOSgcm, it already assumes ng = 2, so NumCol = 10 + NumCol = NumGlobalVars+NumLocalVars*ng + allocate(AVR(ntile, NumCol)) + AVR(:, 1) = iTable_(:,0) + ! for EASE grid, the second collum is replaced by the area + AVR(:, 2) = rTable_(:,3) + AVR(:, 3) = rTable_(:,1) + AVR(:, 4) = rTable_(:,2) + + AVR(:, 5) = iTable_(:,2) + AVR(:, 6) = iTable_(:,3) + AVR(:, 7) = rTable_(:,4) + if (ng == 1) then + AVR(:,8) = iTable_(:,4) + else + AVR(:, 8) = iTable_(:,6) + + AVR(:, 9) = iTable_(:,4) + AVR(:, 10) = iTable_(:,5) + AVR(:, 11) = rTable_(:,5) + AVR(:, 12) = iTable_(:,7) + endif + endif + + if (present(iTable)) then + call move_alloc(iTable_, iTable) + endif + + if (present(rTable)) then + call move_alloc(rTable_, rTable) + do ll = 1, ng + where ( rTable(:,3+ll) /=0.0 ) rTable(:,3+ll) = rTable(:,3)/rTable(:,3+ll) + enddo + endif + _RETURN(ESMF_SUCCESS) + end subroutine MAPL_ReadTilingNC4 + + subroutine MAPL_ReadTilingASCII(layout, FileName, GridName, NT, im, jm, n_Grids, N_PfafCat, AVR,rc) + type(ESMF_DELayout), intent(IN) :: LAYOUT + character(*), intent(IN) :: FileName + character(*), intent(out) :: GridName(:) + integer, intent(out) :: NT + integer, intent(out) :: IM(:), JM(:) + integer, intent(out) :: n_Grids + integer, intent(out) :: N_PfafCat + real, pointer, intent(out) :: AVR(:,:) ! used by GEOSgcm + integer, optional, intent(out) :: rc + + integer :: unit, status, hdr(2), N + real, pointer :: AVR_Transpose(:,:) + + UNIT = GETFILE(FILENAME, form='FORMATTED', RC=status) + _VERIFY(STATUS) + +! Total number of tiles in exchange grid +!--------------------------------------- + + call READ_PARALLEL(layout, hdr, UNIT=UNIT, rc=status) + _VERIFY(STATUS) + NT = hdr(1) + N_PfafCat = hdr(2) + + call READ_PARALLEL(layout, N_GRIDS, unit=UNIT, rc=status) + _VERIFY(STATUS) + + do N = 1, N_GRIDS + call READ_PARALLEL(layout, GridName(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + call READ_PARALLEL(layout, IM(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + call READ_PARALLEL(layout, JM(N), unit=UNIT, rc=status) + _VERIFY(STATUS) + enddo + + if(index(GridNAME(1),'EASE') /=0 ) then + allocate(AVR(NT,9), STAT=STATUS) ! 9 columns for EASE grid + _VERIFY(STATUS) + allocate(AVR_transpose(9,NT)) + else + allocate(AVR(NT,NumGlobalVars+NumLocalVars*N_GRIDS), STAT=STATUS) + _VERIFY(STATUS) + allocate(AVR_transpose(NumGlobalVars+NumLocalVars*N_GRIDS,NT), STAT=STATUS) + _VERIFY(STATUS) + endif + + call READ_PARALLEL(layout, AVR_transpose(:,:), unit=UNIT, rc=status) + AVR = transpose(AVR_transpose) + deallocate(AVR_transpose) + call FREE_FILE(UNIT) + + _RETURN(ESMF_SUCCESS) + + end subroutine MAPL_ReadTilingASCII end module NCIOMod diff --git a/generic/MAPL_Generic.F90 b/generic/MAPL_Generic.F90 index 60019fc96746..1df5c7986940 100644 --- a/generic/MAPL_Generic.F90 +++ b/generic/MAPL_Generic.F90 @@ -6061,7 +6061,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) nwrgt1 = (mpl%grid%num_readers > 1) - isNC4 = -100 + isNC4 = MAPL_FILETYPE_UNK if (on_tiles) mpl%grid%split_restart = .false. if(INDEX(FNAME,'*') == 0) then if (AmIRoot) then @@ -6111,7 +6111,7 @@ subroutine MAPL_ESMFStateReadFromFile(STATE,CLOCK,FILENAME,MPL,HDR,RC) !end if if (FileExists) then - if (isNC4 == 0) then + if (isNC4 == MAPL_FILETYPE_NC4) then filetype = 'pnc4' else if (.not.nwrgt1) then diff --git a/shared/Constants/InternalConstants.F90 b/shared/Constants/InternalConstants.F90 index bd02f1fe51d3..5b9fc763a405 100644 --- a/shared/Constants/InternalConstants.F90 +++ b/shared/Constants/InternalConstants.F90 @@ -183,6 +183,11 @@ module MAPL_InternalConstantsMod enumerator MAPL_MASK_OUT enumerator MAPL_MASK_IN endenum + + integer, parameter :: MAPL_FILETYPE_NC4 = 0 + integer, parameter :: MAPL_FILETYPE_TXT = 1 + integer, parameter :: MAPL_FILETYPE_BIN = 2 + integer, parameter :: MAPL_FILETYPE_UNK = -1 !EOP end module MAPL_InternalConstantsMod