;------------------------------------------------------------------------------ ; NCL4GC: NCL Tools for use with the GEOS-Chem Chemical Transport Model ! ;------------------------------------------------------------------------------ ;BOP ; ; !ROUTINE: read_bpch.ncl ; ; !DESCRIPTION: Reads a binary punch file (which is the tradtional GEOS-Chem ; output file format) and will (eventually) convert it to a netCDF file. ;\\ ;\\ ; !USES: ; ; NCL code load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; Local code load "$NCL4GC/file_io/get_char_from_bpch.ncl" load "$NCL4GC/file_io/get_float_from_bpch.ncl" load "$NCL4GC/file_io/get_string_from_bpch.ncl" ; External Fortran code (built into shared libraries w/ WRAPIT) external BIN "$NCL4GC/fortran/binary.so" ; ; !INTERFACE: ; undef( "read_bpch" ) function read_bpch( inFile:string, num_tracer:integer ) ; ; !INPUT PARAMETERS: ; inFile: Path of the bpch file to read from disk. ; ; !LOCAL VARIABLES: ; local numRec, numDataBlocks, R, fti, toptitle, \ N, bytes, c_modelname, lon, lat, \ halfpolar, center180, modelname, c_category, c_unit, \ c_reserved, tau0, tau1, tracer, NI, \ NJ, NL, I0, J0, L0, \ skip, category, unit, reserved, data ; ; !CALLING SEQUENCE: ; ncl 'inFile="myfile.bpch"' read_bpch.ncl ; !REMARKS: ; (1) Requires the bpch.so shared object file in the ../fortran directory. ; To build this library, issue these commands: ; (a) cd ../fortran ; (b) WRAPIT -in bpch.f # if using the IFORT compiler, or ; WRAPIT -pg bpch.f # if using the PGI compiler ; (2) As of 05 Dec 2013, this file is rough around the edges. It reads ; data properly from the bpch file, but we have not added code to save ; to a netCDF file. ; (3) This file will also eventually be a subroutine but now it is listed ; as a main program. The code to convert this to a subroutine is ; commented out above. ; ; !REVISION HISTORY: ; 05 Dec 2013 - R. Yantosca - Initial version ; 12 Mar 2014 - R. Yantosca - Now get Fortran routines from binary.so ;EOP ;------------------------------------------------------------------------------ ;BOC begin ; Make sure that the file name is passed as an argument if ( .not. isvar( "inFile" ) ) then print( "READ_BPCH: Need to pass an input file!" ) exit end if ; Tell NCL to expect big-endian data setfileoption( "bin","ReadByteOrder", "BigEndian" ) setfileoption( "bin","RecordMarkerSize", 4 ) ; Get # of binary records in the file. ; A binary record is data written in one single WRITE statement. numRec = fbinnumrec( inFile ) ; Exit if there are no data records if ( numRec .eq. 0 ) then print( "READ_BPCH: Could not find any data records in" + inFile ) exit end if ; Get the # of data blocks in the file. Each data block consists of 3 ; records (two header lines and the data array). Also, we need to skip ; over the file type identifier (FTI) and TOPTITLE field, which only ; occur at the very top of the bpch file. numDataBlocks = ( numRec - 2 ) / 3 ; Initialize record counter R = 0 ; Read top-of-file fields fti := get_string_from_bpch( inFile, R, 40 ) toptitle := get_string_from_bpch( inFile, R, 80 ) ; Loop over the number of data blocks in the file ; (there are 3 lines per data block header) do N = 1, numDataBlocks ;======================================================================= ; Parse first line of BPCH data block header ;======================================================================= ; Read the complete first line of each data block header ; Return this as an array of characters (i.e. unsigned bytes) bytes := get_char_from_bpch( inFile, R, -1 ) ; Initialize variables c_modelname := new( 20, character ) lon := 0.0 lat := 0.0 halfpolar := 0 center180 := 0 ; Call a Fortran routine to parse bytes into separate variables BIN::bpch_hdr1( bytes, c_modelname, lon, lat, halfpolar, center180 ) ; Convert character arrays to strings modelname = chartostring( c_modelname ) ; Debug print ;print( modelname ) ;print( lon ) ;print( lat ) ;print( halfpolar ) ;print( center180 ) ;======================================================================= ; Parse second line of BPCH data block header ;======================================================================= ; Read the complete second line of each data block header ; Return this as an array of characters (i.e. unsigned bytes) bytes := get_char_from_bpch( inFile, R, -1 ) ; Initialize c_category := new( 40, character ) c_unit := new( 40, character ) c_reserved := new( 40, character ) tau0 := 0d tau1 := 0d lon := 0.0 lat := 0.0 tracer := 0 NI := 0 NJ := 0 NL := 0 I0 := 0 J0 := 0 L0 := 0 skip := 0 ; Call a Fortran routine to parse bytes into separate variables BIN::bpch_hdr2( bytes, c_category, tracer, c_unit, tau0, tau1, c_reserved,\ NI, NJ, NL, I0, J0, L0, skip ) ; Convert character arrays to strings category = chartostring( c_category ) unit = chartostring( c_unit ) reserved = chartostring( c_reserved ) ; Debug print if (tracer.eq. num_tracer) then print( "category : " + category ) print( "tracer : " + tracer ) print( "unit : " + unit ) print( "tau0 : " + tau0 ) print( "tau1 : " + tau1 ) print( "reserved : " + reserved ) print( "NI NJ NL : " + NI + " : " + NJ + " : " + NL ) print( "I0 J0 L0 : " + I0 + " : " + J0 + " : " + L0 ) print( "skip : " + skip ) end if ;======================================================================= ; Read data array from BPCH file ;======================================================================= if ( NL .eq. 1 ) then ; Read 2D data. NOTE: we have to list the dimensions in ; row-major order, which is the reverse of FORTRAN order. data = get_float_from_bpch( inFile, R, (/ NJ, NI /) ) else ; Read 3D data. NOTE: we have to list the dimensions in ; row-major order, which is the reverse of FORTRAN order. data = get_float_from_bpch( inFile, R, (/ NL, NJ, NI /) ) end if ; Debug if (tracer .eq. num_tracer) then print( min(data) ) print( max(data) ) return (data) end if delete( data ) end do end ;EOC