!> hdftake-para.F !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !> mhd_open_read_hdf !! !! Opens the hdf file and reads the data !! !! if modified julian date mjd > 0, mhd_open_write_hdf !! takes the modified julian date as valid input. !! function MHD_OPEN_READ_HDF(basename, suffix_1, suffix_2, $ mjd, iStep,time_of_step) implicit none ! ... Local variables .................................................. integer NI,NJ,NK,no,li,lj,lk integer llow, lihigh, ljhigh integer num_procs, num_mhd, num_ion integer ixx real iyy integer istep,ierror,ierr,len,iistep,ini,inj,ink,i,j,k integer nsize_i,nsize_j,nsize_k,npad INTEGER mhd_open_read_hdf, hdftake_var, hdftake_close logical first, GRID character*80 fname character*80 dumpfile character*2048 namer_in character*256 run_ident, input_file logical SAMEOLD, SAMERUN #ifdef BCAST #include "mpif.h" #endif #include "hdf.inc" #include "boundx.inc" #include "global_dims.inc" #include "help.inc" #include "param.inc" #include "run-time.inc" #include "dipole.inc" #include "runattr.inc" c declaration (implicit none) real x2,y2,phi,x2ion,y2ion,xinterp,yinterp c dimension x2((ni_global+1),(nj_global+1)), $ y2((ni_global+1),(nj_global+1)),phi(nk_global+1) #ifdef ION_ON dimension x2ion(nj+1,nk/2+3),y2ion(nj+1,nk/2+3) dimension xinterp(mion_i,mion_j),yinterp(mion_i,mion_j) #endif integer*4 iMJDId integer*4, save :: iFileId integer*4 iSecId,iSecIndex, iaDimId, iAttrId, iRank integer iNtype, iNAttr integer sfstart, sfend, sfscatt, sfcreate, sfwcdata integer sfwdata, sfendacc,sffattr,sfrnatt,sfrcatt integer sfn2index,sfginfo,sfrdata,sfselect integer sfsattr, sfrattr integer iLen, iStatus integer iFileStart, iFileStop, iFileInterval integer iaDim4d(4), iaStart4d(4), iaEdge4d(4),iaStride4d(4) integer iaDim3d(3), iaStart3d(3), iaEdge3d(3),iaStride3d(3) integer iaDim2d(2), iaStart2d(2), iaEdge2d(2),iaStride2d(2) integer iaDim1d(1), iaStart1d(1), iaEdge1d(1),iaStride1d(1) character*80 caCont, caDim, caSec, caAttr, dipole character*20 caTmp logical InFile c declarations for hdftake_var entry character *(40) VARNAME real dumvar((-npad+1):(nsize_i+npad),-npad+1:nsize_j+npad, $ -npad+1:nsize_k+npad) * real dumvar(-3:55,-3:53,-3:69) real, save :: realdummy(-3:55,-3:53,-3:69) real*4, allocatable :: dumhdf(:,:,:) ! ... Parameter Variables .............................................. character*128 basename, suffix_1, suffix_2 real*8 mjd ! Additional declarations (implit none!) real*8 time_of_step,ttime real*4 time_of_step4,ttime4 ! ! ... Begin ............................................................ ! #ifdef DEBUG_MODE_ON write(*,*) "DEBUG: in hdftake-para.F::mhd_open_read_hdf(...)" #endif c c check to see if this is the I/O processor c if ( mype .eq. 0 ) then c c !open and setup hdf file c ixx = MPI_INTEGER iyy = MPI_REAL #ifdef DEBUG_MODE_ON write(6,*) 'hdftake-para: Opening HDF file for reading' write(6,*) ' MPI Variables',ixx,iyy write(6,*) MPI_INTEGER, MPI_REAL #endif ilen = index(basename,char(0))-1 if (mjd < 0) then write(input_file, 976) basename(1:ilen), suffix_1(1:7), $ suffix_2(1:3) ilen = ilen + 16 else write(input_file, 976) basename(1:ilen), suffix_1(1:20), $ suffix_2(1:3) ilen = ilen + 29 endif 976 format(A, '_mhd_', A, '.', A) write(*,977) input_file(1:ilen) 977 format("Reading, '", A, "'") iFileId = sfstart(input_file(1:ilen),DFACC_RDONLY) if (iFileId .lt. 0 ) then write(6,979) input_file(1:ilen) 979 format('Unable to open ',A,' for read!') call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif !determine who wrote the file iSecId = sffattr(iFileId,'run_descriptor') iStatus = sfrcatt(iFileId,iSecId,run_label) #ifdef DEBUG_MODE_ON write (6,*) run_label #endif !determine dipole moment iSecId = sffattr(iFileId,'dipole_moment') iStatus = sfrcatt(iFileId,iSecId,dipole) read (dipole,FMT=101) geoqmu 101 format('Geoqmu (cgs) = ',1pe16.7) #ifdef DEBUG_MODE_ON write (6,*) ' Dipole moment', geoqmu #endif iSecId = sffattr(iFileId,'written_by') iStatus = sfrcatt(iFileId,iSecId,caAttr) if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) if (caAttr(1:1) .eq. 'C') then write(6,*) ' Unable to restart from C files' call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif c c c c check to see if this is the I/O processor c c >jgl test write (9,1911) mype 1911 format ('In mhd_open_read ',i3) c c !open and setup hdf file c !determine the current run history iSecId = sffattr(iFileId,'run_history') iStatus = sfrcatt(iFileId,iSecId,run_history) * write (6,*) run_history !determine the past INPUT1 file iSecId = sffattr(iFileId,'run_namelists') iStatus = sfrcatt(iFileId,iSecId,namer_in) #ifdef DEBUG_MODE_ON write (6,*) namer_in #endif iSecId = sffattr(iFileId,'written_by') iStatus = sfrcatt(iFileId,iSecId,caAttr) if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) if (caAttr(1:1) .eq. 'C') then write(6,*) ' Unable to restart from C files' call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif if ( mjd > 0 ) then iMJDId = sffattr(iFileId,'mjd') iStatus = sfrcatt(iFileId,iMJDId,mjd) if (iStatus .eq. -1) then write(*,*) "Error reading MJD from mhd file." call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif endif iSecId = sffattr(iFileId,'time_8byte') !! iSecId will fail if old version of hdf file, check for old "time" attribute if ( iSecId .eq. -1) then iSecId = sffattr(iFileId,'time') iStatus = sfrcatt(iFileId,iSecId,time_of_step4) time_of_step = time_of_step4 else iStatus = sfrcatt(iFileId,iSecId,time_of_step) endif #ifdef DEBUG_MODE_ON write (6,*) ' time',istatus,iFileId,isecid,caAttr #endif if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) iAttrId = sffattr(iFileId,'time_step') iStatus = sfrnatt(iFileId,iAttrId,iStep) lstep = iStep !! check for old style time iAttrId = sffattr(iFileId,'time_8byte') if (iAttrId .eq. -1 ) then iAttrId = sffattr(iFileId,'time') iStatus = sfrnatt(iFileId,iAttrId,time_of_step4) time_of_step = time_of_step4 else iStatus = sfrnatt(iFileId,iAttrId,time_of_step) endif time = time_of_step c tzero (defined in boundx.inc) iAttrId = sffattr(iFileId, 'tzero') iStatus = sfrnatt(iFileId,iAttrId,tzero) if ((iStatus .eq. -1) .or. (iAttrId .eq. -1)) then write(*,*) "Error reading tzero from mhd file." call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif c endif ! ends work on I/O processor ( = 0 ) #ifdef BCAST c broadcast the dipole moment to everyone call MPI_BCAST(geoqmu,1,MPI_REAL,0, $ MPI_COMM_WORLD,ierror) #ifdef DEBUG_MODE_ON if ( mype .eq. 0 ) $ write (6,*) ' MPI broadcast geoqmu ',ierror #endif c c broadcast the time and time step to everyone call MPI_BCAST(lstep,1,MPI_INTEGER,0, $ MPI_COMM_WORLD,ierror) call MPI_BCAST(time,1,MPI_DOUBLE_PRECISION,0, $ MPI_COMM_WORLD,ierror) call MPI_BCAST(iStep,1,MPI_INTEGER,0, $ MPI_COMM_WORLD,ierror) call MPI_BCAST(time_of_step,1,MPI_DOUBLE_PRECISION,0, $ MPI_COMM_WORLD,ierror) c tzero (defined in boundx.inc) call MPI_BCAST(tzero,1,MPI_REAL,0, $ MPI_COMM_WORLD,ierror) #ifdef DEBUG_MODE_ON if (mype .eq. 0) $ write (6,*) ' MPI broadcast step and time ',ierror #endif c broadcast the modified julian date (mjd) to everyone call MPI_BCAST(mjd,1,MPI_DOUBLE_PRECISION,0, $ MPI_COMM_WORLD,ierr) #ifdef DEBUG_MODE_ON if (mype .eq. 0) $ write(6,*) ' MPI broadcast mjd ',ierror #endif #endif c mhd_open_read_hdf = iStatus ! ... End FUNCTION mhd_open_read_hdf ................................... ! return * * ************ read in variables ! ! ... Entry statement .................................................. ! Entry HDFTAKE_VAR(varname, dumvar,nsize_i,nsize_j,nsize_k, $ npad,iiStep,ttime,GRID) c c check for I/O processor c if ( mype .ne. 0 ) return c c >jgl test write (9,1912) mype,varname 1912 format ('In hdfdump_var ',i3,a40) c ! Read in variable for timestep len = index(VARNAME,char(0))-1 call rststr(caSec) caSec(1:len) = VARNAME(1:len) if ( GRID ) then caSec(len+1:len+4) = 'grid' len = len+4 endif * write(6,*) istart,icatmp iSecIndex = sfn2index(iFileId,caSec(1:len)) iSecId = sfselect(iFileId, iSecIndex) iAttrId = sffattr(iFileId,'time_step') iStatus = sfrnatt(iFileId,iAttrId,iiStep) #ifdef DEBUG_MODE_ON write (6,*) ' time_step',istatus,iFileId,isecid,caAttr #endif if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) iAttrId = sffattr(iFileId,'time_8byte') if ( iAttrId .eq. -1 ) then iAttrId = sffattr(iFileId,'time') iStatus = sfrnatt(iFileId,iAttrId,ttime4) ttime = ttime4 else iStatus = sfrnatt(iFileId,iAttrId,ttime) endif if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) iStatus = sfginfo(iSecId,caSec,iRank,iaDim3d,iNtype,iNAttr) if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) iAttrId = sffattr(iSecId,'ni') iStatus = sfrnatt(iSecId,iAttrId,iNi) iAttrId = sffattr(iSecId,'nj') iStatus = sfrnatt(iSecId,iAttrId,iNj) iAttrId = sffattr(iSecId,'nk') iStatus = sfrnatt(iSecId,iAttrId,iNk) allocate( dumhdf(iaDim3d(1),iaDim3d(2),iaDim3d(3))) ! ! Note: LFM-RCM has extended the MHDManager class to use HDFDUMP_VAR ! and HDFTAKE_VAR to read/write the accumulate & bleed arrays ! necessary for doing an LFM-RCM restart. The bleed arrays are of ! length {ni,nj,nk}, but this code will print errors to stdout ! because those arrays aren't of length {ni+1,nj+1,nk+1}. ! ! I've extended this check to pass for variables of size {ni,nj,nk} ! -- schmitt ucar edu; Oct 12, 2009 ! if (((iaDim3d(1) .ne. nsize_i) .or. $ (iaDim3d(2) .ne. nsize_j) .or. $ (iaDim3d(3) .ne. nsize_k) ) .and. $ ((iaDim3d(1) .ne. nsize_i-1) .or. $ (iaDim3d(2) .ne. nsize_j-1) .or. $ (iaDim3d(3) .ne. nsize_k-1) ) ) then write (6,*) ' dimensions in hdf file might be ', $ ' inconsistent with code dims for variable ', $ varname write (6,*) ' index code file' write(6,*) ' 1 ',nsize_i,iaDim3d(1) write(6,*) ' 2 ',nsize_j,iaDim3d(2) write(6,*) ' 3 ',nsize_k,iaDim3d(3) call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) endif do i=1,3 iaStride3d(i) = 1 iaStart3d(i) = 0 iaEdge3d(i) = iaDim3d(i) enddo iStatus = sfrdata(iSecId,iaStart3d,iaStride3d,iaEdge3d,dumhdf) if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) iStatus = sfendacc(iSecId) if (iStatus .eq. -1) $ call MPI_ABORT(MPI_COMM_WORLD,MPI_ERR_OTHER) do k=1,iaDim3d(3) do j = 1, iaDim3d(2) do i = 1,iaDim3d(1) dumvar(i,j,k) = dumhdf(i,j,k) enddo enddo enddo c deallocate(dumhdf) #ifdef DEBUG_MODE_ON write(6,917) varname(1:10),iiStep,ttime 917 format('HDFTAKE: ',A10,' ',i7,' ',f10.1) * write(6,*) ' Done with HDFTAKE' #endif HDFTAKE_VAR = iStatus ! ! ... End FUNCTION mhd_open_read_hdf ................................... ! return ! ! ... Entry Statement .................................................. entry hdftake_close !close file on last dump c c check to see if this is the I/O processor c if ( mype .ne. 0 ) return c iStatus = sfend(iFileId) hdftake_close = iStatus ! ... End FUNCTION mhd_open_read_hdf ................................... return end !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!