C Copyright(c) 2004, Space Science and Engineering Center, UW-Madison C Refer to "McIDAS Software Acquisition and Distribution Policies" C in the file mcidas/data/license.txt C *** $Id: imgmake.f,v 1.2 2012/03/20 16:46:23 russd Tst $ *** C ? IMGMAKE -- Creates a McIDAS area from a binary file C ? IMGMAKE sfile ddataset lines elems C ? Parameters C ? sfile | source image file name C ? ddataset | destination ADDE dataset name and position; specify as C ? alias.position or group/descriptor.position; only positive C ? integers (absolute positions) or the default are valid for C ? the position number; to use the default position, omit the C ? .position portion (no def for alias or group/descriptor, for C ? position def=the position following the most recent image) C ? * lines | number of lines in the source image file C ? * elems | number of elements per line in the source image file C ? Keywords C ? ********************* Source File Keywords *********************** C ? * FLIp=line elem | flip source lines and/or elements (def=NO NO) C ? * FORmat= | source data format; INTEGER or FLOAT; (def=INTEGER) C ? * GAIN= | multiplicative factor to apply to each source element (def=1.0) C ? * HEAder= | byte length of source file header (def=0) C ? * HFIle= | name of source header file (def=no header file) C ? * INTerleave= | data interleave format; (def=NON) C ? 'NON' = no interleave C ? 'BIL' = Band Interleave Line C ? * MAG=lmag emag | line and element magnification factores (def=1 1) C ? * MASk= | sfile mask used to extract day/time info (def=no mask) C ? * MISs=in1 out1 | missing input and output data values C ? (def=all data values are valid ) C ? * NBPE= | Number of Bytes Per Element in source image (def=1) C ? * NFILe=file format | source LALO navigation file name and format C ? BINary or ASCii; see Remarks (def=X BIN) C ? * OFFset= | additive constant to apply to each source element (def=0.0) C POST=OFF | Turn OFF post process (def=ON) C ? PROfile= | name of a context file parameter/keyword default profile C ? (def=use internal keyword defaults) C ? * SCAle=slo shi dlo dhi | linear scaling of the source data range to C ? a destination data range (def=0 255 0 255) C ? * SBAND= | source band number (def=destination band number) C ? * SWByte=data nav | switch bytes of DATA and NAVigation file contents C ? from big endian to little endian format (def=NO YES) C ? * TABle=file | file containing mapping of source BRIT to TEMP (K) C ? used for conversion to SSEC standard BRIT (def=no file) C ? * TYPE=type sub_type | source image file type and sub_type; see Remarks C ? (def=RAS ) C ? ***************** Destination Dataset Keywords ******************* C ? * ADD= | additive constant to apply to each destination element (def=0.0) C ? * AUX= offset | byte offset to start of Auxillary block (def=0) C ? * CAL= | byte offset to the CAL block (def=768) C ? * CTYpe= | calibration type (def=BRIT) C ? * BANd= | destination band number (def=1) C ? * DATa= | byte offset to the DATA block (def=1280) C ? * DAY= | image day (def=current day) C ? * ELEm= | beginning image element number (def=1) C ? * EREs= | element resolution (def=line resolution) C ? * LINe= | beginning image line number (def=1) C ? * LREs= | line resolution (def=1) C ? * MEMo='text' | memo field (def=source file name) C ? * MUL= | multiplicative factor to apply to each destination element (def=1.0) C ? * NBYte= | number of bytes per data element; 1, 2 or 4 (def=1) C ? * MAKCAL=type prodlo prodhi britlo brithi unit scale miss | make a PRD C ? cal entry for the destination PRD image; see Remarks (def=no cal) C ? * MAKNAV=proj | make a navigation block based on a geographic C ? projection, see Remarks (def=no projection) C ? * NAV= | byte offset to the NAV block (def=256) C ? * PAD=nline value | pad the top of the destination image with lines C ? containing a constant value (def=0 0) C ? * PLAnet= | planet for geographic naviagtion (def=EARTH) C ? * SS= | sensor source number (def=0) C ? * STYpe= | source type (def=PRD) C ? * TIMe= | image time (def=current hhmm00) C ? C ? Remarks C ? The input image must be decompressed. C ? C ? * Denotes parameters and keywords that can be set by a profile within C ? a IMGMAKE.CORE, IMGMAKE.SITE or IMGMAKE.USER context file. C ? C ? Valid image source TYPEs and Sub_Types are: C ? RASter - 2D array of image lines and elements C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? KALPANA KALPANA-1 images C ? C ? NSSl - 3D array with self describing header; values for LINES, C ? ELEMS, HEADER, DAY and TIME are extracted from header. C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? C ? ENVi - 2D or 3D array with self describing header; values for LINES, C ? ELEMS, HEADER, and INTERLEAVE are extracted from header. C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? AIRS AIRS Relative Humidity images C ? DRS DRS Radar images C ? C ? TIF - TIF format image file; values for LINES, ELEMS and C ? HEADER are extracted from the header C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? C ? RUN - 2D RUN length encoded source file; a temporary RASTER type C ? image is created to hold the decoded elements; Values for C ? LINES and ELEMS are defined internally. C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? MOSIAC UNISYS NEXRAD national composite C ? C ? VMC - 2D array of image lines and elements; elements are C ? 2bytes; nlines, nelems, header are read from file C ? header C ? Sub_Type Description C ? -------- ---------------------------------------------- C ? C ? Predefined profiles are contained in the context files IMGMAKE.CORE, C ? IMGMAKE.SITE and IMGMAKE.USER. The .USER file overrides the .SITE C ? file, which overrides the .CORE file. The .CORE file is supplied with C ? McIDAS-X and should not be modified. You can, however, change the C ? settings of its products and/or create your own products. To do so, C ? copy the .CORE file from ~mcidas/data to your $HOME/mcidas/data or C ? another MCPATH directory, rename it to IMGMAKE.USER, then edit the C ? existing products and/or add new ones. C ? C ? The MASK= keyword defines a text string that marks the location of C ? day and time values embedded within the source file name. If specified, C ? this string will be used to extract the DAY and TIME of the image C ? from the source file name. The format for a MASK entry is as follows: C ? MASK=XXXX(tag)XXXX(tag)... C ? XXXX - used to pad the position of day/time values C ? within the source file name string. C ? (tag) - predefined set of day/time match strings: C ? (CC) = 2 digit century C ? (YY) = 2 digit year C ? (DD) = 2 digit day of month C ? (MN) = 2 digit month of year C ? (HH) = 2 digit hour of day C ? (MM) = 2 digit minute of hour C ? (SS) = 2 digit second of minute C ? (DDD) = 3 digit Julian day C ? *char = variable length field truncated by char C ? Example: for the source files named: C ? MET7_VIS_2004_03_30_1200 C ? MET7_IR_2004_03_30_1200 C ? use the MASK keyword entry: C ? MASK=XXXXX*_(CC)(YY)X(MN)X(DD)X(HH)(MM) C ? C ? SWBYTE= is a no-op on a big endian platform. On a little endian C ? platform this will reverse the byte order of all 4-byte words within the C ? data block and/or navigation block. C ? C ? MAG= magnifies the number of source file lines and elements before writing C ? the data to the ADDE image object. Values other than 1 will modify the C ? values for number of image lines, number of pixels per line, line C ? resolution and element resolution to provide correct Earth coordinate C ? transforms. Values specified for line/element size and resolution are C ? relative to the contents of the source file prior to a magnification. C ? C ? MAKCAL= creates a product calibration for the destination image based C ? on a linear mapping of the image data values to matching linear output C ? range. Parameters are: C ? type | name of product defined by the segment, 1-4 characters, C ? e.g., P, RH, TEMP (def=CVAL) C ? prodlo | minimum output product value (def=0) C ? prodhi | maximum output product value (def=255) C ? britlo | input brightness value assigned to 'prodlo' (def=0) C ? brithi | input brightness value assigned to 'prodhi' (def=255) C ? unit | name of units for product value defined by the segment; 1-4 C ? characters, e.g., F, MB, % C ? scale | divisor to apply to the product values when they are accessed C ? by another application, e.g., D or IMGPROBE; thus C ? output_value=product_value/scale (def=0, meaning no scaling) C ? miss | missing value to assign to product when a brightness value is C ? outside the range defined with 'britlo' and 'brithi'; use C ? decimal value -2139062144 for HEX80 (def=0) C ? C ? MAKNAV= creates navigation entry for the destination image based on a C ? geographic map projection. Available navigation types with associated C ? parameters are: C ? proj parms C ? ---- ------------------------------------------------------------- C ? LAMB npl npe slat1 slat2 qlon size pole C ? MERC eql eqe slat qlon size C ? PS npl npe slat qlon size pole C ? RADR cenl cene clat clon size rot C ? RECT lin lat ele lon deglin degele C ? SIN cenl cene clat clon size C ? TANC npl npe slat qlon size C ? MOLL eql eqe qlon size C ? NOAA kepler day time tmsec orient flip C ? where: C ? cenl,cene - line (cenl) and element (cenl) of radar site C ? clat,clon - lat (clat) and lon (clon) of radar site C ? deglin,degele - degrees/line (deglin) and element (degele) C ? ele,lon - image element (ele) of longitude (lon) C ? eql,eqe - image line (eql) and element (eqe) of equator C ? lin,lat - image line (lin) of latitude (lat) C ? npl,npe - image line (npl) and element (npe) of pole C ? pole - N(orth) or S(outh) pole (def=N) C ? qlon - normal longitude (center) C ? rot - rotation of radar picture (def=0) C ? size - size of one pixel (km) at standard latitude C ? - or map scale factor (ratio) C ? slat - standard latitude (PS def=60, MERC def=0) C ? slat1- standard latitude (LAMB def=30) C ? slat2- standard latitude (LAMB def=60) C ? kepler - name of 3-line Keplerian element file C ? day - nominal DAY for image C ? time - nominal TIME for image C ? tmsec - time of first line to navigate [msec from DAY] C ? orient - ascending(+1)/descending(-1) pass (def=+1) C ? flip - inverted image flag (0-normal/1-flipped) (def=0) C ? C ? Use the NFILE keyword to create a LALO navigation from an input C ? file which contains a latitude and longitude (west negative) for C ? for each image element. C ? ASCII type nav file: 1 lat/lon per text line C ? NFILE= file ASC lap lal lop lol skip res C ? lap - latitude starting char position on a text line C ? lal - latitude max char length C ? lop - longitude starting char position on a text line C ? lol - longitude max char length C ? skip - number of text header lines to skip C ? res - resolution reduction factor (def=1) C ? BINARY type nav file: 2x2D (lat/lon) array of 4byte floats C ? NFILE= file BIN C ? file - name of file containing Lat/Lon pairs C ? BIN - identifies files as BINARY format C ? NFILE= Latfile Lonfile BIN lines elems imgLine imgElem C ? Latfile - name of file containing Latitude values C ? Lonfile - name of file containing Longitude values C ? BIN - identifies files as BINARY format C ? lines - line dimension of nav arrays C ? elems - element diension of nav arrays C ? imgLine - image line coordinate of first nav element C ? imgElem - image element coordinate of first nav element C ? C ? IMGMAKE will automatically PAD fill the bottom of any "short" C ? image with the pad_value. A short image is defined as any C ? source file that does not contain the stated number of source C ? bytes (lines x elems x NBPE). C ? ---------- subroutine main0 implicit none INCLUDE 'areaparm.inc' c --- constants integer MAXNAVSIZE parameter (MAXNAVSIZE = 40000) integer MAXCALSIZE parameter (MAXCALSIZE = 10000) integer MAXDATSIZE parameter (MAXDATSIZE = 80000) integer MAXLALOSIZE parameter (MAXLALOSIZE = 5000000) integer MAXCALAUXSIZE parameter (MAXCALAUXSIZE = 5000000) c --- external functions character*4 clit character*12 cff character*12 cfz character*12 cfe character*12 cfu character*12 cfi character*12 cfj integer lbi, lwi integer len_trim integer lit integer lwfile integer mcgetdaytime integer mcstrtodbl integer mccmd integer mccmdnum integer mccmdiyd integer mccmdint integer mccmddbl integer mccmdihr integer mccmdirng integer mccydtoiyd integer mccydtodmy integer mccmdstr integer mcskeyin integer mcstrtoint integer m0dsnnam integer m0arasize integer m0setdefaultfiles integer m0getdefaultkeynumparms integer m0getdefaultargstr integer m0getdefaultargint integer m0getdefaultargdbl integer m0getdefaultargiyd integer m0getdefaultargihr integer mcaput integer mcaout integer mcacou c --- internal functions integer gryscl integer list_TIF integer list_NSSL integer list_ENVI integer list_VMC integer dec_RUN integer filedaytime integer filemask integer maknav integer read_lookup_cal integer read_ascii_lalo integer vmcread integer readENVInav c --- internal variables character*12 imgorbit character*12 imgnumber character*12 imgchannel character*12 swbyte_data character*12 cal_num character*12 cal_type character*12 cal_prodlo character*12 cal_prodhi character*12 cal_britlo character*12 cal_brithi character*12 cal_unit character*12 cal_scale character*12 cal_miss character*12 cvalue character*12 cvalue2 character*12 cval1, cval2, cval3 character*12 ctype, sub_type character*12 stype character*12 cday character*12 clat, clon character*12 ctime character*12 proj character*12 cmiss, cmoss character*12 cmiss2, cmoss2 character*12 format character*12 interleave character*12 context_envvar character*12 context_file(3) character*12 ntype character*32 cmemo character*80 profile character*80 sorts(20) character*80 mask character*80 link character*80 linkstring character*80 linkprofile character*80 fields(20) character*80 card character*256 cmd character*256 sfile character*256 tfile character*256 sheader character*256 nfile, nfile2 character*256 lookup_file character*256 ddataset character*256 ddataset_nam character*256 ddataset_pos character*256 ddataset_ind double precision dvalue double precision gain double precision offset double precision mul double precision add double precision input_min, input_max double precision output_min, output_max double precision input_low, input_high double precision output_low, output_high double precision scale double precision input_miss, output_miss double precision input_miss2, output_miss2 double precision cal_dscale double precision cal_doffset integer nmiss integer nfields integer nblank integer len integer band integer sband integer n_file_bands, file_bands(9999) integer ndbands, dbands(64) integer nsbands, sbands(64) integer nlines_sent integer nsort integer i, ii integer j, jj integer k integer isign integer ivalue, jvalue integer ptr, iptr, optr integer day integer time, mstime integer line integer lines, lalo_lines integer pad_lines, pad_value integer lat_offset, lat_pos, lat_len integer lon_offset, lon_pos, lon_len integer skip, nav_res, lalo_res integer nlalo, lalo_size integer lbeg, lend, linc, lmag, lbu, lbd integer ebeg, eend, einc, emag, ebu, ebd integer elems, lalo_elems integer status integer header integer nbpe integer nbyte integer num integer num_links integer navsize integer in_low, in_high, in_inc integer ibuf(MAXDATSIZE) integer jbuf(MAXDATSIZE) integer nbuf(MAXNAVSIZE) integer calaux(MAXCALAUXSIZE) integer ddir(64) integer icard(20), ncard integer num_test, max_test integer mislat, mislon integer mx_value integer*2 ibuf2(MAXDATSIZE*2) real rlat, rlats(MAXLALOSIZE) real rlon, rlons(MAXLALOSIZE) real minlat, minlon real maxlat, maxlon real rvalue real rtemp real table(65535) real table2(65535) logical CONTEXT logical LINE_FLIP logical ELEM_FLIP logical TMP_FILE logical POST logical TEST logical MCIDAS_NAV c --- seperate Lat/Lon files integer lalox2_line integer lalox2_elem integer lalox2_lines integer lalox2_elems integer lalox2_size integer lalox2_lres integer lalox2_eres integer lalox2_lat integer lalox2_lon integer lalox2_latptr integer lalox2_lonptr logical LALOX2 c --- post process command integer ipos character*160 ppcmd common /PPCOM/ & ipos, & ppcmd c --- navigation common character*4 NAV_proj integer NAV_lin integer NAV_ele double precision NAV_lat double precision NAV_lon double precision NAV_deglin double precision NAV_degele common /NAVCOM/ & NAV_proj, & NAV_lin, & NAV_ele, & NAV_lat, & NAV_lon, & NAV_deglin, & NAV_degele equivalence (card, icard) equivalence (ibuf, ibuf2) c --- initialize application status to FAILED call mccodeset( 1 ) call sdest('IMGMAKE: Version 1.13',0) c --- initialize the MicroSecond Time mstime = 0 c --- initialize the navigation common call initnavcom c --- initialize the postprocess common ppcmd = ' ' c --- check for test TEST = .FALSE. status = mccmdstr('TEST',1,'N',cvalue) if( cvalue(1:1).eq.'Y' ) then TEST=.TRUE. num_test = 0 status = mccmdint('TEST',2,'MAX',1,0,-1,max_test) endif c --- get current day and time status = mcgetdaytime( day, time) time = (time/100) * 100 c --- source file TMP_FILE = .FALSE. status = mccmdstr(' ',1,' ',sfile) if( status.lt.0 ) then call edest('Source File name must be specified',0) return endif c --- check if file exists status = lwfile(sfile) if( status.eq.0 ) then call edest('Source file does not exist: '//sfile,0) return else call sdest('Source File: '//sfile,0) endif c --- initialize context file context_envvar = 'IMGMAKE' context_file(1) = 'IMGMAKE.USER' context_file(2) = 'IMGMAKE.SITE' context_file(3) = 'IMGMAKE.CORE' CONTEXT = .FALSE. status=m0setdefaultfiles (context_envvar, context_file, 3) if( status.eq.0 ) CONTEXT = .TRUE. c --- test for filename -> profile links profile = ' ' if( CONTEXT ) then status = m0getdefaultkeynumparms( & 'LINK', & ' ', & num_links & ) if( num_links.gt.0 ) then do i = 1, num_links status=m0getdefaultargstr( 'LINK',' ',i,link) len = len_trim( link ) ptr = index( link, ' -> ') linkstring = link(1:ptr-1) linkprofile = link(ptr+4:len) len = len_trim( linkstring) if( index(sfile,linkstring(1:len)).ne.0 ) then profile = linkprofile goto 1 endif enddo endif endif 1 continue c --- check for a named profile status = mccmdstr('PRO.FILE', 1, profile, profile) if( status.lt.0 ) return if( profile.eq.' ' ) then CONTEXT = .FALSE. else call sdest('PROFILE: '//profile,0) endif c --- destination dataset group/descritor.position status = mccmdstr(' ',2,' ',ddataset) if( status.le.0 ) then call edest('Invalid destination dataset name='//ddataset,0) return else status=m0dsnnam(ddataset,ddataset_nam,ddataset_pos, & ddataset_ind,0) if( status.lt.0 ) then call edest('Invalid dataset name ='//ddataset,0) return endif endif c --- make integer variable of dataset position status = mcstrtoint( ddataset_pos, ipos ) if( status.lt.0 ) then call sdest('FAILED -- DDataset position',0) return endif c --- get the image TYPE ctype = 'RAS' if( CONTEXT ) then status = m0getdefaultargstr( profile,'TYPE',1,ctype) endif status = mccmdstr('TYP.E',1,ctype,ctype) c --- get the image SUB-TYPE sub_type = ' ' if( CONTEXT ) then status = m0getdefaultargstr( profile,'TYPE',2,sub_type) endif status = mccmdstr('TYP.E',2,sub_type,sub_type) c --- name of source header file sheader = sfile if( CONTEXT ) then status = m0getdefaultargstr( profile,'HFILE',1,sheader) endif status = mccmdstr('HFI.LE',1,sheader,sheader) if( status.lt.0 ) return c --- initialize image parameters to most common settings cmiss = ' ' cmoss = ' ' cmiss2 = ' ' cmoss2 = ' ' cal_dscale = 1.0D0 cal_doffset = 0.0D0 cal_type = ' ' interleave = 'NON' format = 'INTEGER' nbpe = 1 n_file_bands = 1 file_bands(1) = 1 c --- RASTER if( ctype(1:3).EQ.'RAS' ) then lines = 0 elems = 0 header = 0 c --- TIFF elseif( ctype(1:3).EQ.'TIF' ) then status = list_TIF( & sfile, & lines, & elems, & header & ) if( status.lt.0 ) return c --- NSSL elseif( ctype(1:3).eq.'NSS' ) then status = list_NSSL( & sfile, & lines, & elems, & n_file_bands, & header, & day, & time & ) if( status.lt.0 ) return nbpe = 2 if( n_file_bands.gt.1 ) then call edest('Multi-level NSSL image file not supported',0) return endif c --- RUN elseif( ctype(1:3).eq.'RUN' ) then status = dec_RUN( & sfile, & tfile, & lines, & elems & ) if( status.lt.0 ) return sfile = tfile nbpe = 1 TMP_FILE = .TRUE. c --- ENVI elseif( ctype(1:3).eq.'ENV' ) then c ------ check header file name if( sheader(1:5).eq.'*.hdr' ) then len = len_trim(sfile) sheader = sfile(1:len-3)//'hdr' endif status = list_ENVI( & sheader, & lines, & elems, & n_file_bands, & header, & nbpe, & format, & interleave, & cmiss, & cmiss2, & cal_dscale, & cal_doffset & ) if( status.lt.0 ) return c ------ fill the array of file bands do i = 1, n_file_bands file_bands(i) = i enddo c ------ settings supplied from header call sdest('ENVI Header Info:',0) call sdest(' file='//sheader,0) call sdest(' lines=',lines) call sdest(' elems=',elems) call sdest(' source bands=',n_file_bands) call sdest(' header offset=',header) call sdest(' Bytes/element=',nbpe) call sdest(' format='//format,0) call sdest(' interleave='//interleave,0) call sdest(' missing 1 ='//cmiss,0) call sdest(' missing 2 ='//cmiss2,0) cmoss = cmiss cmoss2 = cmiss cvalue = cff( cal_dscale, 4 ) call sdest(' CAL scale ='//cvalue,0) cvalue = cff( cal_doffset, 4 ) call sdest(' CAL offset ='//cvalue,0) c --- VMC elseif( ctype(1:3).eq.'VMC' ) then call sdest('Decode VMC Header block',0) c ------ read the VMC header status = list_VMC( & sfile, & lines, & elems, & header, & day, & time, & imgorbit, & imgnumber, & imgchannel, & mstime & ) if( status.lt.0 ) return c --- SATURN ENVI elseif( & ctype(1:3).eq.'SAT'.or. & ctype(1:3).eq.'VIR' & ) then continue else call edest('Unsupported image type: '//ctype,0) return endif c --- source image lines if( CONTEXT ) then status = m0getdefaultargint( profile,'LINES',1,lines) endif status = mccmdint(' ',3,'Source Lines',lines,0,-1,lines) if( status.lt.0 ) return c --- source image elements if( CONTEXT ) then status = m0getdefaultargint( profile,'ELEMS',1,elems) endif status = mccmdint(' ',4,'Source Elements',elems,0,-1,elems) if( status.lt.0 ) return c --- get lalo nav lines, elems, resolution and blowdown status = mccmdint('LALO',1,'LaLo Lines',lines,0,-1,lalo_lines) status = mccmdint('LALO',2,'LaLo Elems',elems,0,-1,lalo_elems) c --- byte length of source file header if( CONTEXT ) then status = m0getdefaultargint( profile,'HEADER',1,header) endif status = mccmdint('HEA.DER',1,'Header',header,0,-1,header) if( status.lt.0 ) return c --- data format if( CONTEXT ) then status = m0getdefaultargstr( profile,'FORMAT',1,format) endif status = mccmdstr('FOR.MAT',1,format,format) if( status.lt.0 ) return c --- data interleave if( CONTEXT ) then status = m0getdefaultargstr( profile,'INTERLEAVE',1,interleave) endif status = mccmdstr('INT.ERLEAVE',1,interleave,interleave) if( status.lt.0 ) return c --- check for invalid dimensions if( lines.le.0 ) then call edest('Number of image lines must be greater than 0',0) return endif if( elems.le.0 ) then call edest('Number of image elements must be greater than 0',0) return endif if( header.lt.0 ) then call edest('Size of image header cannot be less than 0',0) return endif c --- get the file day and time from the file name mask = ' ' if( CONTEXT ) then status = m0getdefaultargstr( profile,'MASK',1,mask) endif status = mccmdstr('MAS.K',1,mask,mask) if( mask.ne.' ' ) then status = filemask ( & sfile, & mask, & day, & time, & nfields, & fields & ) if( status.lt.0 ) then call edest('Failed to decode source file mask',0) return endif call sdest('Name MASK:',0) call sdest(' day =',day) call sdest(' time=',time) if( nfields.gt.0 ) then do i = 1, nfields call sdest(' fields='//fields(i),0) enddo endif endif c --- MISS if( CONTEXT ) then status = m0getdefaultargstr( profile,'MISS',1,cvalue) if( status.ge.0 ) then if( cvalue.ne.' ' ) then cmiss = cvalue cmoss = cmiss cmoss2 = cmiss endif endif status = m0getdefaultargstr( profile,'MISS',2,cvalue) if( status.ge.0 ) then if( cvalue.ne.' ' ) then cmoss = cvalue cmoss2 = cmoss endif endif status = m0getdefaultargstr( profile,'MISS',3,cvalue) if( status.ge.0 ) then if( cvalue.ne.' ' ) then cmiss2 = cvalue endif endif status = m0getdefaultargstr( profile,'MISS',4,cvalue) if( status.ge.0 ) then if( cvalue.ne.' ' ) then cmoss2 = cvalue endif endif if( cmiss2.eq.' ' ) cmoss2 = ' ' endif status = mccmdstr('MIS.S',1,cmiss,cmiss) status = mccmdstr('MIS.S',2,cmoss,cmoss) status = mccmdstr('MIS.S',3,cmiss2,cmiss2) status = mccmdstr('MIS.S',4,cmoss2,cmoss2) if( cmiss.ne.' ' ) then status = mcstrtodbl( cmiss, input_miss ) if( status.lt.0 ) then call edest('FAILED -- Invalid Input MISSING DATA value',0) return endif endif if( cmiss2.ne.' ' ) then status = mcstrtodbl( cmiss2, input_miss2 ) if( status.lt.0 ) then call edest('FAILED -- Invalid Input MISSING DATA value2',0) return endif endif if( cmoss.ne.' ' ) then status = mcstrtodbl( cmoss, output_miss ) if( status.lt.0 ) then call edest('FAILED -- Invalid Output MISSING DATA value',0) return endif else output_miss = input_miss endif if( cmoss2.ne.' ' ) then status = mcstrtodbl( cmoss2, output_miss2 ) if( status.lt.0 ) then call edest('FAILED -- Invalid Output MISSING DATA value2',0) return endif else output_miss2 = output_miss endif call sdest('MISSing data:',0) call sdest(' value1 ='//cmiss//' -> '//cmoss,0) call sdest(' value2 ='//cmiss2//' -> '//cmoss2,0) c --- NBPE keyword if( CONTEXT ) then status = m0getdefaultargint( profile,'NBPE',1,nbpe) endif status = mccmdint('NBP.E',1,'NBPE',nbpe,1,4,nbpe) if( status.lt.0 ) return c --- MAG keyword lmag = 1 emag = lmag if( CONTEXT ) then status = m0getdefaultargint( profile,'MAG',1,lmag) status = m0getdefaultargint( profile,'MAG',2,emag) endif if( mccmdnum( 'MAG' ).ne.0 ) then status = mccmdint('MAG',1,'Line Mag',lmag,0,-1,lmag) status = mccmdint('MAG',2,'Elem Mag',lmag,0,-1,emag) endif if( lmag.eq.0 ) lmag = 1 if( lmag.lt.0 ) then lbu = 1 lbd = abs( lmag ) else lbu = lmag lbd = 1 endif if( emag.eq.0 ) emag = 1 if( emag.lt.0 ) then ebu = 1 ebd = abs( emag ) else ebu = emag ebd = 1 endif c --- FLIP keyword LINE_FLIP = .FALSE. cvalue = 'NO' if( CONTEXT ) then status = m0getdefaultargstr( profile,'FLIP',1,cvalue) endif status = mccmdstr('FLI.P',1,cvalue,cvalue) if( cvalue(1:1).eq.'Y' ) LINE_FLIP = .TRUE. ELEM_FLIP = .FALSE. cvalue = 'NO' if( CONTEXT ) then status = m0getdefaultargstr( profile,'FLIP',2,cvalue) endif status = mccmdstr('FLI.P',2,cvalue,cvalue) if( cvalue(1:1).eq.'Y' ) ELEM_FLIP = .TRUE. c --- PAD keyword pad_lines = 0 pad_value = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'PAD',1,pad_lines) status = m0getdefaultargint( profile,'PAD',2,pad_value) endif status = mccmdint('PAD',1,'Pad Lines',pad_lines,0,-1,pad_lines) if( status.lt.0 ) return status = mccmdint('PAD',2,'Pad Value',pad_value,0,-1,pad_value) if( status.lt.0 ) return c --- POST process (ON/OFF) POST = .FALSE. c cvalue = 'ON' c status = mccmdstr('POST',1,cvalue,cvalue) c if( status.lt.0 ) return c if( cvalue(1:2).eq.'OF' ) POST = .FALSE. c --- GAIN keyword gain = 1.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'GAIN',1,gain) endif status = mccmddbl('GAI.N',1,'Gain',gain,0.0D0,-1.0D0,gain) if( status.lt.0 ) return c --- OFFSET keyword offset = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'OFFSET',1,offset) endif status = mccmddbl('OFF.SET',1,'Offset',offset,0.0D0,-1.0D0,offset) if( status.lt.0 ) return c --- MUL keyword mul = 1.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MUL',1,mul) endif status = mccmddbl('MUL',1,'Mul',mul,0.0D0,-1.0D0,mul) if( status.lt.0 ) return c --- ADD keyword add = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'ADD',1,add) endif status = mccmddbl('ADD',1,'Add',add,0.0D0,-1.0D0,add) if( status.lt.0 ) return c --- initialize destination image directory call zerow(64,ddir) c --- image version number ddir(2)=4 c --- sensor source number ivalue = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'SS',1,ivalue) endif status = mccmdint('SS',1,'Sensor Source',ivalue,0,999,ddir(3)) if( status.lt.0 ) return c --- image day if( CONTEXT ) then status = m0getdefaultargiyd( profile,'DAY',1,day) endif status = mccmdiyd('DAY',1,'Day',day,1970001,day,ddir(4)) if( status.lt.0 ) return c --- image time if( CONTEXT ) then status = m0getdefaultargihr( profile,'TIME',1,time) endif status = mccmdihr('TIM.E',1,'Time',time,0,235959,ddir(5)) if( status.lt.0 ) return c --- Upper Left line in satellite coordinates ivalue = 1 if( CONTEXT ) then status = m0getdefaultargint( profile,'LINE',1,ivalue) endif status=mccmdint('LIN.E',1,'Beginning Line',ivalue,0,-1,ddir(6)) if( status.lt.0 ) return NAV_lin = ddir(6) c --- Upper Left element in satellite coordinates ivalue = 1 if( CONTEXT ) then status = m0getdefaultargint( profile,'ELEM',1,ivalue) endif status=mccmdint('ELE.M',1,'Beginning Element',ivalue,0,-1,ddir(7)) if( status.lt.0 ) return NAV_ele = ddir(7) c --- number of image lines ddir(9) = (lines*lbu) / lbd c --- number of destination image elements (inforce multiple of 4 rule) ddir(10) = (elems*ebu) / ebd if( mod(ddir(10),4).ne.0 ) ddir(10) = (ddir(10)/4+1)*4 c --- bytes per element nbyte = 1 if( CONTEXT ) then status = m0getdefaultargint( profile,'NBYTE',1,nbyte) endif status=mccmdint('NBY.TE',1,'Bytes per Element',nbyte,1,4,nbyte) if( status.lt.0 ) return ddir(11) = nbyte c --- line resolution ivalue = 1 if( CONTEXT ) then status = m0getdefaultargint( profile,'LRES',1,ivalue) endif status=mccmdint('LRES',1,'Line resolution',ivalue,0,-1,ddir(12)) if( status.lt.0 ) return ddir(12) = (ddir(12)*lbd)/lbu if( ddir(12).eq.0 ) then call edest('Line resolutions less than 1 not allowed',0) return endif c --- element resolution ivalue = ddir(12) if( CONTEXT ) then status = m0getdefaultargint( profile,'ERES',1,ivalue) endif status=mccmdint('ERES',1,'Elem resolution',ivalue,0,-1,ddir(13)) if( status.lt.0 ) return ddir(13) = (ddir(13)*ebd)/ebu if( ddir(13).eq.0 ) then call edest('Element resolutions less than 1 not allowed',0) return endif c --- bands ndbands = 1 dbands(ndbands) = 1 if( CONTEXT ) then num = 0 status = m0getdefaultkeynumparms(profile,'BAND',num) if( num.gt.0 ) then do i = 1, num status = m0getdefaultargint( profile,'BAND',i,dbands(i)) if( status.lt.0 ) return enddo ndbands = num endif endif num = mccmdnum( 'BAN.D' ) if( num.gt.0 ) then ndbands = 0 do i = 1, num status = mccmdirng( & 'BAN.D', & i, & 'Band', & 0, & 0, & 1, & 1, & 64, & in_low, & in_high, & in_inc & ) if( status.lt.0 ) return do band = in_low, in_high, in_inc ndbands = ndbands + 1 dbands(ndbands) = band enddo enddo endif do i = 1, ndbands if( dbands(i).le.32 ) then ddir(19) = ddir(19) + 2**(dbands(i)-1) elseif( dbands(i).le.64 ) then ddir(20) = ddir(20) + 2**(dbands(i)-32) else call edest('Invalid BAND specified = ',dbands(i)) return endif enddo ddir(14) = ndbands c --- source band nsbands = ndbands do i = 1, nsbands sbands(i) = dbands(i) enddo sband = sbands(1) if( CONTEXT ) then status = m0getdefaultargint( profile,'SBAND',sband,sband) endif status=mccmdint('SBA.ND',1,'Source Band',sband,0,9999,sband) if( status.lt.0 ) return if(sband.ne.sbands(1)) then j = sband do i = 1, nsbands sbands(i) = j j = j+1 enddo endif c --- length of prefix if( ndbands.gt.1 ) then ddir(15) = ndbands/4*4 if( ddir(15).lt.ndbands ) ddir(15) = ddir(15)+4 ddir(51) = ddir(15) endif c --- memo cmemo = sfile if( CONTEXT ) then status = m0getdefaultargstr( profile,'MEMO',1,cmemo) endif status = mccmdstr('MEM.O',1,cmemo,cmemo) call movcw(cmemo,ddir(25)) c --- offset to the nav ivalue = 256 if( CONTEXT ) then status = m0getdefaultargint( profile,'NAV',1,ivalue) endif status = mccmdint('NAV',1,'NAV offset',ivalue,0,-1,ddir(35)) if( status.lt.0 ) return c --- offset to the cal ivalue = 768 if( CONTEXT ) then status = m0getdefaultargint( profile,'CAL',1,ivalue) endif status = mccmdint('CAL',1,'CAL offset',ivalue,0,-1,ddir(63)) if( status.lt.0 ) return c --- sensor type stype = 'PRD ' if( CONTEXT ) then status = m0getdefaultargstr( profile,'STYPE',1,stype) endif status = mccmdstr('STY.PE',1,stype,stype) if( status.lt.0 ) return ddir(52) = lit(stype(1:4)) c --- calibration type cvalue = 'BRIT' if( CONTEXT ) then status = m0getdefaultargstr( profile,'CTYPE',1,cvalue) endif status = mccmdstr('CTY.PE',1,cvalue,cvalue) if( status.lt.0 ) return ddir(53) = lit(cvalue(1:4)) c --- actual image start date and time ddir(46) = ddir(4) ddir(47) = mstime c --- BRIT to TEMP lookup table lookup_file = ' ' if( CONTEXT ) then status = m0getdefaultargstr( profile,'TABLE',1,lookup_file) endif status = mccmdstr('TAB.LE',1,lookup_file,lookup_file) if( status.lt.0 ) return if( lookup_file.ne.' ' ) then status = read_lookup_cal( & lookup_file, & ctype, & sub_type, & sbands(1), & table, table2 & ) if( status.lt.0 ) then len = len_trim( lookup_file ) call edest(lookup_file(1:len)//' does not exist',0) return endif endif c --- SCALE keyword if( nbpe.eq.2 ) then input_low = 0.D0 input_high = 65535.0D0 mx_value = 65535 else input_low = 0.D0 input_high = 255.0D0 mx_value = 255 endif if( nbyte.eq.2 ) then output_low = 0.0D0 output_high = 65535.0D0 else output_low = 0.0D0 output_high = 255.0D0 endif if( CONTEXT ) then status = m0getdefaultargdbl( profile,'SCALE',1,input_low) status = m0getdefaultargdbl( profile,'SCALE',2,input_high) status = m0getdefaultargdbl( profile,'SCALE',3,output_low) status = m0getdefaultargdbl( profile,'SCALE',4,output_high) endif status = & mccmddbl('SCA.LE',1,'Scale',input_low,0.0D0,-1.0D0,input_low) if( status.lt.0 ) return status = & mccmddbl('SCA.LE',2,'Scale',input_high,0.0D0,-1.0D0,input_high) if( status.lt.0 ) return status = & mccmddbl('SCA.LE',3,'Scale',output_low,0.0D0,-1.0D0,output_low) if( status.lt.0 ) return status = & mccmddbl('SCA.LE',4,'Scale',output_high,0.0D0,-1.0D0,output_high) if( status.lt.0 ) return scale = (output_high-output_low) / (input_high-input_low) c --- byte flipping of DATA swbyte_data = 'NO' if( CONTEXT ) then status = m0getdefaultargstr( profile,'SWBYTE',1,swbyte_data) endif status = mccmdstr('SWB.YTE',1,swbyte_data,swbyte_data) c --- zero the nav and cal arrays call zerow( MAXNAVSIZE, nbuf ) call zerow( MAXCALSIZE, calaux ) c --- checking for an external NAV file nfile = ' ' nfile2 = ' ' ntype = 'BIN' lalo_size = 0 LALOX2 = .FALSE. if( CONTEXT ) then status = m0getdefaultargstr( profile,'NFILE',1,nfile) status = m0getdefaultargstr( profile,'NFILE',2,ntype) endif status = mccmdstr('NFI.LE',1,nfile,nfile) if( status.lt.0 ) return status = mccmdstr('NFI.LE',2,ntype,ntype) if( status.lt.0 ) return c --- check if ntype is the name of the Longitude file status = 0 if( ntype(1:3).eq.'VMC' ) status = 1 if( ntype(1:3).eq.'BIN' ) status = 1 if( ntype(1:3).eq.'MCN' ) status = 1 if( ntype(1:3).eq.'ASC' ) status = 2 if( ntype(1:3).eq.'SAT' ) status = 2 if( ntype(1:3).eq.'ENV' ) status = 2 if( status.eq.0 ) then nfile2 = ntype status = mccmdstr('NFI.LE',3,ntype,ntype) if( status.lt.0 ) return LALOX2 = .TRUE. endif if( nfile.ne.' ' ) then c ------ Seperate Binary Latitude/Longitude files if( LALOX2 ) then c --------- what are the array dimensions? status = mccmdint('NFI.LE',4,'LALO lines', & lines,0,lines,lalox2_lines) status = mccmdint('NFI.LE',5,'LALO elems', & elems,0,elems,lalox2_elems) c --------- will it fit into the LALO size limit? lalox2_size = lalox2_lines*lalox2_elems if( lalox2_size.gt.MAX_AUXBLOCK_SIZE ) then call edest('STOP -- LALO NAV is too big',0) return endif c --------- what is the resolution of the lalo WRT image lalox2_lres = lines/lalox2_lines lalox2_eres = elems/lalox2_elems c --------- AUX block entry: LATITUDE ptr = 129 if( mccmdnum('AUX').ne.0 ) then status=mccmdint('AUX',1,'AUX offset',ptr*4,0,-1,ptr) ptr = (ptr-768)/4+1 endif calaux(ptr) = z'04030201' calaux(ptr+1) = 24 calaux(ptr+2) = 8 calaux(ptr+3) = lalox2_size*4 calaux(ptr+4) = lit('LATI') calaux(ptr+5) = lit('TUDE') call zerow(lalox2_size,calaux(ptr+6)) lat_offset = (ptr+5) * 4 lalox2_latptr = ptr+6 c --------- AUX block entry: LONGITUDE ptr = (ptr+6) + lalox2_size calaux(ptr) = z'04030201' calaux(ptr+1) = 28 calaux(ptr+2) = 9 calaux(ptr+3) = lalox2_size*4 calaux(ptr+4) = lit('LONG') calaux(ptr+5) = lit('ITUD') calaux(ptr+6) = lit('E ') call zerow(lalox2_size,calaux(ptr+7)) lon_offset = (ptr+6) * 4 lalox2_lonptr = ptr+7 c --------- read lalo and get min/max do i = 1, lalox2_lines do j = 1, lalox2_elems ptr = (i-1)*lalox2_elems + j-1 status = lwi( nfile, ptr,1,lalox2_lat) if( status.lt.0 ) return status = lwi( nfile2,ptr,1,lalox2_lon) call movb( 4, lalox2_lat, rlat, 0 ) call movb( 4, lalox2_lon, rlon, 0 ) if( ptr.eq.1 ) then minlat = rlat minlon = rlon maxlat = rlat maxlon = rlon else minlat = MIN( rlat, minlat ) minlon = MIN( rlon, minlon ) maxlat = MAX( rlat, maxlat ) maxlon = MAX( rlon, maxlon ) endif call movb( 4, rlat, calaux(lalox2_latptr), 0 ) call movb( 4, rlon, calaux(lalox2_lonptr), 0 ) lalox2_latptr = lalox2_latptr+1 lalox2_lonptr = lalox2_lonptr+1 enddo enddo c --------- compute the size of the LALO navigation (AUX) block lalo_size = (lalox2_size*4*2)+52 c --------- what are the image line and element of first nav point? status = mccmdint('NFI.LE',6,'LALO line', & 1,1,-1,lalox2_line) status = mccmdint('NFI.LE',7,'LALO elem', & 1,1,-1,lalox2_elem) c --------- fill navigation array nbuf(1) = lit('LALO') nbuf(2) = lit('AREA') nbuf(66) = lalox2_lines nbuf(67) = lalox2_elems call movb( 4, minlat, nbuf(68), 0 ) call movb( 4, minlon, nbuf(69), 0 ) call movb( 4, maxlat, nbuf(70), 0 ) call movb( 4, maxlon, nbuf(71), 0 ) nbuf(72) = ddir(12) * lalox2_lres nbuf(73) = ddir(13) * lalox2_eres nbuf(74) = 512 nbuf(75) = 512 + (lalox2_size*4) nbuf(76) = lalox2_line nbuf(77) = lalox2_elem nbuf(78) = (lalox2_size*4)*2+52 nbuf(79) = 24 nbuf(80) = (lalox2_size*4)+52 c ------ VMC Planetary navigation file elseif( ntype(1:3).eq.'VMC' ) then status = vmcread( & nfile, & imgorbit, & imgnumber, & imgchannel & ) if( status.lt.0 ) then call edest('FAILED -- No navigation for image',0) return endif c --------- mark this image as PLANetary navigation do i = 1,640 nbuf(i) = 0 enddo c ------ BINARY navigation file elseif( & ntype(1:3).eq.'BIN' .or. & ntype(1:3).eq.'MCN' & ) then navsize = m0arasize( 'NAV', ddir ) call sdest('right here ...size=',navsize) if( navsize.gt.0 ) then status = lbi( nfile, 0, navsize, nbuf ) if( status.lt.0 ) then call edest('Failed to read NAV file: '//nfile,0) return endif c ------------ byte flipping cvalue = 'YES' if( CONTEXT ) then status=m0getdefaultargstr(profile,'SWBYTE',2,cvalue) endif status = mccmdstr('SWB.YTE',2,cvalue,cvalue) if( cvalue(1:1).eq.'Y' ) then call swbyt4(nbuf(2),(navsize/4)-1) endif c------------- flag indicates McIDAS navigation block MCIDAS_NAV = .FALSE. if( ntype(1:3).eq.'MCN' ) MCIDAS_NAV = .TRUE. endif c ------ ASCII navigation file elseif( & ntype(1:3).eq.'ASC' .or. & ntype(1:3).eq.'SAT' .or. & ntype(1:3).eq.'ENV' & ) then call sdest('ASCII Navigation info: type='//ntype,0) c --------- get other ASC parameters lat_pos = 1 lat_len = 8 lon_pos = 10 lon_len = 9 skip = 0 nav_res = 1 if( CONTEXT ) then status=m0getdefaultargint( profile,'NFILE',3,lat_pos) status=m0getdefaultargint( profile,'NFILE',4,lat_len) status=m0getdefaultargint( profile,'NFILE',5,lon_pos) status=m0getdefaultargint( profile,'NFILE',6,lon_len) status=m0getdefaultargint( profile,'NFILE',7,skip) status=m0getdefaultargint( profile,'NFILE',8,nav_res) endif status=mccmdint('NFI.LE',3,'Lat_pos',lat_pos,0,-1,lat_pos) status=mccmdint('NFI.LE',4,'Lat_len',lat_len,0,-1,lat_len) status=mccmdint('NFI.LE',5,'Lon_pos',lon_pos,0,-1,lon_pos) status=mccmdint('NFI.LE',6,'Lon_len',lon_len,0,-1,lon_len) status=mccmdint('NFI.LE',7,'skip',skip,0,-1,skip) status=mccmdint('NFI.LE',8,'nav_res',nav_res,0,-1,nav_res) c --------- seperate lalo_res status = mccmdint('LALO',3,'LaLo Res',nav_res,0,-1,lalo_res) c --------- check navigation file name if( nfile(1:5).eq.'*.nav' ) then len = len_trim(sfile) nfile = sfile(1:len-3)//'nav' endif c --------- AUX block entry: LATITUDE ptr = 129 if( mccmdnum('AUX').ne.0 ) then status=mccmdint('AUX',1,'AUX offset',ptr*4,0,-1,ptr) ptr = (ptr-768)/4+1 endif calaux(ptr) = z'04030201' calaux(ptr+1) = 24 calaux(ptr+2) = 8 calaux(ptr+3) = ((lines/lalo_res)*(elems/lalo_res))*4 calaux(ptr+4) = lit('LATI') calaux(ptr+5) = lit('TUDE') call zerow((lines/lalo_res)*(elems/lalo_res),calaux(ptr+6)) lat_offset = (ptr+5) * 4 c --------- AUX block entry: LONGITUDE ptr = (ptr+6) + (lines/lalo_res)*(elems/lalo_res) calaux(ptr) = z'04030201' calaux(ptr+1) = 28 calaux(ptr+2) = 9 calaux(ptr+3) = ((lines/lalo_res)*(elems/lalo_res))*4 calaux(ptr+4) = lit('LONG') calaux(ptr+5) = lit('ITUD') calaux(ptr+6) = lit('E ') call zerow((lines/lalo_res)*(elems/lalo_res),calaux(ptr+7)) lon_offset = (ptr+6) * 4 c --------- Planetary navigation files (ENVI format) if(ntype(1:3).eq.'ENV'.or.ntype(1:3).eq.'SAT' ) then c ------------ make the name of the longitude file if( ntype(1:3).eq.'SAT' ) then status = index( nfile, 'latitude' ) if( status.le.0 ) then call edest('Invalid file name for ENVI nav',0) return else nfile2 = nfile(1:status-1)//'longitude.txt' endif else status = index( nfile, 'Latitude' ) if( status.le.0 ) then call edest('Invalid file name for ENVI nav',0) return else nfile2 = nfile(1:status-1)//'Longitude.txt' endif endif c ------------ get the ENVI lat/lon arrays call sdest(' ENVI Lat file='//nfile,0) call sdest(' ENVI Lon file='//nfile2,0) status = readENVInav( & nfile, & nfile2, & lines, & elems, & rlats, & rlons & ) if( status.lt.0 ) then call edest('FAILED - read of ENVI nav file=',0) return endif c ----------- fill the LALO arrays mislat=0 mislon=0 nlalo = lines*elems do i = 1, lines do j = 1, elems iptr = (i-1)*elems+j rlat = rlats(iptr) rlon = rlons(iptr) ii = i jj = j if( LINE_FLIP ) ii = lines-(i-1) if( ELEM_FLIP ) jj = elems-(j-1) optr = (ii-1)*elems+jj ptr = lat_offset + (optr-1)*4 call movb( 4, rlat, calaux, ptr ) ptr = lon_offset + (optr-1)*4 call movb( 4, rlon, calaux, ptr ) if( i.eq.1 ) then minlat = rlat minlon = rlon maxlat = rlat maxlon = rlon else if( rlat.ne.-999.0 ) then minlat = MIN(rlat,minlat) maxlat = MAX(rlat,maxlat) else mislat = mislat+1 endif if( rlon.ne.-999.0 ) then minlon = MIN(rlon,minlon) maxlon = MAX(rlon,maxlon) else mislon = mislon+1 endif endif enddo enddo call sdest( & ' Min Lat='//cff( dble(minlat),4 ),0) call sdest( & ' Max Lat='//cff( dble(maxlat),4 ),0) call sdest( & ' #Missing Lat=',mislat) call sdest( & ' Min Lon='//cff( dble(minlon),4 ),0) call sdest( & ' Max Lon='//cff( dble(maxlon),4 ),0) call sdest( & ' Missing Lon=',mislon) c --------- read the lats and lons else call sdest('LAT OFFSET = ',lat_offset) call sdest('LON OFFSET = ',lon_offset) nlalo = 0 101 status = read_ascii_lalo( & nfile, & skip, nav_res, & lalo_lines, lalo_elems, & lat_pos, & lat_len, & lon_pos, & lon_len, & rlat, & rlon & ) if(status.lt.0) then call edest('FAILED - read of ascii nav file=',0) return elseif(status.eq.1) then nlalo = nlalo+1 if( nlalo.gt.MAXLALOSIZE ) then call edest('FAILED -- MAXLALO overflow',0) return endif ptr = lat_offset + (nlalo-1)*4 call movb( 4, rlat, calaux, ptr ) ptr = lon_offset + (nlalo-1)*4 rlon = -1*rlon call movb( 4, rlon, calaux, ptr ) c -------------- min/max lat/lon if( nlalo.eq.1 ) then minlat = rlat minlon = rlon maxlat = rlat maxlon = rlon else minlat = MIN( rlat, minlat ) minlon = MIN( rlon, minlon ) maxlat = MAX( rlat, maxlat ) maxlon = MAX( rlon, maxlon ) endif goto 101 endif endif c --------- compute the size of the LALO navigation (AUX) block lalo_size = (nlalo*4*2)+52 c --------- fill navigation array call sdest(' Number of lalo points=',nlalo) call sdest(' size of LALO navigation=',lalo_size) nbuf(1) = lit('LALO') nbuf(2) = lit('AREA') nbuf(66) = lalo_lines/nav_res nbuf(67) = lalo_elems/nav_res call movb( 4, minlat, nbuf(68), 0 ) call movb( 4, minlon, nbuf(69), 0 ) call movb( 4, maxlat, nbuf(70), 0 ) call movb( 4, maxlon, nbuf(71), 0 ) nbuf(72) = ddir(12) * lalo_res nbuf(73) = ddir(13) * lalo_res nbuf(74) = 512 nbuf(75) = 512 + (nlalo*4) nbuf(76) = 1 nbuf(77) = 1 nbuf(78) = (nlalo*4)*2+52 nbuf(79) = 24 nbuf(80) = (nlalo*4)+52 endif endif c --- number of comment cards ncard = 1 c if( POST .and. ppcmd.ne.' ' ) ncard = 3 !Post Process command ddir(64) = ncard c --- cal array cvalue = clit( ddir(52) ) if( cvalue(1:3).eq.'PRD' ) then calaux(1) = lit('PROD') calaux(2) = nint(input_low) calaux(3) = nint(input_high) calaux(4) = nint(output_low) calaux(5) = nint(output_high) endif c --- offset to the aux ivalue = 0 jvalue = 0 if( lalo_size.gt.0 ) then ivalue = 1280 jvalue = lalo_size endif if( CONTEXT ) then status = m0getdefaultargint( profile,'AUX',1,ivalue) status = m0getdefaultargint( profile,'AUX',2,jvalue) endif status = mccmdint('AUX',1,'AUX offset',ivalue,0,-1,ddir(60)) status = mccmdint('AUX',2,'AUX size' ,jvalue,0,-1,ddir(61)) if( status.lt.0 ) return c --- offset to the data if( ivalue.gt.0 ) then call sdest('ddir(63)=',ddir(63)) call sdest('ddir(60)=',ddir(60)) call sdest('lalo_size=',lalo_size) if( ddir(60).gt.0 ) then ivalue = ddir(60) + lalo_size else ivalue = ddir(63) + lalo_size endif else ivalue = 1280 + lalo_size endif if( CONTEXT ) then status = m0getdefaultargint( profile,'DAT',1,ivalue) endif status = mccmdint('DAT.A',1,'DAT offset',ivalue,0,-1,ddir(34)) if( status.lt.0 ) return call sdest('DAT=',ddir(34)) c --- checking for a MAKNAV entry proj = NAV_proj if( CONTEXT ) then status = m0getdefaultargstr( profile, 'MAKNAV', 1, proj ) endif status = mccmdstr('MAKNAV', 1, proj, proj) if( proj.ne.' ' ) then status = maknav( profile, ddir, nbuf ) if( status.lt.0 ) then call edest('Failed to MAKNAV: '//proj ,0) return endif endif c --- dat arrays call zerow( MAXDATSIZE, ibuf ) call zerow( MAXDATSIZE, jbuf ) c --- construct sort conditions to be sent to the PUT server len = len_trim( ddataset_pos ) nsort = 1 sorts(nsort) = 'POS'// & ' '//ddataset_pos(1:len)// & ' '//ddataset_pos(1:len) c --- SPECIAL CODE FOR KALPANA CAL if( sub_type(1:4).eq.'KALP' ) then call sdest('LOAD KALPANA Calibration arrays',0) c ------ REFLECTANCE if( dbands(1).eq.1 ) then do i = 1,1024 calaux(i) = NINT( table2(i)*10000.0 ) enddo c ------ RADIANCE else do i = 1,1024 calaux(i) = INT( table2(i)*10000.0 ) calaux(i+1024) = INT( table(i)*10000.0 ) enddo endif endif c --- check for tracing if( mccmdnum( 'TRA.CE' ).ne. 0 ) then nsort = nsort+1 sorts(nsort) = 'TRACE=1' endif c --- initiate a transaction to write an image file status= mcaput( & ddataset_nam, & nsort, & sorts, & ddir, & nbuf, & calaux & ) if( status.lt.0 ) then call edest('Failed to initiate data write',0) return endif c --- initialize the min/max values input_min = 99999.0D0 input_max = -99999.0D0 output_min = 99999.0D0 output_max = -99999.0D0 c --- initialize number of missing data values nmiss = 0 c --- initialize number of blank data lines nblank = 0 c --- initialize the number of line sent nlines_sent = 0 c --- check for flip lines if( LINE_FLIP ) then lbeg = lines lend = 1 linc = -lbd else lbeg = 1 lend = lines linc = lbd endif c --- check for flip elemnets if( ELEM_FLIP ) then ebeg = elems eend = 1 einc = -ebd else ebeg = 1 eend = elems einc = ebd endif c --- pad the image lines if( pad_lines.gt.0 ) then c ------ make a line with all pad values do i = 1, MAXDATSIZE jbuf(i) = pad_value enddo call mpixel( MAXDATSIZE, 4, nbyte, jbuf) c ------ send to ADDE destination do line = 1, pad_lines if( nlines_sent.lt.lines ) then status = mcaout( jbuf ) if( status.lt.0 ) then call edest('Failed to write data to: '//ddataset_nam,0) return endif nlines_sent = nlines_sent+1 endif enddo endif do line = lbeg, lend, linc c ------ compute the starting byte of the specified line if( interleave(1:3).eq.'BIL' ) then ptr = ((line-1)*(elems*nbpe*n_file_bands)) + header ptr = ptr + ((sbands(1)-1)*elems*nbpe) else ptr = ((line-1)*(elems*nbpe))+header endif c ------ reading from the source file status = lbi( sfile, ptr, elems*nbpe, ibuf ) if( TEST ) then call sdest('TEST -- file='//sfile,0) call sdest('TEST -- ptr=',ptr) call sdest('TEST -- nbyte=',(elems*nbpe)) call sdest('TEST -- status=',status) endif c ------ FAILED reading image line ... if( status.lt.0 ) then c --------- ... pad fill the line array do i = 1, elems ibuf(i) = pad_value enddo c --------- ... increment number of blank lines nblank = nblank+1 c ------ SUCCESS reading image line ... else c --------- ... crack 1 byte into 4 bytes if( nbpe.eq.1 ) call mpixel(elems,1,4,ibuf) c --------- ... check for byte flipping if(swbyte_data(1:1).eq.'Y') then if(nbpe.eq.2 ) then call fbyte2(ibuf,elems) else call fbyte4(ibuf,elems) endif endif endif c ------ scale from source to destination j = 0 do i = ebeg, eend, einc c ------ get 1 data value if( nbpe.eq.1 ) then ivalue = ibuf(i) elseif( nbpe.eq.2 ) then ivalue = ibuf2(i) elseif( nbpe.eq.4 ) then ivalue = ibuf(i) endif c ------ FLOAT value, apply gain and offset if( format(1:1).eq.'R' .or. format(1:1).eq.'F' ) then call movb( 4, ivalue, rvalue, 0 ) c if( ivalue.ne.0 ) then c call edest('Ivalue = ',ivalue) c endif dvalue = dble(rvalue) c -------INTEGER value, apply gain and offset else dvalue = dble(ivalue) endif c ------ this will stop execution after first N elements if( TEST ) then cvalue = cfj( ivalue ) call edest('TEST -- i_value='//cvalue,0) cvalue = cff( dvalue, 4 ) call edest('TEST -- d_value='//cvalue,0) num_test = num_test+1 if( num_test.ge.max_test) return endif c ------ source min/max input_min = MIN( input_min, dvalue ) input_max = MAX( input_max, dvalue ) c ------ check for missing data value1 if( cmiss.ne.' ' ) then if( dvalue.eq.input_miss ) then dvalue = output_miss nmiss = nmiss+1 goto 10 endif endif c ------ check for missing data value2 if( cmiss2.ne.' ' ) then if( dvalue.eq.input_miss2 ) then dvalue = output_miss2 nmiss = nmiss+1 goto 10 endif endif c ------ NOT a missing data code ... apply the gain and offset dvalue = dvalue*gain+offset c ------ rescale the source to destination dvalue = ((dvalue-input_low)*scale)+output_low if( output_low.gt.output_high ) then dvalue = MAX(MIN(dvalue,output_low),output_high) else dvalue = MIN(MAX(dvalue,output_low),output_high) endif c ------ apply mul and add dvalue = dvalue*mul+add c ------ destination min/max 10 continue output_min = MIN( output_min, dvalue ) output_max = MAX( output_max, dvalue ) c ------ integerize the value ivalue = nint(dvalue) c ------ lookup table BRIT -> TEMP -> SSEC BRIT if( lookup_file.ne.' ' .and.sub_type(1:4).ne.'KALP' ) then jvalue = MIN( MAX(0, ivalue), mx_value) rtemp = table(jvalue+1) if( rtemp.lt.163.0 ) rtemp = 163.0 if( rtemp.gt.330.0 ) rtemp = 330.0 ivalue = gryscl( rtemp ) endif c ------ element blow up do k = 1, ebu j = j+1 jbuf(j) = ivalue enddo enddo call mpixel( j, 4, nbyte, jbuf) c ------ send to ADDE destination do k = 1, lbu status = mcaout( jbuf ) if( nlines_sent.lt.lines ) then if( status.lt.0 ) then call edest('Failed to write data to: '//ddataset_nam,0) return endif nlines_sent = nlines_sent+1 endif enddo enddo c --- add the comment card status = mcgetdaytime(day,time) cday = cfu(day) ctime = cfu(time) status = mccmd( cmd ) card(1:80) = cday(3:7)//' '//ctime(1:6)//' '//cmd(1:67) c --- post process command c if( POST .and. ppcmd.ne.' ' ) then c card(81:160) = cday(3:7)//' '//ctime(1:6)//' '//ppcmd(1:67) c call sdest( card(81:160),0) c card(161:240) =' '//ppcmd(68:134) c call sdest( card(161:240),0) c endif c --- write comment cards to image object call movb( ncard*80, card, icard, 0 ) status = mcacou( icard) c ***************** END OF IMAGE WRITE *************** c --- Post Process commands (these modify the image after it has been written) c --- 1: PRDUTIL to add a calibration to the image. This is only allowed if the c source type of the image is 'PRD'. c --- source type has to be PRD (Product) if( stype(1:3).eq.'PRD' ) then c ------ length of the dataset name len = len_trim( ddataset ) c ------ DRS sub-type has its own default parameters if( sub_type(1:3).eq.'DRS' ) then cal_type = 'REFL' cal_prodlo = '-320' cal_prodhi = '950' cal_britlo = '0' cal_brithi = '254' cal_unit = 'DBZ' cal_scale= '10' cal_miss = '255' else cal_type = ' ' cal_prodlo = ' ' cal_prodhi = ' ' cal_britlo = ' ' cal_brithi = ' ' cal_unit = ' ' cal_scale= ' ' cal_miss = ' ' endif c ------ get the cal_type if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',1,cal_type) endif status = mccmdstr('MAKCAL',1,cal_type,cal_type) c ------ get the cal_prodlo if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',2,cal_prodlo) endif status = mccmdstr('MAKCAL',2,cal_prodlo,cal_prodlo) c ------ get the cal_prodhi if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',3,cal_prodhi) endif status = mccmdstr('MAKCAL',3,cal_prodhi,cal_prodhi) c ------ get the cal_britlo if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',4,cal_britlo) endif status = mccmdstr('MAKCAL',4,cal_britlo,cal_britlo) c ------ get the cal_brithi if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',5,cal_brithi) endif status = mccmdstr('MAKCAL',5,cal_brithi,cal_brithi) c ------ get the cal_unit if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',6,cal_unit) endif status = mccmdstr('MAKCAL',6,cal_unit,cal_unit) c ------ get the cal_scale if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',7,cal_scale) endif status = mccmdstr('MAKCAL',7,cal_scale,cal_scale) c ------ get the cal_miss if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKCAL',8,cal_miss) endif status = mccmdstr('MAKCAL',8,cal_miss,cal_miss) c ------ create the PRDUTIL command if( cal_type.ne.' ' ) then cmd = 'PRDUTIL'// & ' ADD'// & ' '//ddataset(1:len)// & ' 2'// & ' '//cal_type// & ' '//cal_prodlo// & ' '//cal_prodhi// & ' '//cal_britlo// & ' '//cal_brithi// & ' '//cal_unit// & ' '//cal_scale// & ' '//cal_miss call bsquez( cmd ) call spout( cmd ) status = mcskeyin( cmd ) if( status.ne.0 ) then call edest('FAILED -- PRDUTIL: '//ddataset_nam,0) return endif endif endif c --- 2: Execute command in Post Process in PPCOM c if( POST .and. ppcmd.ne.' ' ) then c len = len_trim(ppcmd) c status = mcskeyin( ppcmd(1:len) ) c if( status.ne.0 ) then c call edest('FAILED -- Post Process of '//ddataset_nam,0) c return c endif c endif c --- set application status to SUCCESS call mccodeset( 0 ) c --- print the source min/max cvalue = cff( input_min, 3 ) cvalue2 = cfe( input_min,3 ) call sdest('Source minimum value = '//cvalue//cvalue2,0) cvalue = cff( input_max, 3 ) cvalue2 = cfe( input_max, 3 ) call sdest('Source maximum value = '//cvalue//cvalue2,0) cvalue = cff( output_min, 3 ) call sdest('Destination minimum value = '//cvalue,0) cvalue = cff( output_max, 3 ) call sdest('Destination maximum value = '//cvalue,0) cvalue = cfi( nmiss ) call sdest('Number of missing data values = '//cvalue,0) cvalue = cfi( nblank ) call sdest('Number of blank filled data lines = '//cvalue,0) c --- if source file is really a temp file -> delete it if( TMP_FILE ) then call lwd( tfile ) endif call edest('Done ...',0) return end integer function list_TIF(cfile,nlin,nele,nhed) implicit none c --- constants integer LREC parameter (LREC=12) c --- parameters character*(*) cfile integer nlin integer nele integer nhed integer ierr c --- external functions integer lbi c --- internal variables character*4 cunit integer file_byteorder integer machine_byteorder integer rat(2) integer unit integer nbit integer ibyte integer istat integer ifd integer nifd integer comprs integer iorder(2) integer itag(3) integer ival integer icurve(128) integer j integer*2 II integer*2 MM integer*2 curve(256) logical FLIP real od(256) real xres real yres real zscal real zxres real zyres integer*2 order integer*2 vsn integer point common/HEADER/ order, & vsn, & point integer*2 tag integer*2 type integer length integer voff common/RECORD/ tag, & type, & length, & voff equivalence (iorder, order) equivalence (itag, tag) equivalence (icurve, curve) data II/2HII/ data MM/2HMM/ unit = 2 xres = 2.54 yres = 2.54 nlin = 0 nele = 0 c --- get file byte order ibyte = 0 istat = lbi( cfile, ibyte, 8, iorder ) if( order.eq.II ) then file_byteorder = 0 call edest('File byte order is little endian',0) elseif( order.eq.MM ) then file_byteorder = 1 call edest('File byte order is big endian',0) else call edest('Invalid TIFF file',0) list_TIF = -1 return endif c --- get machine byte order call mcendn( machine_byteorder ) if( machine_byteorder.eq.0 ) then call edest('Machine byte order is little endian',0) else call edest('Machine byte order is big endian',0) endif c --- test if bytes need to be flipped FLIP = .FALSE. if( machine_byteorder.ne.file_byteorder ) FLIP=.TRUE. if( FLIP ) call edest('Bytes need to be flipped',0) if( FLIP ) call fbyte4( point, 1 ) ibyte = point istat = lbi( cfile, ibyte, 2, nifd) ibyte = ibyte+2 c --- scan the file tags call edest('Number of TIFF tags =',nifd) do 100 ifd = 1, nifd c --- read a tag line from the TIFF file istat = lbi( cfile, ibyte, LREC, itag ) ival = tag c --- gray response unit if( tag.eq.290 ) zscal = 10**voff c --- gray response curve if( tag.eq.291 ) then istat = lbi( cfile, voff, 512, icurve) do j = 1, 256 od(j) = curve(j)/zscal enddo endif c --- lines if( tag.eq.257 ) nlin = voff c --- elements if( tag.eq.256 ) nele = voff c --- strip offsets if( tag.eq.273 ) then if( length.eq.1 ) then nhed = voff else istat = lbi( cfile, voff, 4, nhed ) endif endif c --- number of bits per element if( tag.eq.258 ) then nbit = voff if( nbit.ne.8 ) then call edest('Input file must be an 8-bit image.',0) list_TIF = -1 return endif endif c --- compression if( tag.eq.259 ) then comprs = voff if( comprs.ne.1 ) then call edest('Image must be decompressed',0) list_TIF = -1 return endif endif c --- resolution unit if( tag.eq.296 ) unit = voff c --- xres if( tag.eq.282 ) then istat = lbi( cfile, voff, 8, rat ) xres = real(rat(1))/real(rat(2)) endif c --- yres if( tag.eq.283 ) then istat = lbi( cfile, voff, 8, rat ) yres = real(rat(1))/real(rat(2)) endif c --- increment the byte counter ibyte = ibyte+LREC 100 continue cunit='UNK ' if( unit.eq.2 ) then cunit = 'INCH' zxres = (10000.0) / (xres/2.54) zyres = (10000.0) / (yres/2.54) elseif( unit.EQ.3 ) then CUNIT = 'CM ' zxres = (10000.0) / (xres) zyres = (10000.0) / (yres) endif list_TIF = 0 return end integer function filedaytime( & filename, & namemask, & day, & time, & nfields, & fields & ) implicit none c --- constants integer NTYPES parameter (NTYPES = 9) integer MAX_FIELDS parameter (MAX_FIELDS = 20) c --- parameters character*(*) filename character*(*) namemask integer day integer time integer nfields character*80 fields(MAX_FIELDS) c --- external functions integer len_trim integer mcdmytocyd integer mcstrtoint c --- internal variables character*12 ctypes(NTYPES) character*80 cout character*256 name character*256 mask integer ccyy integer i, j, k integer status integer len integer len_name integer len_mask integer beg integer end integer pos integer bpos(NTYPES) integer epos(NTYPES) integer vals(NTYPES) data ctypes/ & '(CC)', ! 2 digit century & '(YY)', ! 2 digit year & '(DD)', ! 2 digit day of month & '(MN)', ! 2 digit month of year & '(HH)', ! 2 digit hour of day & '(MM)', ! 2 digit minute of hour & '(SS)', ! 2 digit second of minute & '(DDD)', ! 3 digit Julian day & '(*_)' ! variable length field ended by _ & / c --- internal copies of file name and mask name = filename mask = namemask c --- zero the number of fields nfields = 0 c --- lengths of name and mask len_name = len_trim( name ) len_mask = len_trim( mask ) c --- initialize the position arrays do i = 1, NTYPES bpos(i) = 0 epos(i) = 0 enddo c --- looking for place holders in the mask 1 beg = index( mask(1:len_mask), '(' ) end = index( mask(1:len_mask), ')' ) if( beg.gt.0 .and. end.gt.0 ) then do i = 1, NTYPES len = len_trim( ctypes(i) ) if( mask(beg:end).eq.ctypes(i)(1:len) ) then bpos(i) = beg epos(i) = end-2 do pos = bpos(i),epos(i) mask(pos:pos) = 'X' enddo pos = epos(i)+1 do j = pos, len_mask k = j+2 mask(j:j) = mask(k:k) enddo len_mask = len_trim( mask) goto 1 endif enddo endif c --- convert characters to integers do i = 1, NTYPES if( i.le.8 ) then if( bpos(i).ne.0) then status = mcstrtoint( name(bpos(i):epos(i)), vals(i) ) if( status.lt.0 ) then filedaytime = -1 return endif else vals(i) = 0 endif endif enddo c --- make julian day ccyy = (vals(1)*100) + (vals(2)) if( vals(8).gt.0 ) then day = (ccyy*1000) + (vals(8)) elseif( vals(3).gt.0 .and. vals(4).gt.0 ) then status = mcdmytocyd( vals(3), vals(4), ccyy, day ) if( status.lt.0 ) then filedaytime = -2 return endif else day = 0 endif c --- make a time time = (vals(5)*10000) + (vals(6)*100) + (vals(7)) filedaytime = 0 return end integer function filemask( & filename, & namemask, & day, & time, & nfields, & fields & ) implicit none c --- constants integer NTYPES parameter (NTYPES = 8) integer MAX_FIELDS parameter (MAX_FIELDS = 20) c --- parameters character*(*) filename character*(*) namemask integer day integer time integer nfields character*80 fields(MAX_FIELDS) c --- external functions integer len_trim integer mcdmytocyd integer mcstrtoint c --- internal variables character*1 cval character*12 ctypes(NTYPES) character*80 cout character*80 string character*256 name character*256 workname character*256 mask character*256 oldmask character*256 newmask integer ccyy integer i, j, k integer status integer len integer len_name integer len_mask integer len_workname integer len_oldmask integer len_newmask integer beg integer end integer pos integer bpos(NTYPES) integer epos(NTYPES) integer vals(NTYPES) data ctypes/ & '(CC)', ! 2 digit century & '(YY)', ! 2 digit year & '(DD)', ! 2 digit day of month & '(MN)', ! 2 digit month of year & '(HH)', ! 2 digit hour of day & '(MM)', ! 2 digit minute of hour & '(SS)', ! 2 digit second of minute & '(DDD)' ! 3 digit Julian day & / c --- internal copies of file name and mask name = filename mask = namemask c --- zero the number of fields nfields = 0 c --- lengths of name and mask len_name = len_trim( name ) len_mask = len_trim( mask ) c --- initialize the position arrays do i = 1, NTYPES bpos(i) = 0 epos(i) = 0 enddo workname = name len_workname = len_name oldmask = mask len_oldmask = len_mask newmask = ' ' len_newmask = 1 100 continue beg = index( oldmask, '*' ) if( beg.gt.0 ) then end = beg+1 cval = oldmask(end:end) pos = index( workname, cval(1:1) ) if( pos.le.0 ) then filemask = -1 return else c --------- save the field text nfields = nfields+1 fields(nfields) = workname(1:pos-1) c --------- make a string of X's string = ' ' do i = 1, pos string(i:i) = 'X' enddo c --------- reduce the old mask oldmask = oldmask(end+1:len_oldmask) len_oldmask = len_trim( oldmask ) c --------- reduce the work name workname = workname(pos+1:len_workname) len_workname = len_trim( workname ) c --------- create the new mask newmask = newmask(1:len_newmask)// & string(1:pos) call bsquez( newmask ) len_newmask = len_trim( newmask ) goto 100 endif else if( newmask.ne.' ' ) then mask = newmask(1:len_newmask)//oldmask(1:len_oldmask) len_mask = len_trim( mask ) endif endif c --- looking for place holders in the mask 200 continue beg = index( mask(1:len_mask), '(' ) end = index( mask(1:len_mask), ')' ) if( beg.gt.0 .and. end.gt.0 ) then do i = 1, NTYPES len = len_trim( ctypes(i) ) if( mask(beg:end).eq.ctypes(i)(1:len) ) then bpos(i) = beg epos(i) = end-2 do pos = bpos(i),epos(i) mask(pos:pos) = 'X' enddo pos = epos(i)+1 do j = pos, len_mask k = j+2 mask(j:j) = mask(k:k) enddo len_mask = len_trim( mask) goto 200 endif enddo endif c --- convert characters to integers do i = 1, NTYPES if( i.le.8 ) then if( bpos(i).ne.0) then status = mcstrtoint( name(bpos(i):epos(i)), vals(i) ) if( status.lt.0 ) then filemask = -2 return endif else vals(i) = 0 endif endif enddo c --- make julian day if( vals(1).eq.0 ) vals(1) = 20 ccyy = (vals(1)*100) + (vals(2)) if( vals(8).gt.0 ) then day = (ccyy*1000) + (vals(8)) elseif( vals(3).gt.0 .and. vals(4).gt.0 ) then status = mcdmytocyd( vals(3), vals(4), ccyy, day ) if( status.lt.0 ) then filemask = -3 return endif else day = 0 endif c --- make a time time = (vals(5)*10000) + (vals(6)*100) + (vals(7)) filemask = 0 return end integer function maknav( profile, dir, nav ) implicit none c --- constants integer NPLANETS parameter (NPLANETS=8) c --- parameters character*80 profile integer dir(*) integer nav(*) c --- external functions integer lit integer ilalo integer isdgch integer isatnum integer mccmdnum integer mccmdstr integer mccmdint integer mccmddll integer mccmddbl integer mccydtodmy integer mcstrtoint integer mcstrtodbl integer m0getdefaultargstr integer m0getdefaultargint integer m0getdefaultargdbl integer m0getdefaultargdll integer mcaput integer nchars integer numsat integer volnam c --- internal variables character*12 pole character*12 planet character*12 planets(NPLANETS) character*12 cname character*12 proj character*80 cfile character*80 cline character*80 cpath double precision dpar double precision decc(NPLANETS) double precision lin_0 double precision ele_0 double precision slat double precision slon double precision sscale double precision depoctime double precision dincline double precision drtascen double precision deccentr double precision dperigee double precision danomaly double precision dmmotion double precision depocrev integer i integer ib integer ie integer iday integer imon integer iyr integer nsec integer ihr integer imin integer isec integer indx integer irad integer iecc integer iret1 integer iret2 integer iepocdate integer ipowlat integer ipowlon integer ipowdlin integer ipowdele integer ipowrad integer ipowecc integer nch integer status logical CONTEXT real t2 real gm4pi2 real r real rad(NPLANETS) c --- navigation common character*4 NAV_proj integer NAV_lin integer NAV_ele double precision NAV_lat double precision NAV_lon double precision NAV_deglin double precision NAV_degele common /NAVCOM/ & NAV_proj, & NAV_lin, & NAV_ele, & NAV_lat, & NAV_lon, & NAV_deglin, & NAV_degele c --- PLANETARY INFORMATION data planets/ & 'MERCURY', & 'VENUS', & 'EARTH', & 'MARS', & 'JUPITER', & 'SATURN', & 'URANUS', & 'NEPTUNE' & / data rad/ & 2440., & 6120., & 6378.388, & 3393.5, & 71400., & 60330., & 25900., & 24750. & / data decc/ & 0.D0, & 0.D0, & 0.081992D0, & 0.1333386D0, & 0.35412D0, & 0.41019D0, & 0.2086D0, & 0.2D0 & / c --- set the context table logical CONTEXT = .FALSE. if( profile.ne.' ' ) CONTEXT = .TRUE. c --- get the planet constants planet = 'EARTH' if( CONTEXT ) then status = m0getdefaultargstr( profile,'PLANET',1,planet) endif status = mccmdstr('PLA.NET',1,planet,planet) if( status.lt.0 ) return do i = 1,NPLANETS if( planet.eq.planets(i)) then ipowrad = 3 irad = nint(rad(i)*10.**ipowrad) ipowecc = 6 iecc = idnint(decc(i)*10.D0**ipowecc) goto 10 endif enddo call edest('MAKNAV FAILED: Invalid Planet name='//planet,0) maknav = -1 return 10 continue c --- get the projection proj = NAV_proj if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKNAV',1,proj) endif status = mccmdstr('MAKNAV',1,proj,proj) if( status.lt.0 ) then maknav = -2 return endif c --- Polar Sterographic Projection if( proj(1:2).eq.'PS' ) then maknav = -3 nav(1) = lit('PS ') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'NPL',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'NPE',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 60.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = mccmddll('MAKNAV',4,'SLAT',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(4) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = mccmddbl('MAKNAV',6,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return if( dpar.lt.1000.D0 ) then nav(5) = int(dpar*1000.0D0) else nav(5) = int(dpar*.000570D0) endif dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'QLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(6) = ilalo(real(dpar)) nav(7) = irad nav(8) = iecc if( planet(1:5).ne.'EARTH') then nav(9) = -1 nav(10) = -1 else nav(9) = 0 nav(10) = 0 endif pole = 'N' if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKNAV',7,pole) endif status = mccmdstr('MAKNAV',7,pole,pole) if( status.lt.0 ) return if( pole.eq.'S') then nav(11) = -900000 else nav(11) = 900000 endif c --- Lambert Conformal projection elseif( proj(1:4).eq.'LAMB' ) then maknav = -3 nav(1) = lit('LAMB') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'NPL',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'NPE',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 30.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'SLAT1',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(4) = ilalo(real(dpar)) dpar = 60.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'SLAT2',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(5) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',7,dpar) endif status = & mccmddbl('MAKNAV',7,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return if( dpar.lt.1000.D0 ) then nav(6) = int(dpar*1000.0D0) else nav(6) = int(dpar*.000570D0) endif dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',6,dpar) endif status = & mccmddll('MAKNAV',6,'QLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(7) = ilalo(real(dpar)) nav(8) = irad nav(9) = iecc if( planet(1:5).ne.'EARTH') then nav(10) = -1 nav(11) = -1 else nav(10) = 0 nav(11) = 0 endif pole = 'N' if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKNAV',8,pole) endif status = mccmdstr('MAKNAV',8,pole,pole) if( status.lt.0 ) return if( pole.eq.'S') then nav(12) = -900000 else nav(12) = 900000 endif c --- Mercator Projection elseif(proj(1:4).eq.'MERC' ) then maknav = -3 nav(1) = lit('MERC') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'EQL',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'EQE',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'SLAT1',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(4) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = & mccmddbl('MAKNAV',6,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return if( dpar.lt.1000.D0 ) then nav(5) = int(dpar*1000.0D0) else nav(5) = int(dpar*.000570D0) endif dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'QLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(6) = ilalo(real(dpar)) nav(7) = irad nav(8) = iecc if( planet(1:5).ne.'EARTH') then nav(9) = -1 nav(10) = -1 else nav(9) = 0 nav(10) = 0 endif c --- Radar projection elseif( proj(1:4).eq.'RADR' ) then maknav = -3 nav(1) = lit('RADR') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'CENL',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'CENE',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'CLAT',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(4) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'CLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(5) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = & mccmddbl('MAKNAV',6,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(6) = int( dpar*1000.0D0 ) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',7,dpar) endif status = & mccmddbl('MAKNAV',7,'ROT',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(7) = int( dpar*1000.0D0 ) c --- Sinusoidal projection elseif( proj(1:3).eq.'SIN' ) then maknav = -3 nav(1) = lit('SIN ') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'CENL',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'CENE',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'CLAT',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(4) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'CLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(5) = ilalo(real(dpar)) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = & mccmddbl('MAKNAV',6,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(6) = int( dpar*1000.0D0 ) nav(7) = irad nav(8) = iecc if( planet(1:5).ne.'EARTH') then nav(9) = -1 nav(10) = -1 else nav(9) = 0 nav(10) = 0 endif c --- Rectilinear projection elseif( proj(1:4).eq.'RECT' ) then maknav = -3 nav(1) = lit('RECT') nav(2) = NAV_lin if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'LIN',nav(2),0,-1,nav(2)) if( status.lt.0 ) return ipowlat = 4 dpar = NAV_lat if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',3,dpar) endif status = mccmddll('MAKNAV',3,'LAT',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(3)=idnint(10.**ipowlat*dpar) nav(4) = NAV_ele if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',4,nav(4)) endif status = mccmdint('MAKNAV',4,'ELE',nav(4),0,-1,nav(4)) if( status.lt.0 ) return ipowlon = 4 dpar = -NAV_lon if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'LON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(5)=idnint(10.**ipowLON*dpar) ipowdlin = 6 dpar = NAV_deglin if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = & mccmddbl('MAKNAV',6,'DEGLIN',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(6) = idnint(10.**ipowdlin*dpar) ipowdele = 6 dpar = NAV_degele if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',7,dpar) endif status = & mccmddbl('MAKNAV',7,'DEGELE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(7) = idnint(10.**ipowdele*dpar) nav(8) = irad nav(9) = iecc if( planet(1:5).ne.'EARTH') then nav(10) = -1 nav(11) = -1 else nav(10) = 0 nav(11) = 0 endif nav(12) = ipowlat nav(13) = ipowlon nav(14) = ipowdlin nav(15) = ipowdele nav(16) = ipowrad nav(17) = ipowecc c --- Orthographic projection elseif( proj(1:4).eq.'ORTH') then maknav = -3 nav(1) = lit('ORTH') nav(2) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',2,nav(2)) endif status = mccmdint('MAKNAV',2,'MAKNAV(2)',nav(2),0,-1,nav(2)) if( status.lt.0 ) return nav(3) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(3)) endif status = mccmdint('MAKNAV',3,'MAKNAV(3)',nav(3),0,-1,nav(3)) if( status.lt.0 ) return dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',4,dpar) endif status = & mccmddbl('MAKNAV',4,'MAKNAV(4)',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(4) = idnint( dpar*1000.0D0 ) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',5,dpar) endif status = & mccmddbl('MAKNAV',5,'MAKNAV(5)',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(5) = idnint( dpar*1000.0D0 ) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = & mccmddbl('MAKNAV',6,'MAKNAV(6)',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(6) = idnint( dpar*100.0D0 ) nav(7) = irad nav(8) = iecc dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',7,dpar) endif status = & mccmddbl('MAKNAV',7,'MAKNAV(7)',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(9) = idnint( dpar*100.0D0 ) if( planet(1:5).ne.'EARTH') then nav(10) = -1 nav(11) = -1 else nav(10) = 0 nav(11) = 0 endif c --- Tangent Cone projection elseif( proj(1:4).eq.'TANC' ) then maknav = -3 nav(1) = lit('TANC') dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',2,dpar) endif status = mccmddbl('MAKNAV',2,'NPL',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(2) = nint(10000.0D0*dpar) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',3,dpar) endif status = mccmddbl('MAKNAV',3,'NPE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(3) = nint(10000.0D0*dpar) dpar = 1.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',6,dpar) endif status = mccmddbl('MAKNAV',6,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(4) = NINT(10000.0D0*dpar) dpar = 25.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'SLAT',dpar,-90.0D0,90.0D0,dpar) if( status.lt.0 ) return nav(5) = nint(10000.0D0*dpar) dpar = 25.0D0 if( CONTEXT ) then status = m0getdefaultargdll( profile,'MAKNAV',5,dpar) endif status = & mccmddll('MAKNAV',5,'SLON',dpar,-180.0D0,180.0D0,dpar) if( status.lt.0 ) return nav(6) = nint(10000.0D0*dpar) c --- Mollweide projection elseif( proj(1:4).eq.'MOLL') then maknav = -3 nav(1) = lit('MOLL') dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',2,dpar) endif status = mccmddbl('MAKNAV',2,'CENL',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(2) = nint(dpar) dpar = 0.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',3,dpar) endif status = mccmddbl('MAKNAV',3,'CENE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(3) = nint(dpar) dpar = 1.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',5,dpar) endif status = mccmddbl('MAKNAV',5,'SIZE',dpar,0.0D0,-1.0D0,dpar) if( status.lt.0 ) return nav(4) = nint(dpar) dpar = 25.0D0 if( CONTEXT ) then status = m0getdefaultargdbl( profile,'MAKNAV',4,dpar) endif status = & mccmddll('MAKNAV',4,'SLON',dpar,-180.0D0,180.0D0,dpar) nav(5) = nint(dpar) nav(6) = 0 nav(7) = irad nav(8) = iecc if( planet(1:5).ne.'EARTH') then nav(9) = -1 nav(10) = -1 else nav(9) = 0 nav(10) = 0 endif nav(11) = 0 nav(12) = 0 nav(13) = 0 nav(14) = 0 nav(15) = lit('KMPP') c --- TIROs projection elseif( proj(1:4).eq.'NOAA') then maknav = -3 c ------ Keplerian element file cfile = ' ' if( CONTEXT ) then status = m0getdefaultargstr( profile,'MAKNAV',2,cfile) endif status = mccmdstr('MAKNAV',2,cfile,cfile) if( status.lt.0 ) then call edest('MAKNAV FAILED: NOAA file name: '//cfile,0) return else status = volnam( cfile, cpath ) if( status.le.0 ) then call edest('MAKNAV FAILED: NOAA file name: '//cfile,0) return endif endif c ------ Name the navigation type nav(1) = lit('TIRO') c ------ Get the nominal satellite DAY [YYDDD] nav(2) =-1 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',3,nav(2)) endif status = mccmdint( 'MAKNAV',3,'DAY',nav(2),1,0,nav(2) ) if( status.lt.0 .or. status.eq.100 ) then call edest('MAKNAV FAILED: DAY for NOAA projection',0) return endif dir(4) = nav(2) c ------ Get the nominal satellite TIME [HHMMSS] nav(3) =-1 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',4,nav(3)) endif status = mccmdint('MAKNAV',4,'TIME',nav(3),1,0,nav(3) ) if( status.lt.0 .or. status.eq.100 ) then call edest('MAKNAV FAILED: TIME for NOAA projection',0) return endif dir(5) = nav(3) c ------ Orbit type (always 1) nav(4) = 1 c ------ Open the Keplerian element file OPEN(10,FILE=cpath,STATUS='OLD',IOSTAT=status) if( status.lt.0 ) then call edest('MAKNAV FAILED: OPEN element file: '//cpath,0) return endif c ------ Read the Keplerian element file cline = ' ' READ(10,'(80A)') cline nch = nchars(cline,ib,ie) c ------ Check to see if satellite name occupies the first line if( isdgch(cline(1:1)).lt.0 ) then nch = nchars(cline,ib,ie) if( nch.le.0 ) then call edest('MAKNAV FAILED: NOAA satellite name: '//cpath,0) return endif cname = cline(ib:ie) else cname = 'NOAA-14' nch = 7 endif indx = index( cname(1:nch),' ' ) if( indx.gt.0 ) cname(indx:indx) = '-' dir(3) = numsat( cname(1:nch) ) nav(2) = 100000 * dir(3) + nav(2) c ------ Make sure that the line number of element data is '1' 1 read(10,'(80A)') cline if( cline(1:1).ne.'1' ) then call edest('MAKNAV FAILED: extract values from '//cpath,0) return endif c ------ Satellite number status = mcstrtoint( cline(3:7), isatnum ) c ------ Epoch date [YYMMDD] status = mcstrtoint( cline(19:23), iepocdate ) status = mccydtodmy( iepocdate, iday, imon, iyr ) nav(5) = 10000*MOD(iyr,100) + 100*imon + iday c ------ Epoch time [HHMMSS] status = mcstrtodbl( cline(24:32), dpar ) nsec = 86400 * dpar ihr = nsec / 3600 nsec = nsec - 3600 * ihr imin = nsec / 60 isec = nsec - 60 * imin nav(6) = 10000*ihr + 100*imin + isec c ------ Make sure that the line number of element data is '2' 2 READ(10,'(80A)') cline if( cline(1:1).ne.'2' ) then call edest('MAKNAV FAILED: extract values from '//cpath,0) return endif c ------ close the navigation file CLOSE(10) c ------ Inclination [deg * 1000] status = mcstrtodbl( cline(9:16), dpar ) nav(9) = idnint(dpar*1000.0D0) c ------ Right ascension of ascending node [deg * 1000] status = mcstrtodbl( cline(18:25), dpar ) nav(12) = idnint(dpar*1000.0D0) c ------ Eccentricity [ * 1000000] cline(26:26) = '.' status = mcstrtodbl( cline(26:33), dpar ) nav(8) = idnint(dpar*1000000.0D0) c ------ Argument of perigee [deg * 1000] status = mcstrtodbl( cline(35:42), dpar ) nav(11) = idnint(dpar*1000.0D0) c ------ Mean anomaly [deg * 1000] status = mcstrtodbl( cline(44:51), dpar ) nav(10) = idnint(dpar*1000.0D0) c ------ Semi-major axis [km * 100] c c G * M * T**2 / (4 * pi**2) = r**3 c c G - 6.67E-11 nt*m**2/kg**2 -- Gravitational constant c M - 5.9760194E24 kg -- Mass of Earth c T - s -- orbital period c pi - 3.1415926 -- pi c r - m -- length of semi-major axis c c G * M / (4 * pi**2) - 1.0108462E13 c c dmmotion -- revolutions per day status = mcstrtodbl( cline(53:63), dpar ) T2 = (86400 / dpar)**2 GM4PI2 = 6.67E-11 * 5.9760194E24 / (4 * (3.1415926)**2) r = (GM4PI2 * T2) ** (1.0/3.0) nav(7) = 0.1 * r c ------ Number of samples per line (elements) nav(13) = dir(10) c ------ Angular increment between samples [deg * 1000000] nav(14) = 54128 c ------ Fraction of a second at epoch time (fixed to be 560 for NOAA satellites nav(15) = 560 c ------ Direction of satellite pass [-1->descending, +1->ascending] nav(46) =1 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',6,nav(46)) endif status = mccmdint('MAKNAV',6,'A/D',nav(46),-1,1,nav(46) ) if( status.lt.0 ) then call edest('MAKNAV FAILED: A/D for NOAA projection',0) return endif c ------ Number of earliest line to navigate [image coordinates] nav(47) = dir(6) c ------ Time at the start of the earliest line [msec from start of DAY] nav(48) =1 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',5,nav(48)) endif status = mccmdint('MAKNAV',5,'ITIME',nav(48),1,0,nav(48) ) if( status.lt.0 .or. status.eq.100 ) then call edest('MAKNAV FAILED: TMSEC for NOAA proection',0) return endif c ------ Time interval between lines [msec] nav(49) = 167 c ------ Flag for inverted image [0->normal, 1->inverted] nav(50) = 0 if( CONTEXT ) then status = m0getdefaultargint( profile,'MAKNAV',7,nav(50)) endif status = mccmdint('MAKNAV',7,'FLIP',nav(50),0,1,nav(50) ) if( status.lt.0 ) then call edest('MAKNAV FAILED: FLIP for NOAA projection',0) return endif c ------ Number of lines in the image nav(51) = dir(9) c ------ Number of samples per line (elements) nav(52) = dir(10) c ------ Time interval between lines [usec] nav(53) = 166667 c ------ Time between individual pixels [*100000000] nav(54) = 8138 c --- Invalid projection else maknav = -4 call edest('MAKNAV FAILED: Invalid projection',0) return endif maknav = 0 return end integer function numsat(cname) implicit none c --- parameters character*(*) cname c --- external functions integer lwi integer mcstrtoint c --- internal variables character*12 cfile character*80 ctext integer status integer strtwd integer indx integer line(20) data cfile/'SATANNOT'/ equivalence (line, ctext) strtwd = 0 10 status = lwi(cfile,strtwd,20,line) if( status.ne.0 ) then if( strtwd.eq.0) THEN call edest('MAKNAV FAILED: Unable to read SATANNOT',0) endif numsat = -1 return endif indx = index( ctext, cname ) if( indx.le.0 ) then strtwd = strtwd+20 goto 10 endif status= mcstrtoint(ctext(29:32),numsat) return end integer function list_NSSL( & sfile, & lines, & elems, & bands, & header, & cyd, & hms & ) implicit none c --- parameters character*(*) sfile integer lines integer elems integer bands integer header integer cyd integer hms c --- external functions character*12 cff integer lbi integer mcdmytocyd c --- internal variables character*12 cvalue integer i integer ptr integer status c --- NSSLCOM holds all values decoded from the header integer year integer month integer day integer hour integer minute integer second integer nlines integer nelems integer nlevels character*4 projection integer map_scale integer slat1 integer slat2 integer slon integer nw_lon integer nw_lat integer ll_scale integer lon_res integer lat_res integer res_scale integer height integer height_scale integer spare(10) character*20 name character*6 unit integer factor integer missing integer nradars character*4 radars(100) common /NSSLCOM/ & year, & month, & day, & hour, & minute, & second, & nelems, & nlines, & nlevels, & projection, & map_scale, & slat1, & slat2, & slon, & nw_lon, & nw_lat, & ll_scale, & lon_res, & lat_res, & res_scale, & height, & height_scale, & spare, & name, & unit, & factor, & missing, & nradars, & radars c --- navigation common character*4 NAV_proj integer NAV_lin integer NAV_ele double precision NAV_lat double precision NAV_lon double precision NAV_deglin double precision NAV_degele common /NAVCOM/ & NAV_proj, & NAV_lin, & NAV_ele, & NAV_lat, & NAV_lon, & NAV_deglin, & NAV_degele list_NSSL = -1 c --- read the image header status = lbi( sfile, 0, 154, year ) if( status.lt.0 ) return c --- second read is used to get over byte boundary problem status = lbi( sfile, 154, 12, factor ) if( status.lt.0 ) return c --- read the radar ids ptr = 166 do i = 1, nradars status = lbi( sfile, ptr, 4, radars(i) ) if( status.lt.0 ) return ptr = ptr+4 enddo list_NSSL = -2 c --- convert year, month and day to cyd status = mcdmytocyd( day, month, year, cyd ) if( status.lt.0 ) return c --- convert hour, minute, second to hms hms = (hour*10000) + (minute*100) + second c --- Number of image lines lines = nlines c --- Number of image elements elems = nelems c --- Number of image bands bands = nlevels c --- Size of header header = 166 + (nradars*4) c --- setup navigation common if( projection.eq.'LL ' ) then NAV_proj = 'RECT' NAV_lat = dble(nw_lat) / dble(ll_scale) NAV_lon = dble(nw_lon) / dble(ll_scale) NAV_deglin = dble(lat_res) / dble(res_scale) NAV_degele = dble(lon_res) / dble(res_scale) call sdest(' ',0) call sdest('NSSL: Default Navigation Parameters:',0) call sdest(' Projection = '//NAV_proj,0) cvalue = cff( NAV_lat, 4 ) call bsquez( cvalue ) call sdest(' NW Latitude = '//cvalue,0) cvalue = cff( NAV_lon, 4 ) call bsquez( cvalue ) call sdest(' NW Longitude = '//cvalue,0) cvalue = cff( NAV_deglin, 4 ) call bsquez( cvalue ) call sdest(' Degrees Per Line = '//cvalue,0) cvalue = cff( NAV_degele, 4 ) call bsquez( cvalue ) call sdest(' Degrees Per Elem = '//cvalue,0) call sdest(' ',0) endif list_NSSL = 0 return end integer function list_ENVI( & sheader, & lines, & elems, & bands, & header, & nbpe, & format, & interleave, & missing, & missing2, & cal_scale, & cal_offset & ) implicit none c --- parameters character*(*) sheader integer lines integer elems integer bands integer header integer nbpe character*(*) format character*(*) interleave character*(12) missing, missing2 double precision cal_scale double precision cal_offset c --- external functions character*12 cff integer lbi integer mcdmytocyd integer mcstrtoint integer mcstrtodbl integer getline c --- internal variables character*12 form character*12 cvalue character*80 string character*256 cline character*256 sfile integer length integer stat integer i integer code integer ptr integer status sfile = sheader lines = 0 elems = 0 bands = 0 header = 0 interleave = 'NON' format = 'INTEGER' nbpe = 1 missing = ' ' missing2 = ' ' cal_scale = 1.0D0 cal_offset = 0.0D0 1 continue status = getline( & sfile, & length, & stat, & cline & ) if ( status.lt.0 ) then call edest('FAILED -- unable to read header file',0) list_ENVI = -1 return elseif( status.eq.1 ) then c ------ check for "lines" ptr = index( cline(1:length), 'lines = ') if( ptr.gt.0 ) then string = cline((ptr+7):length) call bsquez(string) status = mcstrtoint( string, lines) endif c ------ check for "samples" ptr = index( cline(1:length), 'samples = ') if( ptr.gt.0 ) then string = cline((ptr+9):length) call bsquez(string) status = mcstrtoint( string, elems) endif c ------ check for "bands" ptr = index( cline(1:length), 'bands = ') if( ptr.gt.0 ) then string = cline((ptr+7):length) call bsquez(string) status = mcstrtoint( string, bands) endif c ------ check for "header offset" ptr = index( cline(1:length), 'header offset = ') if( ptr.gt.0 ) then string = cline((ptr+15):length) call bsquez(string) status = mcstrtoint( string, header) endif c ------ check for "data type" ptr = index( cline(1:length), 'data type = ') if( ptr.gt.0 ) then string = cline((ptr+11):length) call bsquez(string) status = mcstrtoint( string, code) if( code.eq.1 ) then format = 'INTEGER' nbpe = 1 elseif( code.eq.2 ) then format = 'INTEGER' nbpe = 2 elseif( code.eq.3 ) then format = 'INTEGER' nbpe = 2 elseif( code.eq.4 ) then format = 'FLOAT' nbpe = 4 elseif( code.eq.5 ) then format = 'FLOAT' nbpe = 8 endif endif c ------ check for "interleave" ptr = index( cline(1:length), 'interleave = ') if( ptr.gt.0 ) then form = cline((ptr+12):length) call bsquez(form) call mcupcase( form ) interleave = form endif c ------ check for "non-data" ptr = index( cline(1:length), 'non-data = ') if( ptr.gt.0 ) then missing = cline((ptr+10):length) call bsquez(missing) if( missing(1:2).eq.'-1' ) missing = ' ' endif c ------ check for "non-data2" ptr = index( cline(1:length), 'non-data2 = ') if( ptr.gt.0 ) then missing2 = cline((ptr+11):length) call bsquez(missing2) if( missing2(1:2).eq.'-1' ) missing2 = ' ' endif c ------ check for "scale" ptr = index( cline(1:length), 'scale = ') if( ptr.gt.0 ) then string = cline((ptr+7):length) call bsquez(string) status = mcstrtodbl( string, cal_scale) if( status.lt.0 ) cal_scale = 1.0D0 endif c ------ check for "offset" ptr = index( cline(1:length), 'offset = ') if( ptr.gt.0 ) then string = cline((ptr+8):length) call bsquez(string) status = mcstrtodbl( string, cal_offset) if( status.lt.0 ) cal_offset = 0.0D0 endif goto 1 endif list_ENVI = 0 return end subroutine initnavcom implicit none c --- navigation common character*4 NAV_proj integer NAV_lin integer NAV_ele double precision NAV_lat double precision NAV_lon double precision NAV_deglin double precision NAV_degele common /NAVCOM/ & NAV_proj, & NAV_lin, & NAV_ele, & NAV_lat, & NAV_lon, & NAV_deglin, & NAV_degele NAV_proj = ' ' NAV_lin = 0 NAV_ele = 0 NAV_lat = 0.0 NAV_lon = 0.0 NAV_deglin = 0.0 NAV_degele = 0.0 return end integer function read_lookup_cal( & file_name, & type, & sub_type, & band, & table, table2 & ) implicit none c --- parameters character*(*) file_name character*(*) type character*(*) sub_type integer band real table(1024) real table2(1024) c --- external functions character*12 cff integer lwfile integer mcstrtodbl integer getline c --- internal functions integer read_ASCII_cal integer read_KALPANA_cal c --- internal variables character*12 cvalue character*256 line character*256 file double precision dvalue integer i integer status integer length integer stat c --- test the file type if( sub_type(1:3).eq.'KAL' ) then call sdest('Kalpana read_cal',0) read_lookup_cal = read_KALPANA_cal( file_name, table, table2 ) else call sdest('Standard read_cal',0) read_lookup_cal = read_ASCII_cal( file_name, table ) endif return end ** Name: ** getline - get a line of text from the source ** ** Interface: ** integer function ** getline(character*(*) file, integer length, integer stat, ** character*(*) cline) ** ** Input: ** file - file name ** ** Input and Output: ** none ** ** Output: ** length - length of returned text line ** stat - line type (0=title, 1=text) ** cline - text string ** ** Return values: ** -1 - failure, invalid file type ** 0 - success, done with text read ** 1 - success, text line returned ** ** Remarks: ** none ** ** Categories: ** text integer function getline( infile, length, stat, cline ) implicit none include 'fileparm.inc' character*(MAXPATHLENGTH) path_file character*(*) infile character*(*) cline integer length integer index integer len_trim integer lwi integer mccmdstr integer mccmdquo integer volnam character*12 file_type character*160 text character*256 last_file character*256 file integer type integer ptr integer nword integer iline(20) integer init integer len integer len_file integer unit integer stat integer status integer title data init/0/ data unit/10/ data last_file/' '/ c --- status getline = -1 c --- initialization len_file = len_trim(infile) file = infile(1:len_file) if( last_file(1:len_file).ne.file(1:len_file) ) then last_file = file(1:len_file) call sdest(' NEW FILE= '//file(1:len_file),0) c ------ look for a title title = 0 if( mccmdquo( text ).ne.0 ) title = 1 c ------ pick up the file type c if( mccmdstr('FTYPE',1,'ASC',file_type).lt.0 ) return c call mcupcase( file_type ) file_type = 'ASC' call sdest(' FILE TYPE='//file_type,0) c ------ initialization for file type type = -1 if( file.ne.' ' ) then if ( file_type(1:1).eq.'A' ) then type = 1 status = volnam(file, path_file) OPEN (unit,FILE=path_file) elseif ( file_type(1:1).eq.'L' ) then type = 0 ptr = 0 nword = 20 else return endif else if( title.eq.0 ) return endif endif c --- if there is a title, read it first if( title.eq.1 ) then if( text.eq.' ' ) goto 100 len = index( text, '\\' )-1 if( len.gt.0 ) then cline = text(1:len) text = text(len+2:) else cline = text text = ' ' endif stat = 0 goto 500 endif c --- read a line of text from the file 100 if( type.eq.0 ) then if( lwi(file, ptr, nword, iline).ne.0 ) goto 1000 call movwc(iline,cline) ptr = ptr+nword stat = 1 elseif( type.eq.1 ) then READ (unit, 102, END=1000, ERR=1000, IOSTAT=status ) cline 102 FORMAT(A) stat = 1 else goto 1000 endif 500 continue call cleanc( cline ) c --- length of the line length = len_trim( cline ) getline = 1 return 1000 continue c --- reset the pointers if( type.eq.0 ) ptr = 0 if( type.eq.1 ) REWIND( unit ) if( title.eq.1 ) status = mccmdquo( text ) c --- done getline = 0 init = 0 return end C $ FUNCTION GRYSCL(TEMPK) (RCD) C $ CONVERT A BRIGHTNESS TEMPERATURE TO A GREY SCALE C $ TEMPK = (R) INPUT TEMPERATURE IN DEGRESS KELVIN C $$ GRYSCL = CONVERT INTEGER FUNCTION GRYSCL (TEMPK) IMPLICIT NONE REAL TEMPK ! initialized variables INTEGER CON1 INTEGER CON2 REAL TLIM DATA CON1 / 418 / DATA CON2 / 660 / DATA TLIM / 242.0 / C CONVERT A BRIGHTNESS TEMPERATURE TO A GREY SCALE C TEMPK---TEMPERATURE IN DEGREES KELVIN IF (TEMPK .GE. TLIM) GOTO 100 GRYSCL = MIN0 (CON1 - IFIX (TEMPK), 255) GOTO 200 100 CONTINUE GRYSCL = MAX0 (CON2 - IFIX (2 * TEMPK), 0) 200 CONTINUE RETURN END integer function read_ascii_lalo( & snav, & skip_lines, & nav_res, & nav_lines, & nav_elems, & lat_pos, & lat_len, & lon_pos, & lon_len, & rlat, & rlon & ) implicit none c --- parameters character*(*) snav integer skip_lines integer nav_res integer nav_lines integer nav_elems integer lat_pos integer lat_len integer lon_pos integer lon_len real rlat real rlon c --- external functions character*12 cfi integer getline integer mcstrtodll integer len_trim c --- internal variables character*80 cvalue character*12 clin character*12 cele character*12 clat character*12 clon character*256 sfile character*256 cline double precision dvalue integer line integer length integer stat integer status integer init integer slines integer this_line integer this_elem integer use_line integer use_elem integer blank integer len_cline data init/0/ c --- initialization if( init.eq.0 ) then c ------ local copy of the file name sfile = snav c ------ line counter for resolution reduction slines = 0 c ------ skip text lines if( skip_lines.gt.0 ) then do line = 1, skip_lines status = getline(sfile, length, stat, cline) call sdest('NAV LINE='//cline,0) if( status.lt.0 ) then read_ascii_lalo = -1 return elseif( status.eq.0 ) then read_ascii_lalo = -2 return endif enddo endif c ------ set initialization complete init = 1 endif c --- read a line of text 1 status = getline( sfile, length, stat, cline ) if( status.lt.0 ) then read_ascii_lalo = -3 return elseif( status.eq.0 ) then init = 0 read_ascii_lalo = 0 return else c ------ increment the line counter slines = slines+1 c ------ compute 2D coordinate from 1D array this_line = (slines/nav_elems) +1 this_elem = slines - ((this_line-1)*nav_elems) c ------ use Resolution factor to compute if we keep this lat/lon pair use_line = MOD( (this_line-1),nav_res ) use_elem = MOD( (this_elem-1),nav_res ) if( use_line.ne.0 .or. use_elem.ne.0 ) goto 1 clin = cfi( this_line ) cele = cfi( this_elem ) c ------ compress out unneeded blanks call bsquez( cline ) len_cline = len_trim( cline ) c ------ find first blank blank = index( cline(1:len_cline), ' ' ) c ------ read the lat parameter cvalue = cline(1:blank-1) status = mcstrtodll( cvalue, dvalue ) if( status.lt.0 ) then call sdest('FAILED: LAT='//cvalue,0) read_ascii_lalo = -4 return endif rlat = REAL( dvalue ) clat = cvalue c ------ read the lon parameter cvalue = cline(blank+1:len_cline) status = mcstrtodll( cvalue, dvalue ) if( status.lt.0 ) then call sdest('FAILED: LON='//cvalue,0) read_ascii_lalo = -5 return endif rlon = -REAL( dvalue ) clon = cvalue c call sdest('Lat='//clat//' Lon='//clon,0) read_ascii_lalo = 1 return endif end integer function read_KALPANA_cal( infile, temp, rad ) implicit none c --- parameters character*(*) infile real rad(*) real temp(*) c --- external functions character*12 cfi, cff integer lwi integer mccmdstr integer mcargparse integer mcargnum integer mcargint integer mcargdbl integer mcargfree integer getline c --- internal variables character*12 cindx character*12 crad character*12 ctemp character*80 file character*256 line integer handle integer i, j integer brit integer status integer length integer stat integer indx integer len double precision dvalue c --- check file name file = infile call sdest('Cal file = '//file,0) if( file.eq.' ' ) then read_KALPANA_cal = -1 return endif c --- read a line of text 1 status = getline( file, length, stat, line ) if( status.lt.0 ) then read_KALPANA_cal = -2 return elseif( status.eq.0 ) then read_KALPANA_cal = 0 return else c ------ parse the line of text handle = mcargparse( line, 0, len ) if( handle.lt.0 ) then read_KALPANA_cal = -3 return endif c ------ first value is the index status = mcargint(handle, ' ', 0, -1, 0, 1023, indx, cindx) if( status.lt.0 ) then read_KALPANA_cal = -4 return endif c ------ second value is the radiance status = mcargdbl(handle, ' ', 1, -999.0D0, & 0.0D0, 100.0D0, dvalue, crad) if( status.lt.0 ) then rad(indx+1) = -999.0 else rad(indx+1) = REAL(dvalue) endif c ------ third value is the temperature (kelvin) status = mcargdbl(handle, ' ', 2, -999.0D0, & 100.0D0, 400.0D0, dvalue, ctemp) if( status.lt.0 ) then temp(indx+1) = -999.0 else temp(indx+1) = REAL(dvalue) endif c ------ free the parsed line buffers status = mcargfree( handle ) c ------ go get another line goto 1 endif return end integer function read_ASCII_cal( & file_name, & table & ) implicit none c --- parameters character*(*) file_name real table(*) c --- external functions character*12 cff integer lwfile integer mcstrtodbl integer getline c --- internal variables character*12 cvalue character*256 line character*256 file double precision dvalue integer i integer status integer length integer stat c --- internal copy of file name file = file_name if( lwfile(file).eq.0 ) then read_ASCII_cal = -1 return endif c --- read the ASCII text i = 0 1 status = getline( file, length, stat, line ) if( status.lt.0 ) then read_ASCII_cal = -2 return elseif( status.eq.0 ) then read_ASCII_cal = 0 return else status = mcstrtodbl( line(1:length), dvalue ) if( status.ge.0 ) then i = i+1 table(i) = real( dvalue ) cvalue = cff( dvalue, 4) endif goto 1 endif return end integer function dec_RUN( & sfile, & tfile, & lines, & elems & ) implicit none c --- constants integer MSIZE parameter (MSIZE = 200000) integer MAXSITE parameter (MAXSITE = 1250) c --- parameters character*256 sfile ! source file name character*256 tfile ! temporary file name integer lines integer elems c --- external functions character*12 cfi integer lbi, lbo integer len_trim integer lwfile integer mccmdstr integer mccmdint integer mcsectodaytime integer mcskeyin c --- internal variables character*12 clines character*12 celems character*80 ddataset ! destination ADDE dataset group/descriptor.pos character*512 cmd ! command line buffer integer next integer run integer byte integer start integer i, j, k ! loop indicies integer site_array(MAXSITE) ! holds site layer data integer status ! function status integer year ! year of image integer month ! month of image integer day ! day of image integer time ! nominal hour of the image integer len_sfile ! length of the source file name integer len_tfile ! length of the temporary file name integer len_ddataset ! length of the destination dataset integer len_cmd ! length of the command string integer jstat(50) integer maxlvl integer i_value integer row, col ! product row and column coordinates integer hh, mm, ss ! Product hour, minute and second integer ibuf(MSIZE) integer lvl(1792,2560) integer rmax(1792,2560) integer*2 s_value(16) c --- MOSIAC common integer product_id ! product id number integer product_day ! product julian date (ccyyddd) integer product_hms ! product time (hhmmss) integer product_totalLength ! product length (in bytes) integer product_levels(16) ! product vip levels integer product_maxlvl ! product maximum vip level integer product_nRows ! product number of rows integer product_nColumns ! product number of columns integer product_blockLength ! product length of block integer product_nLayers ! product number of layers integer product_dataLength ! product length of data layer integer product_siteLength ! product length of site layer common /MOSIAC/ & product_id, & product_day, & product_hms, & product_totalLength, & product_levels, & product_maxlvl, & product_nRows, & product_nColumns, & product_blockLength, & product_nLayers, & product_dataLength, & product_siteLength c --- set function status to FAILED dec_RUN = -1 c --- initialize output rmax array do row = 1, 1792 do col = 1, 2560 rmax(row,col) = 0 enddo enddo maxlvl = 0 c --- does source file exist? len_sfile = len_trim( sfile ) status = lwfile(sfile) if( status.eq.0 ) then call edest('FAILED - '//sfile(1:len_sfile)// & ' does not exist',0) return endif c --- create a temporary file name tfile = sfile(1:len_sfile)//'.tmp' len_tfile = len_trim( tfile ) c --- read the image year status = mccmdint(' ',3,'Year',0,0,99,year) if( status.lt.0 ) return c --- read the image month status = mccmdint(' ',4,'Month',1,1,12,month) if( status.lt.0 ) return c --- read the image day status = mccmdint(' ',5,'Day',1,1,31,day) if( status.lt.0 ) return c --- read the image hour status = mccmdint(' ',6,'Hour',0,0,23,time) if( status.lt.0 ) return c --- product id status = lbi( sfile, 0, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, product_id, 0 ) call ddest(' Product Id = ', product_id) c -- days since 1/1/70 status = lbi( sfile, 2, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, i_value, 0 ) call ddest(' days since 1970001 = ',i_value) i_value = i_value * 86400 status = mcsectodaytime( i_value, product_day, product_hms ) call ddest(' Product Date = ',product_day) c --- seconds since midnight status = lbi( sfile, 4, 4, i_value ) call swbyt4( i_value, 1) hh = i_value / 3600 i_value = i_value - (hh*3600) mm = i_value / 60 ss = i_value - (mm*60) product_hms = (hh*10000) + (mm*100) + ss call ddest(' Product Time = ',product_hms) c --- total length status = lbi( sfile, 8, 4, i_value ) call swbyt4( i_value, 1) product_totalLength = i_value call ddest(' Product Total Length = ',product_totalLength) c --- data levels call ddest(' Product Levels:',0) status = lbi( sfile, 60, 32, s_value ) call swbyt2( s_value, 16) do i = 1, 16 call movb(2, s_value(i), product_levels(i), 0 ) if( i.eq.1 .and. product_levels(i).eq.32770 ) product_levels(i)=0 call ddest(' Level = ',product_levels(i)) enddo c --- maximum data level status = lbi( sfile, 92, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, product_maxlvl, 0 ) if( product_maxlvl.gt.maxlvl ) maxlvl = product_maxlvl call ddest(' Product Maximum Level = ',product_maxlvl) c --- number of rows status = lbi( sfile, 154, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, product_nRows, 0 ) call ddest(' Product Number of Rows = ',product_nRows) c --- number of columns status = lbi( sfile, 104, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, product_nColumns, 0 ) call ddest(' Product Number of Columns = ',product_nColumns) c --- length of block status = lbi( sfile, 124, 4, product_blockLength ) call swbyt4( product_blockLength, 4) call ddest(' Product Block Length = ', product_blockLength) c --- number of layers (2 means site ident layer is present) status = lbi( sfile, 128, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, product_nLayers, 0 ) call ddest(' Product Number of Layers = ',product_nLayers) c --- length of data layer status = lbi( sfile, 132, 4, product_dataLength ) call swbyt4( product_dataLength, 4) call ddest(' Product Length Data Layer = ',product_dataLength) c --- length of site layer product_siteLength = product_blockLength - product_dataLength - 16 call ddest(' Product Length Site Layer = ',product_siteLength) c --- fill the site layer i_value = product_totalLength - product_sitelength status = lbi( sfile, i_value, product_siteLength, site_array ) call swbyt4( site_array, product_siteLength/4 ) c --- read in the run-length encoded data start = 158 do row = 1, product_nRows status = lbi( sfile, start, 2, s_value ) call swbyt2( s_value, 1) call movb(2, s_value, next, 0 ) i=0 status = lbi( sfile, (start+2), next, ibuf ) call mpixel( next, 1, 4, ibuf ) do k=1,next byte = ibuf(k) run = byte/16 col = mod( byte,16 ) do j = 1, run i = i+1 lvl(row,i) = col if( lvl(row,i).gt.rmax(row,i) ) then rmax(row,i) = lvl(row,i) endif enddo enddo start = start + next + 2 enddo c --- write the destination image 1 row at a time byte = 0 do row = 1, product_nRows c --- make 1 line of data do col = 1, product_nColumns ibuf(col) = rmax(row,col) enddo c --- pack the array call mpixel( product_nColumns, 4, 1, ibuf ) c --- write the line to the destination file status = lbo( tfile, byte, product_nColumns, ibuf ) c --- increment the word pointer byte = byte + product_nColumns enddo c --- send back image dimensions lines = product_nRows elems = product_nColumns dec_RUN = 0 return end integer function list_VMC( & infile, & numlin, & numele, & lblbyt, & imgjday, & imgtime, & imgorbit, & imgnumber, & imgchannel, & imgmstime & ) implicit none c --- parameters character*(*) infile integer numlin integer numele integer lblbyt integer imgjday integer imgtime, imgmstime character*12 imgorbit character*12 imgnumber character*12 imgchannel c --- external functions character*12 cfj character*12 cfu integer len_trim integer lwi integer mccmdint integer mccmdstr integer mcargparse integer mcargstr integer mcargint integer mcargdbl integer mcargfree integer mcdmytocyd c --- internal variables character*2 cfilter character*3 CKWD_TYP (4) character*4 cmask character*4 csubdir character*4 csubnames (3) character*11 cprefix character*12 cnumlin character*12 cnumele character*12 clblbyt character*12 cjday character*12 ctime character*12 corbit character*12 cimage character*12 cversion character*12 csuffix character*30 PDS_KWD (255) character*30 CKWVAL (255) character*30 CDEFVAL (255) character*30 CKEYWORD character*30 CIMAGETIME character*30 CPERITIME character*40 cvalue character*80 cmes character*255 cfile character*255 cname character*255 cmd character*8000 cstring integer status integer IKWTYP (255) integer LENKWD (255) integer IDEFVAL (255) integer IMGYEAR integer IMGMONTH integer IMGDAY integer IMGH integer IMGM integer IMGS, IMGMS integer numkwd integer iorbit integer image_num integer image_start integer image_stop integer iversion integer handle integer len integer len_prefix integer len_subdir integer len_file integer len_suffix integer len_name integer len_path integer len_cmd double precision DEFVAL (255) double precision DEFMIN (255) double precision DEFMAX ( 255) double precision DKEYVAL (255) double precision ORBINC double precision ORBECC double precision ARGPER double precision PERIOD double precision ASNLON double precision PERALT double precision PERLAT double precision PERLON double precision SEMIMJ data DEFVAL / 255 * 0.0D0/ data DEFMIN / 255 * -1.0d32/ data DEFMAX / 255 * 1.0D32/ data DKEYVAL / 255 * 0.0D0 / data csubnames /'vmc1', 'vmc2', 'vmc3'/ c --- initialize function status to FAILED condition list_vmc = -1 c --- types of keyword formats CKWD_TYP (1) = 'STR' CKWD_TYP (2) = 'INT' CKWD_TYP (3) = 'SNG' CKWD_TYP (4) = 'DBL' c --- Keyword #1 PDS_KWD(1) = 'ASCENDING_NODE_LONGITUDE' IKWTYP(1) = 4 DEFVAL(1) = 90.0D0 LENKWD(1) = len_trim( PDS_KWD(1) ) c --- Keyword #2 PDS_KWD(2) = 'ORBITAL_INCLINATION' IKWTYP(2) = 4 DEFVAL(2) = 90.0D0 LENKWD(2) = len_trim( PDS_KWD(2) ) c --- Keyword #3 PDS_KWD(3) = 'ORBITAL_ECCENTRICITY' IKWTYP(3) = 4 DEFVAL(3) = 90.0D0 LENKWD(3) = len_trim( PDS_KWD(3) ) c --- Keyword #4 PDS_KWD(4) = 'PERIAPSIS_ARGUMENT_ANGLE' IKWTYP(4) = 4 DEFVAL(4) = 90.0D0 LENKWD(4) = len_trim( PDS_KWD(4) ) c --- Keyword #5 PDS_KWD(5) = 'PERIAPSIS_ALTITUDE' IKWTYP(5) = 4 DEFVAL(5) = 300.0D0 LENKWD(5) = len_trim( PDS_KWD(5) ) c --- Keyword #6 PDS_KWD(6) = 'PERIAPSIS_SEMIMAJOR_AXIS' IKWTYP(6) = 4 DEFVAL(6) = 39500.0D0 LENKWD(6) = len_trim( PDS_KWD(6) ) c --- ORBIT number parameter c NUMORB = IPP ( 3, IKWP('ORBIT',1,33)) c status = mccmdint(' ',1,'Orbit',0,1,9999,iorbit) c if( status.lt.0 ) return c call ddest('VMC Orbit Number = ', iorbit) c --- IMAGE keyword c ISTART = IPP ( 1, IKWP('IMAGE', 1, 0 )) c ISTOP = IPP ( 5, IKWP ('IMAGE', 2, IMAGE + 1 )) c status = mccmdint('IMA.GE',1,'Image',0,0,-1,image_start) c if( status.lt.0 ) return c status = mccmdint('IMA.GE',2,'Image',image_start,0,-1,image_stop) c if( status.lt.0 ) return c call ddest('Image Start Number = ',image_start) c call ddest('Image Stop Number = ',image_stop) c --- FILTER keyword c cfilter = CPP ( 2,CKWP('FILTER',1,'uv')) c status = mccmdstr('FIL.TER',1,'UV',cfilter) c if( status.lt.0 ) return c call mclocase( cfilter ) c call ddest('Filter Name = '//cfilter,0) c --- VERSION keyword c IVERSION = IPP ( 4, IKWP ('VERSION', 1, 2)) c status = mccmdint('VER.SION',1,'Version',2,0,-1,iversion) c if( status.lt.0 ) return c call ddest('Version Number = ',iversion) c --- use the version number to set the sub-directory name c csubdir = csubnames(IVERSION) c len_subdir = len_trim( csubdir ) c call ddest ('Sub-directory = ' //csubdir(1:len_subdir), 0) c --- use the orbit number to set the file name prefix c corbit = cfj( iorbit ) c cprefix =corbit(9:12)//'/v'//corbit(9:12)//'_' c len_prefix = len_trim( cprefix ) c call ddest('File Name Prefix = '//cprefix(1:len_prefix),0) c --- glue together the sub-directory and prefix c cfile = './'//csubdir(1:len_subdir)// c & '/'//cprefix(1:len_prefix) c len_file = len_trim( cfile ) c call ddest ('File Name = '//cfile(1:len_file), 0) c --- image number loop c do image_num = image_start, image_stop c ------ construct the file name suffix c cimage = cfj( image_num ) c cversion = cfj( iversion ) c csuffix = cimage(9:12)//'.'// c & cfilter(1:2)// c & cversion(12:12)//'.00' c len_suffix = len_trim( csuffix ) c call ddest('File Name Suffix = '//csuffix(1:len_suffix),0) c ------ construct the entire file name c cname = cfile(1:len_file)//csuffix(1:len_suffix) cname = infile len_name = len_trim( cname ) c ------ break the file name into orbit/number/channel components imgorbit = cname(2:5) imgnumber = cname(7:10) if( cname(12:13).eq.'n1' ) then imgchannel = 'IR1' elseif( cname(12:13).eq.'n2' ) then imgchannel = 'IR2' elseif( cname(12:13).eq.'uv' ) then imgchannel = 'UV' elseif( cname(12:13).eq.'vi' ) then imgchannel = 'VIS' else call edest('FAILED -- unrecognized file name',0) return endif call sdest( ' VMC orbit = '//imgorbit,0) call sdest( ' VMC number = '//imgnumber,0) call sdest( ' VMC channel = '//imgchannel,0) c ------ read the first 1024 words of the file status = lwi(cname(1:len_name), 0, 1024, cstring) if( status.ne.0 ) then call ddest('Cannot read '//cname, status ) goto 1 endif c ------ parse tokens handle = mcargparse ( cstring, 0, len ) c ------ LBLSIZE (size of the HEADER) status = mcargint( handle, 'LBLSIZE', 1, 80, 80, 2048, & LBLBYT, cvalue) call sdest(' VMC header = ', LBLBYT) clblbyt = cfu( lblbyt ) c ------ NL (number of image lines) status = mcargint( handle, 'NL', 1, 512, 512, 512, & NUMLIN, cvalue ) call sdest(' VMC lines = ', numlin) cnumlin = cfu( numlin ) c ------ NS (number of elements per line) status = mcargint( handle, 'NS', 1, 512, 512, 512, & NUMELE, cvalue ) call sdest(' VMC elems = ', numele) cnumele = cfu( numele ) c ------ ORBIT_NUMBER status = mcargint( handle, 'ORBIT_NUMBER', 1, 512, 512, 512, & iorbit, cvalue ) call ddest(' VMC ORBIT_NUMBER = ', iorbit) c ------ ORBITAL_ECCENTRICITY status = mcargdbl( handle, 'ORBITAL_ECCENTRICITY', 1, 0.0D0, & 0.0d0, 0.0d0, ORBECC, cvalue) write( cmes, 103) ORBECC 103 format(' VMC Orbital Eccentricity: ', G10.5) call ddest( cmes, 0 ) c ------ ORBITAL_INCLINATION status = mcargdbl( handle, 'ORBITAL_INCLINATION', 1, 0.0D0, & 0.0d0, 0.0d0, ORBINC, cvalue) write ( cmes, 104) ORBINC 104 format(' Orbital Inclinationy: ', G10.5) call ddest ( cmes, 0 ) c ------ ASCENDING_NODE_LONGITUDE STATUS = mcargdbl( handle, 'ASCENDING_NODE_LONGITUDE', 1, 0.0D0, & 0.0d0, 0.0d0, ASNLON, cvalue) write ( cmes, 105) ASNLON 105 format(' Longitude of Ascending Node: ', G10.5, 'Deg') call ddest ( cmes, 0 ) c ------ PERIAPSIS_SEMIMAJOR_AXIS status= mcargdbl( handle, 'PERIAPSIS_SEMIMAJOR_AXIS', 1, 0.0D0, & 0.0D0, 0.0D0, SEMIMJ, cvalue) write ( cmes, 106) SEMIMJ 106 format(' Semi-Major Axis: ', G10.5, 'Km') call ddest ( cmes, 0 ) c ------ PERIAPSIS_ARGUMENT_ANGLE status= mcargdbl( handle, 'PERIAPSIS_ARGUMENT_ANGLE', 1, 0.0D0, & 0.0D0, 0.0D0, ARGPER, cvalue) write ( cmes, 107) ARGPER 107 format(' Argument of Periapsis: ', G10.5, 'Deg') call ddest ( cmes, 0 ) c ------ PERIAPSIS_ARGUMENT_ANGLE status= mcargdbl( handle, 'PERIAPSIS_ARGUMENT_ANGLE', 1, 0.0D0, & 0.0D0, 0.0D0, ARGPER, cvalue) write ( cmes, 108) SEMIMJ 108 format(' Argument of Periapsis: ', G10.5, 'Deg') call ddest ( cmes, 0 ) c ------ PERIAPSIS_ALTITUDE STATUS = mcargdbl( handle, 'PERIAPSIS_ALTITUDE', 1, 0.0D0, & 0.0D0, 0.0D0, PERALT, cvalue) write ( cmes, 109) PERALT 109 format(' Periapsis Altitude: ', G10.5, 'Km') call ddest ( cmes, 0 ) c ------ PERIAPSIS_TIME STATUS = mcargstr( handle, 'PERIAPSIS_TIME', 1, ' ', CPERITIME) call ddest (' PERIAPSIS_TIME: '//CPERITIME, 0) c ------ IMAGE_TIME STATUS = mcargstr( handle, 'IMAGE_TIME', 1, ' ', CIMAGETIME) c call sdest (' IMAGE_TIME: '//CIMAGETIME, 0) read (CIMAGETIME,110) IMGYEAR, & IMGMONTH, & IMGDAY, & IMGH, & IMGM, & IMGS, & IMGMS 110 format(i4, 1x, i2, 1x, i2, 1x, i2, 1x, i2, 1x, i2, 1x, i3) call ddest (' Year: ', IMGYEAR) call ddest (' HH: ', IMGH) call ddest (' MM: ', IMGM) call ddest (' SS: ', IMGS) call ddest (' MS: ', IMGMS) imgtime = IMGH*10000 + IMGM*100 + IMGS imgmstime = IMGH*10000000 + IMGM*100000 + IMGS*1000 + IMGMS write (cmes, 111) IMGDAY, IMGMONTH, IMGYEAR, IMGTIME, IMGMS 111 format ('Image acquired on ', i3, i3, i6, ' at ', i8, '.', i3) c call ddest (cmes, 0) status=mcdmytocyd( imgday,imgmonth,imgyear,imgjday ) cjday = cfu( imgjday ) cvalue = cfj( imgtime ) ctime = cvalue(7:8)//':'//cvalue(9:10)//':'//cvalue(11:12) c ------ IMGMAKE command c cmd = 'IMGMAKE'// c & ' '//cname(1:len_name)// c & ' A/A.1'//cimage(10:12)// c & ' '//cnumlin// c & ' '//cnumele// c & ' HEADER='//clblbyt// c & ' NBPE=2'// c & ' DAY='//cjday// c & ' TIME='//ctime// c & ' NBYT=2'// c & ' PLANET=VENUS' c call bsquez( cmd ) c len_cmd = len_trim( cmd ) c call sdest( cmd(1:len_cmd), 0 ) c ------ free the parse handle status = mcargfree( handle ) 1 continue c enddo ! End of image numebr loop c --- set function status to SUCCESS list_vmc = 0 return end integer function vmcread( & navfile, & imgorbit, & imgnumber, & imgchannel & ) implicit none C VMCREAD | Reads VMC navigation text file for matching orbit/image/channel C VMCREAD (navfile) (orbit) (image) (channel) C Parameters: C navfile | name of the navigation text file (def=none) C orbit | orbit number to search for (def=none) C image | image number to search for (def=none) C channel | image channel to search for C ALL = match on all chennels matching orbit and image (def) C VIS = match VIS entry only C UV = match UV entry only C IR1 = match NIR-1 entry only C IR2 = match NIR-2 entry only C Keywords: C none C Remarks: C If there is no entry in the file for the specified image number, the C next entry for the orbit/channel will be listed. C c --- parameters character*256 navfile character*12 imgorbit character*12 imgnumber character*12 imgchannel c --- external functions character*12 cfi character*12 cfu character*12 cff integer lbi integer lwd integer lwfile integer lwFileLength integer len_trim integer mcargparse integer mcargstr integer mcargint integer mcargdbl integer mcargdll integer mcargfree integer mcgetdaytime integer mcstrtoint integer read_record c --- internal variables character*12 c_value character*12 c_lat, c_lon character*12 c_alt character*12 c_cemay, c_cemax character*12 c_polean character*12 c_susolat, c_susolon character*12 c_area character*12 d_value character*12 find_channel character*12 file_date character*12 file_time character*80 file character*80 file_datetime character*80 file_channel character*80 file_sclktime character*1000 line character*512 cmd integer len integer len_channel integer len_ppcmd integer ptr integer handle integer len_line integer num_parm integer nrec integer nbytes integer test integer nchar integer status integer read_status integer i integer len_record integer find_orbit integer find_image integer test_image integer file_orbit integer file_image integer file_hit double precision file_lon double precision file_lat double precision file_dven double precision file_centan double precision file_polean double precision file_susoan double precision file_susolon double precision file_susolat double precision file_cemax double precision file_cemay double precision file_polex double precision file_poley double precision file_susox double precision file_susoy logical FOUND c --- post process command integer ipos character*160 ppcmd common /PPCOM/ & ipos, & ppcmd c --- set application status vmcread = -1 FOUND = .FALSE. c --- name of the VMC navigation file file = navfile c --- find this orbit number status = mcstrtoint( imgorbit, find_orbit ) if( status.lt.0 ) return c --- find this image number status = mcstrtoint( imgnumber, find_image ) if( status.lt.0 ) return c --- find this image channel find_channel = imgchannel if( & find_channel.eq.'VIS' .or. & find_channel.eq.'UV' .or. & find_channel.eq.'IR1' .or. & find_channel.eq.'IR2' .or. & find_channel.eq.'ALL' & ) then if( find_channel.eq.'IR1' ) then find_channel = 'NIR-1' elseif( find_channel.eq.'IR2' ) then find_channel = 'NIR-2' endif len_channel = len_trim( find_channel ) else call edest('Invalid Channel specified = '//find_channel,0) return endif c --- length of the VMC nav file nbytes = lwFileLength( file )*4 50 status = lbi( file,nbytes,1,test) if( status.eq.0 ) then nbytes = nbytes+1 goto 50 else nbytes = nbytes-1 endif c --- read a line of text from the navigation file 100 continue line = ' ' read_status = read_record ( & file, & nbytes, & nrec, & nchar, & line & ) if( read_status.lt.0 ) then call edest('FAILED - reading text lines',0) return endif c --- parse the line handle = mcargparse( line, 0, len_line ) c --- check for valid integer status = mcargstr(handle,' ',0,' ',c_value) if( status.lt.0 ) return status = mcstrtoint( c_value, file_orbit) if( status.ge.0 ) then c ------ image number status = mcargint(handle,' ',1,-1,0,999,file_image,c_value) if( status.lt.0 ) return c ------ check for orbit and image number for match if( & file_orbit.eq.find_orbit .and. & file_image.ge.find_image & ) then c --------- image channel status = mcargstr(handle,' ',3,' ',file_channel) if( status.lt.0 ) return c --------- check for channel match if( find_channel.ne.'ALL' ) then status=index( file_channel, find_channel(1:len_channel) ) if( status.le.0 ) goto 666 endif c --------- set FOUND flag if( .NOT.FOUND ) then FOUND = .TRUE. test_image = find_image endif c --------- check the test_image number for a match if( file_image.ne.test_image ) goto 666 call sdest('VMC Planetary Navigation Elements:',0) call sdest(' Orbit= ',file_orbit) call sdest(' Image= ',file_image) call sdest(' Channel='//file_channel,0) c --------- image date/time status = mcargstr(handle,' ',2,' ',file_datetime) if( status.lt.0 ) return len = len_trim( file_datetime ) ptr = index( file_datetime, 'T' ) if( ptr.gt.0 ) then file_date = file_datetime(1:ptr-1) file_time = file_datetime(ptr+1:len) call sdest(' Date= '//file_date,0) call sdest(' Time= '//file_time,0) endif c --------- image lon status = mcargdll(handle,' ',4,999.0D0,-180.0D0,180.0D0, & file_lon,c_lon) if( status.lt.0 ) return len = len_trim( c_lon ) call sdest(' Lon='//c_lon(1:len),0) c --------- image lat status = mcargdll(handle,' ',5,999.0D0,-90.0D0,90.0D0, & file_lat,c_lat) if( status.lt.0 ) return len = len_trim( c_lat ) call sdest(' Lat='//c_lat(1:len),0) c --------- image d_ven status = mcargdbl(handle,' ',6,9999.0D0,0.0D0,-1.0D0, & file_dven,c_alt) if( status.lt.0 ) return len = len_trim( c_alt ) call sdest(' D_ven='//c_alt(1:len),0) c --------- image cent_an status = mcargdbl(handle,' ',7,9999.0D0,360.0D0,-1.0D0, & file_centan,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Cent_an='//c_value(1:len),0) c --------- image pole_an status = mcargdbl(handle,' ',8,9999.0D0,-360.0D0,360.0D0, & file_polean,c_polean) if( status.lt.0 ) return len = len_trim( c_polean ) call sdest(' Pole_an='//c_polean(1:len),0) file_polean = file_polean + 90.0D0 if( file_polean.gt.360.0D0 ) then file_polean = 360.0D0 - file_polean endif c_polean = cff( file_polean, 4 ) c --------- image suso_an status = mcargdbl(handle,' ',9,9999.0D0,0.-360D0,360.0D0, & file_susoan,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Suso_an='//c_value(1:len),0) c --------- image suso lon status = mcargdll(handle,' ',10,999.0D0,-180.0D0,180.0D0, & file_susolon,c_susolon) if( status.lt.0 ) return len = len_trim( c_susolon ) call sdest(' Suso Lon='//c_susolon(1:len),0) c --------- image suso lat status = mcargdll(handle,' ',11,999.0D0,-90.0D0,90.0D0, & file_susolat,c_susolat) if( status.lt.0 ) return len = len_trim( c_susolat ) call sdest(' Suso Lat='//c_susolat(1:len),0) c --------- image cema_x status = mcargdbl(handle,' ',12,9999.0D0,0.0D0,-1.0D0, & file_cemax,c_cemax) if( status.lt.0 ) return len = len_trim( c_cemax ) call sdest(' Cema X='//c_cemax(1:len),0) c --------- image cema_y status = mcargdbl(handle,' ',13,9999.0D0,0.0D0,-1.0D0, & file_cemay,c_cemay) if( status.lt.0 ) return len = len_trim( c_cemay ) call sdest(' Cema Y='//c_cemay(1:len),0) c --------- image pole_x status = mcargdbl(handle,' ',14,9999.0D0,0.0D0,-1.0D0, & file_polex,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Pole X='//c_value(1:len),0) c --------- image pole_y status = mcargdbl(handle,' ',15,9999.0D0,0.0D0,-1.0D0, & file_poley,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Pole Y='//c_value(1:len),0) c --------- image suso_x status = mcargdbl(handle,' ',16,9999.0D0,0.0D0,-1.0D0, & file_susox,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Suso X='//c_value(1:len),0) c --------- image suso_y status = mcargdbl(handle,' ',17,9999.0D0,0.0D0,-1.0D0, & file_susoy,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Suso Y='//c_value(1:len),0) c --------- image hit status = mcargint(handle,' ',18,0,0,-1, & file_hit,c_value) if( status.lt.0 ) return len = len_trim( c_value ) call sdest(' Hit='//c_value(1:len),0) c --------- image sclk_time status = mcargstr(handle,' ',19,' ',file_sclktime) if( status.lt.0 ) return len = len_trim( file_sclktime ) call sdest(' Sclk Time='//file_sclktime(1:len),0) endif endif c --- free the handle 666 status = mcargfree( handle ) c --- read another record if( read_status.eq.0 ) goto 100 c --- Error if no match was found if( .NOT.FOUND ) then call sdest('No Match was found',0) return endif c --- make post process VMCNAV string c_area = cfu( ipos ) cmd = & ' VMCNAV'// & ' AREA='//c_area(1:4)// & ' PLANET=2999'// & ' VIEWPOINT='//c_lat//' '//c_lon//' '//c_alt//' '// & c_cemay//' '//c_cemax// & ' ANGLE='//c_polean// & ' SUN='//c_susolat//' '//c_susolon// & ' REPLACE=NAV' call bsquez( cmd ) len = len_trim( cmd ) if( len.gt.160 ) len=160 ppcmd = cmd(1:len) c --- write the command into the string table call lbput( 'VMCNAV', ppcmd(1:len) ) c --- print the Naviagtion command line call sdest(' SAMPLE VMCNAV command:',0) i = index( ppcmd, 'ANGLE' ) call sdest(' '//ppcmd(1:i-1),0) call sdest(' '//ppcmd(i:len),0) c --- set application status 9999 vmcread = 0 return end integer function readENVInav( * lat_file, * lon_file, * nlines, * nelems, * rlats, * rlons * ) implicit none c --- parameters character*(*) lat_file character*(*) lon_file integer nlines integer nelems real rlats(*) real rlons(*) c --- external functions character*12 cff integer lbi integer read_record integer lwFileLength integer len_trim integer mcargparse integer mcargdbl integer mcargstr integer mcargfree integer mcstrtodbl c --- internal variables character*12 cvalue character*256 file character*5000 line double precision dvalue integer pass integer i, k integer status integer read_status integer len_file integer len_line integer nbytes integer handle integer nline integer nchar integer test c --- set application status to FAILED readENVInav = -1 c --- file name do pass = 1, 2 if( pass.eq.1 ) then file = lat_file else file = lon_file endif len_file = len_trim( file ) c --- length of the temporary file nbytes = lwFileLength( file(1:len_file) )*4 50 status = lbi( file,nbytes,1,test) if( status.eq.0 ) then nbytes = nbytes+1 goto 50 else nbytes = nbytes-1 endif call sdest(' readENVI: FileLength=',nbytes) c --- read a line of text from the temporary file k = 0 1 continue line = ' ' read_status = read_record( file, nbytes, nline, nchar, line ) if( read_status.lt.0 ) then call edest('FAILED -- Reading text file',0) return else c ------ length of line string call bsquez( line ) len_line = len_trim( line ) c ------ parse the line handle = mcargparse( line, 0, len_line ) c ------ check for numbers status = mcargstr(handle,' ',1,' ',cvalue) status = mcstrtodbl( cvalue, dvalue ) if( status.lt.0 ) goto 9999 if( k.eq.0 ) then call sdest(' readENVI: line='//line(1:40),0) endif c ------ read the lat/lon arrays do i = 1,nelems k = k+1 status = mcargstr(handle,' ',i-1,' ',cvalue) if( cvalue(1:3).eq.'NaN' ) then dvalue = -999.0D0 else status = mcargdbl(handle,' ',i-1,-999.0D0, & -360.0D0,360.0D0,dvalue,cvalue) endif if( k.eq.1 ) then call sdest(' readENVI: cvalue(1)='//cvalue,0) endif if( pass.eq.1 ) then rlats(k) = -999.0 if( dvalue.le.90.0D0.and.dvalue.ge.-90.0D0 ) then rlats(k) = real(dvalue) endif if( k.eq.1 ) then cvalue = cff(dble(rlats(k)),4) call sdest(' readENVI: rlat='//cvalue,0) endif else rlons(k) = -999.0 if( dvalue.le.360.0D0.and.dvalue.ge.0.0D0 ) then if( dvalue.gt.180.0D0 ) dvalue=dvalue-360.0D0 rlons(k) = real(dvalue) endif if( k.eq.1 ) then cvalue = cff(dble(rlons(k)),4) call sdest(' readENVI: rlon='//cvalue,0) endif endif enddo c ------ free the handle 9999 status = mcargfree( handle ) endif c --- read another line if( read_status.eq.0 ) goto 1 call sdest(' readENVI: pass=',pass) call sdest(' readENVI: values read=',k) enddo readENVInav = 0 return end integer function read_record( file, nbytes, nline, nchar, line ) implicit none c --- parameters character*80 file integer nbytes integer nline integer nchar character*5000 line c --- external functions integer lbi integer len_trim c --- internal variables character*5000 ctext character*80 this_file character*80 last_file integer lf integer text(1250) integer ptr integer byte integer rbyte, ibyte integer init integer status integer bytes_left integer nread equivalence (text, ctext) data init/0/ data LF/Z'0A'/ data last_file/' '/ this_file = file if( last_file.ne.this_file ) then ibyte = 0 nline = 0 init = 1 last_file = this_file endif c --- read the next record rbyte = 5000 if( (ibyte+rbyte).gt.nbytes ) then rbyte = nbytes - ibyte endif ctext= ' ' status =lbi( file, ibyte, rbyte, text) if( status.lt.0 ) then ptr = index( ctext, char(LF) ) if( ptr.gt.0 ) then nline = nline+1 nchar = ptr-1 line(1:nchar) = ctext(1:nchar) ibyte = ibyte + (nchar+1) endif if( ibyte.ge.nbytes ) then read_record = 1 init = 0 else nline = nline+1 nchar = len_trim( ctext ) line(1:nchar) = ctext(1:nchar) read_record = 0 endif else ptr = index( ctext, char(LF) ) if( ptr.gt.0 ) then nline = nline+1 nchar = ptr-1 line(1:nchar) = ctext(1:nchar) ibyte = ibyte + (nchar+1) read_record = 0 else nline = nline+1 nchar = len_trim( ctext ) line(1:nchar) = ctext(1:nchar) read_record = 1 endif endif return end INTEGER FUNCTION lwFileLength(NAME) IMPLICIT NONE CHARACTER*(*) NAME C --- external functions INTEGER M0CACHEOPEN ! returns file descriptor INTEGER M0FILESIZE ! returns size of file integer len_trim C --- local variables character*256 file integer len_file integer l ! size of the file in bytes file = name len_file = len_trim( file ) call sdest(' lwFileLength: file='//file,0) L = M0FILESIZE(M0CACHEOPEN(file(1:len_file))) call sdest(' lwFileLength: file size=',l) IF (L.LT.0) THEN lwFileLength = -1 ELSE IF (L.EQ.0) THEN lwFileLength = 0 ELSE lwFileLength = (L-1)/4 END IF RETURN END