#!/usr/bin/perl #==================================================== # # This code implements read routines for ionosphere # files from OpenGGCM meant for IPE as input, as well # as for plotting and other purposes. # # The file naming conventions are as for other OpenGGCM # output, i.e., RUN.middlefix.NNNNNN # where: RUN is the unique id of the simulation run, # middlefix for these files is 'ioa' # NNNNNN is the number of seconds since the start of the run # # File format: # for each variable: # 1. 1024 byte header # The header is all readable ascii, filled with char(0) # The header is meant to contain meta data # metadata variables are written as 'var=value', one perl line # metadata of the first record can be easily viewed with # 'head -c 1024 filename' # giving something like: # variable=prec_e_fe_1 # coordinates=SM # nphi=121 # ntheta=361 # The 4 variables shown above are mandatory and necessary # to read the files. We'll add other metadata as needed, # like UT time, creator, OpenGGCM code version, .... # 2. binary data, 4*nphi*nthe bytes # # This repeats for all variables in the file. # # The code below has an API for reading individual fields by variable name. # # Jimmy Raeder, UNH, June 2018 # # $infile="cir07_19970227_ipefile.ioa.000360"; make_f77(); system("./a.out"); sub make_f77{ $t=" c------------------------------------------ program test_ioa c------------------------------------------ real prec_e_fe_1(121*361) real ftv_rcm(121*361) !...... just show a list of variables and return call ioa_read_field('$infile', 'showfields', 1, + prec_e_fe_1, nphi, nthe, ierr) !...... read and return variable ftv_rcm call ioa_read_field('$infile', 'ftv_rcm', 0, + ftv_rcm, nphi, nthe, ierr) !...... read and return variable prec_e_fe_1 call ioa_read_field('$infile', 'prec_e_fe_1', 0, + prec_e_fe_1, nphi, nthe, ierr) stop end c--------------------------------------------------------- subroutine ioa_read_field(fin,what,iverb, + out,nphi,nthe,ierr) c--------------------------------------------------------- c c reads OpenGGCM ioa files c c fin (input) file name c what (input) desired variable c for what=showfields no data will be returned c but a list of all variables written to STDERR c iverb (input) verbosity level, 0 means quiet c out (output) real*4 array to hold the returned field c nphi (output) first dimension of out, azimuthal direction c in SM coordinates, 0 is towards the sun c ntheta (output) second dimension of out, SM co-latitude direction c ierr (output) should be checked for zero c not used right now. All error condirions cause abort. c Jimmy Raeder, UNH June 2018 c character*(*) fin,what real out(*) integer iverb byte header(1024) character*1024 c character*128 cdum,what2 c...... open file ierr=0 open(10,file=fin, status='old',access='stream',err=901) c...... loop until we have the field we want kpos=1 !.... initial file postion: byte 1 1000 continue c...... read header c=' ' k=0 read(10,pos=kpos,err=902,end=902) header kpos=kpos+1024 !.... advance file position do i=1,1024 c...... write out header if needed ! if(iverb.gt.0) then ! write(0,'(1x,i3.3,\$)') header(i) ! if(mod(i,16).eq.0) write(0,*)' ' ! endif ib=header(i) if((ib.eq.10).or.(ib.gt.31)) then k=k+1 c(k:k)=char(ib) endif enddo call ioa_parse(c,'variable',what2 ,idum,rdum) call ioa_parse(c,'nphi',cdum,nphi,rdum) call ioa_parse(c,'ntheta',cdum,nthe,rdum) if(iverb.gt.0)write(0,*)'ioa_read_field parsed: ',TRIM(what2),nphi,nthe if(nphi*nthe.le.0) stop if( TRIM(what) .eq. TRIM(what2) ) then read(10,pos=kpos,end=902,err=902) (out(i),i=1,nphi*nthe) return endif kpos=kpos+nphi*nthe*4 goto 1000 901 continue write(0,*)'error:ioa_read_field: file ',TRIM(fin),' does not exist' stop 902 continue if( TRIM(what) .eq. 'showfields' ) return !... we ran into EOF, but that's expected write(0,*)'error:ioa_read_field: file ',TRIM(fin),' read error, likely EOF' write(0,*)'error:ioa_read_field: file ',TRIM(fin),' likely field not found' stop end c--------------------------------------------------------- subroutine ioa_parse(c,what,cout,iout,rout) c--------------------------------------------------------- character*(*) c,what,cout character*1024 w1,w2 integer ista(3) ! write(0,*)'ioa_parse start ',TRIM(what) isnum=1 isig=1 k1=0 k2=0 isleft=1 w1=' ' w2=' ' do i=1,len(TRIM(c)) ib=ichar(c(i:i)) ! write(0,*)'next char: ',ib if(ib.eq.10) then !..... end of line, see what we got ! write(0,*)'parse found |',TRIM(what),'| |',TRIM(w1),'| --> |',TRIM(w2),'|' if(TRIM(w1).eq.TRIM(what)) then !..... we found it cout=w2 lr=len(TRIM(w2)) isnot=0 !... see if it is an integer iout=0 do j=1,lr ic=ichar(w2(j:j)) if(ic.le.32) cycle !..... skip blanks and ctrl chars if((j.eq.1).and.(ic.eq.45)) then !.... minus in first place isig=-1 else if((ic.lt.48).or.(ic.ge.58)) isnot=1 !..... not a digit !... pretend its a digit but avoid overflow iout=10*min(iout,99999999)+ic-48 endif ! write(0,*)'looking for int ',j,w2(j:j),' ',isnot,iout,isig enddo if(isnot.eq.0) then !...... Yup, its an integer iout=isig*iout rout=iout !...... could be meant to be a real cout=w2 ! write(0,*)'found int ',iout,rout,cout return endif !..... maybe a real read(w2,*,err=900)rout !..... should jump if not a proper real iout=rout ! not tested cout=w2 return 900 continue !..... character if none of the above iout=0 rout=0 cout=w2 return endif !... we found it k1=0 !..... didnt find it, reset everything k2=0 isig=1 isleft=1 w1=' ' w2=' ' else if(ib.eq.61) then !..... equal isleft=0 else if(isleft.eq.1)then !..... left or right of equal sign k1=k1+1 w1(k1:k1)=char(ib) else k2=k2+1 w2(k2:k2)=char(ib) endif endif enddo write(0,*)'error:ioa_parse:cant find ',TRIM(what) stop end "; system("/bin/rm -vf a.out tmp.f"); open(FF,">tmp.f"); print FF $t; close(FF); $k=1; foreach(split('\n',$t)){ printf "%6.6d %s\n",$k++,$_; } system("gfortran -w -ffixed-line-length-512 tmp.f"); }