CTITLE 'CASEFOLD' subroutine casefold( string ) * Routine to convert a string to upper case. * i - Loop counter * len_string - used length of string * trimlen - returns string length integer*4 i, len_string, trimlen * string - Used portion of string character*(*) string ***** Get length of string len_string = trimlen(string) * Now loop and convert do i = 1, len_string if( string(i:i).ge.'a' .and. string(i:i).le.'z') then string(i:i) = CHAR( ICHAR(string(i:i))-32 ) end if end do **** Thats all return end CTITLE 'PROPER_RUNSTRING' subroutine proper_runstring( help_file, program, kill) * Routine to ouput the runstring for nutcorr * To use this subroutine the user needs to set the enviroment * variable HELP_DIR. This directory is then preprended to * the file name * ierr - IOSTAT error (opening help file) * trimlen - Length of string * kill - indicates iof program should be killed. * 0 means do not stop program * >0 kill program * <0 kill program but list help in 24 line * segment (with q to quit) * len_dir - Length of directory name * cnt - Keeps track of number of lines written * so that user can be quizzed to continue if * kill set < 0 integer*4 ierr, trimlen, kill, len_dir, cnt * line - Line read from input. character*100 line * help_file - Name of help file * program - Name of program. character*(*) help_file, program * Full_help_file - Help file name with directory added. character*256 Full_help_file * ans - Answer to query for more lines character*4 ans ***** open the help file. First construct name call getenv('HELP_DIR',Full_help_file) len_dir = trimlen(Full_help_file) if( len_dir.gt.0 ) then Full_help_file(len_dir+1:) = help_file else Full_help_file = help_file end if open(999, file=full_help_file, iostat=ierr, status='old') * NOD TAH 980813: If file not found tell user that to do if( ierr.ne.0 ) then write(*,110) 110 format(/,'**ERROR** Opening help file. Check that', . ' enviroment variable HELP_DIR points to',/, . ' directory with help files.') end if call report_error('IOSTAT',ierr,'open',full_help_file, 0, . program) * ! Loop until EOF cnt = 0 do while ( ierr.eq.0 ) read(999,'(a)', iostat=ierr) line cnt = cnt + 1 if( ierr.eq.0 ) write(*,'(a)') line(1:max(trimlen(line),1)) if( mod(cnt,24).eq.0 .and. kill.lt.0 ) then * see if we should continue list write(*,120) 120 format('More? (q to quit) ',$) read(*,'(a)') ans if( ans(1:1).eq.'q' .or. ans(1:1).eq.'Q' ) then ierr = -1 end if end if end do ***** close the input close(999) **** Now see if we should stop if( kill.ne.0 ) then write(*,200) program(1:max(trimlen(program),1)) 200 format(' PROGRAM ',a,' terminating') stop ' Incomplete runstring' end if **** Thats all return end CTITLE 'REPORT_ERROR' subroutine report_error(type,ierr,action,file, kill, prog) * Routine to Report that an error has occurred and possibly * stop the program (depending on the value of Kill: no kill if * 0, otherwize kill) * ierr - IOSTAT error * kill - Kill program if value is not zero. * len_type - Used length of type * len_action - Used length of action * len_file - Used length of file * len_prog - Used length of prog * trimlen - Get length of string integer*4 ierr, kill, len_type, len_action, len_file, len_prog, . trimlen * type - Type of error * action - Thing being done when error occurred * file - Name of file being acted on * prog - Program/subroutine name of calling * - routine. character*(*) type, action, file, prog * Additional variables needed for use of report_stat with this * routine * rcpar - Get the name of the program calling this routine * len_mod - Length of the module name * jerr - IOSTAT error creating message line * message - The concatinated error message * module - name of program integer*4 rcpar , len_mod, jerr character*128 module character*128 message **** If error is not zero, print message if( ierr.ne.0 ) then call check_ascii( type ) call check_ascii( action ) call check_ascii( file ) call check_ascii( prog ) len_type = max(1,trimlen(type)) len_action = max(1,trimlen(action)) len_file = max(1,trimlen(file)) len_prog = max(1,trimlen(prog)) if( len_file+len_action+len_type+len_prog.lt.55 ) then write(*,100) type(1:len_type), ierr, . action(1:len_action), . file(1:len_file), prog(1:len_prog) 100 format(1x,a,' error ',i6,' occurred ',a,'ing ',a, . ' in ',a) else write(*,120) type(1:len_type), ierr, . action(1:len_action), . file(1:len_file), prog(1:len_prog) 120 format(1x,a,' error ',i6,' occurred ',a,'ing',/, . a,/,' in ',a) end if **** Set up the report_stat string len_mod = rcpar(0, module) if( len_mod.le.0 ) then module = 'MAIN' len_mod = 4 end if write(message,140, iostat=jerr) type(1:len_type), . action(1:len_action) 140 format(a,' error ',a,'ing file') if( kill.ne.0 ) then call report_stat('fatal',module, prog, file, . message,ierr) else call report_stat('warning',module, prog,file, . message,ierr) end if **** Now see if we should kill if( kill.ne.0 ) then write(*,200) prog(1:len_prog) 200 format(' Program terminating in routine ',a) stop ' Stopped from report_error' end if end if ***** Thats all return end CTITLE RCPAR integer*4 function rcpar( iel, arg ) * Routine to emulate RCPAR using the igetarg UNIX subroutine * Modified to use getarg * iel - Element of runstring to get * igetarg - UNIX routine to read runstring * len_arg - Length of arg. * trimlen - Get length of string integer*4 iel, len_arg, trimlen * arg - Arg of runstring character*(*) arg **** Get length of argument and runstring len_arg = LEN(arg) call getarg( iel, arg ) rcpar = trimlen( arg ) ***** Thats all return end CTITLE SUB_CHAR subroutine sub_char ( buffer, old, new ) * Routine to search string buffer and replace all occurrences * of string OLD with string NEW * end - End position of moving buffer * i,j - Loop counter * len_buffer - length of buffer (full dimensioned length) * len_old - Length of old string * len_new - Length of new string * len_used - Used portion of buffer * offset - Difference in length between new and old string * pos - Position of old in string * start - start position for moving buffer * trimlen - Length of used portion of string integer*4 end, i, len_buffer, len_old, len_new, len_used, . offset, pos, start, trimlen * finished - Indicates we have run out of string logical finished * buffer - Buffer to be manipulated * old - String to be replaced * new - string to replace old with character*(*) buffer, old, new ***** Start by get sizes of strings len_old = len( old ) len_new = len( new ) len_buffer = len( buffer ) * Check to make sure old is not a substring of new. If it is the * algorithm used here will not work pos = index ( new, old ) if( pos.gt.0 ) RETURN finished = .false. do while ( .not. finished ) pos = index( buffer, old ) * ! Make substiution if( pos.gt.0 ) then * Now copy string after pos len_used = trimlen( buffer ) offset = len_new - len_old start = pos - offset + 1 end = min(len_used,len_buffer+offset) * ! Move string starting at top if( len_new.gt.len_old ) then do i = end, start, -1 buffer(i+offset:i+offset) = buffer(i:i) end do end if * ! Move string starting a bottom if( len_new.lt.len_old ) then do i = start, end buffer(i+offset:i+offset) = buffer(i:i) end do * ! Clear the characters left at buffer(end+offset+1:) = ' ' * ! the end end if **** Now put in new string buffer(pos:pos+len_new-1) = new * ! we have finished ELSE finished = .true. ENDIF * ! Searching through string END DO **** Thats all return end CTITLE 'dvdot' subroutine dvdot(r, v1, inc1, v2,inc2, iter) c c Routine for a double precision dot product. Result returned c in r c * inc1,inc2 - increment on V1, and v2 * i - loop counter * i1,i2 - position in V1 and v2 * iter - number of values to be searched integer*4 inc1,inc2, i, i1,i2, iter * c c v1(inc1,1),v2(inc2,1) ! vectors to be dotted c r - result of dot product real*8 v1(1), v2(1), r * c c **** Loop swapping values i1 = 1 i2 = 1 r = 0.d0 do i = 1, iter r = r + v1(i1)*v2(i2) i1 = i1 + inc1 i2 = i2 + inc2 end do c ***** Thats all return END CTITLE 'TRIMLEN' integer*4 function trimlen(string) * Routine to return the length of the used portion of string. * Length is with trailing blanks removed. * len_string - declared length of string * i - Loop counter integer*4 len_string, i * blanks - Indicates blanks are still being * - found logical blanks * string - the string whose length we want character*(*) string ***** Get full length of string, and work backwards len_string = LEN(string) i = len_string blanks = .true. * Scan string from end to front do while ( blanks ) if( i.gt.0 ) then if( string(i:i).eq.' ' ) then i = i - 1 else blanks = .false. end if else * ! Force exit at end of string blanks = .false. end if end do ***** Save length of string trimlen = i ***** Thats all return end CTITLE 'YMDHMS_TO_JD' SUBROUTINE YMDHMS_to_JD ( date, seconds, epoch ) * ------------------------- * ., 26JUL85 <871210.1317> * Author: T.Herring 3:53 PM THU., 10 JULY, 1986 * *$S "YMDHMS_to_JD" *----------------------------------------------------------------------- * Routine to convert a calender date with hours, minutes and * seconds to a Modified Julian date. The calender date is * ordered as year, month, day of month, hours, and * minutes. These values are stored in a single I*2 array. * * If the year is greater than 1000 then the it is assumed to * contain the centuries. * * This routine is only valid for dates after 1600 Jan 0. * * CALLING SEQUENCE: * ================= * CALL ymdhms_to_JD ( date, seconds, epoch) * * WHERE: * date is a 5 element array containing year, month, day of month, * hour, and minutes. * (I*4 5 element array INPUT) * seconds is the seconds part of the epoch * (REAL*8 INPUT) * epoch is the JD with fractional days corresponding to date and * seconds. * (REAL*8 OUTPUT) * *---------------------------------------------------------------------- *$E * day - day of month * date(5) - 5 element date array with year, month, day, * - hours and minutes. * days_to_month(12) - number of days from start of year to * - each month * leap_days - the number of leap days we need to include * month - month of year * year - the years since 1900.0 * years_from_1600 - number of years since 1600 Jan 0. * days_from_1600 - number of days since 1600 Jan 0. integer*4 day, date(5), days_to_month(12), leap_days, month, . year, years_from_1600, days_from_1600 * epoch - The JD corresponding to date and seconds * fraction - Fraction of a day. We compute this separately * - to avoiod some rounding error when added to MJD * mjd - the computes MJD with fractional days * seconds - the seconds part of the epoch. real*8 epoch, fraction, mjd, seconds * leap_year - Indicates that this is a leap year logical leap_year data days_to_month / 0, 31, 59, 90, 120, 151, . 181, 212, 243, 273, 304, 334 / ***** START, Make sure year is from 1900 year = date(1) month = date(2) day = date(3) if( year.lt.1000 ) year = year + 1900 ***** Compute number of years from 1600 years_from_1600 = year - 1600 ***** Now compute number of leap days upto the start of the year * we are in (i.e., if the year we are in is a leap year do not * add the leap day yet) leap_days = (years_from_1600 - 1)/ 4 . - (years_from_1600 + 99)/100 . + (years_from_1600 + 399)/400 + 1 if( years_from_1600.eq.0 ) leap_days = leap_days - 1 ***** Now see if we are in leap year leap_year = .false. if( mod(years_from_1600, 4).eq.0 .and. . (mod(years_from_1600,100).ne.0.or. . mod(years_from_1600,400).eq.0) ) leap_year = .true. ***** Now compute number of days sinec 1600 days_from_1600 = years_from_1600*365 + leap_days + . days_to_month(month) + day ***** Add extra day if we are after Februrary and this is a leap year if( month.gt.2 .and. leap_year ) then days_from_1600 = days_from_1600 + 1 end if ***** Compute the mjd and add in the fraction of a day part * The MJD of 1600 Jan 0 is -94554 fraction = seconds/86400.d0 + date(5)/1440.d0 + date(4)/24.d0 mjd = -94554.d0 + days_from_1600 + fraction epoch = mjd + 2 400 000.5d0 ***** THATS ALL RETURN end CTITLE 'JD_TO_YMDHMS' SUBROUTINE JD_to_YMDHMS ( epoch, date, seconds ) * ------------------------ * ., 10JUL86 <871210.1317> * Author: T.Herring 11:28 AM THU., 10 JULY, 1986 * *$S "JD_to_YMDHMS" *---------------------------------------------------------------------- * Routine to convert a Julian date (with fractions of a day) * or a full Julian date, to a calender data with hours and minutes, * and floating-point seconds (Real*8). * * NOTE: If a full Julian date is used the resolution of the seconds * will only be about 10 microseconds, a MJD should yield a resolution * of about 0.1 microseconds in the seconds. * * This routine is only valid for dates after 1600 Jan 0 (MJD -94554) * * CALLING SEQUENCE: * ================= * CALL JD_to_YMDHMS ( epoch, date, seconds ) * * WHERE: * epoch is a modified or full julian date with possibly a fractional * part of a day. The two type of input are destinguished by * the full Julian date being greater than 2,000,000. * (REAL*8 INPUT) * date is an array containing the calender date with (full) year, * month of year, day of month, hours, minutes. Note the year * will be returned with the centuries added. * (I*4 5 element array OUTPUT) * seconds is the floating point seconds part of the MJD * (REAL*8 OUTPUT) * *---------------------------------------------------------------------- * *$E * century - Number of century from 1600 Jan 0. * date(5) - the calender date corresponding to * - 'epoch' resolution to the minute * day - day of month * day_of_year - Number of days from start of year * days_to_month(13) - number of days to start of each month * - in a non-leap-year * month - month of year * year - years since start of century * years_from_1600 - Number of years since 1600 Jan 0. * days_from_1600 - Number of days elapsed since 1600 Jan 0. * - (MJD -94554.0 Julian date 2305447.0) integer*4 century, date(5), day, day_of_year, days_to_month(13), . month, year, years_from_1600, days_from_1600 * epoch - the julian date or modified julian date * - to be converted to calender date * fraction - the fraction of a day part of MJD * mjd - epoch converted to a MJD * mjd_day - the whole number days in the mjd * seconds - the seconds part of the MJD (<60) real*8 epoch, fraction, mjd, mjd_day, seconds * leap_year - Indicates that this days is a leap year logical leap_year data days_to_month / 0, 31, 59, 90, 120, 151, 181, . 212, 243, 273, 304, 334, 365 / ***** START, convert full jd to mjd if we have too. * ! epoch must be Full julian date if( epoch.gt.2000000.d0 ) then * ! convert to MJD. mjd = epoch - 2 400 000.5d0 else mjd = epoch end if ***** Remove the fractional part of the mjd mjd_day = aint ( mjd ) fraction = mod ( mjd,1.d0 ) if( mjd.lt.0 .and. fraction.ne.0.d0 ) then mjd_day = mjd_day - 1 fraction = fraction + 1.d0 end if ***** Now convert MJD (even day) to date (year, month, day ) * Get number of days since 1600. days_from_1600 = mjd_day - (-94554.0d0) years_from_1600 = days_from_1600/365.d0 ***** Now compute day_of_year and years_from_1600 accounting for * leap years * ! Just to get us into loop day_of_year = 0 do while ( day_of_year.le.0 ) century = years_from_1600/100 day_of_year = days_from_1600 - years_from_1600*365 . - (years_from_1600 - 1)/ 4 . + (years_from_1600 + 99)/100 . - (years_from_1600 + 399)/400 - 1 * If we are 1600 then add one day if( years_from_1600.eq.0 ) then day_of_year = day_of_year + 1 end if * See if the leap days have taken us to a earlier year if( day_of_year.le.0 ) then years_from_1600 = years_from_1600 - 1 end if end do ***** We now have number of days from start of year and the year * Convert years back to start of century year = mod( years_from_1600, 100) ***** See if this is a leap year leap_year = .false. * ! we are at beginning of century if( year.eq.0 ) then if( mod(century,4).eq.0 ) leap_year = .true. else if( mod(year,4) .eq.0 ) leap_year = .true. end if ***** If the day of year is less than 60 then the leap years do no not * matter, if we are greater than or equal to 60, need to account * for the leap years * ! Dont worry about leap years IF( day_of_year.lt.60 ) THEN if( day_of_year.le.31 ) then month = 1 day = day_of_year * ! we are in February else month = 2 day = day_of_year - 31 end if * ! Need to account for leap years ELSE if( leap_year .and. day_of_year.eq.60 ) then month = 2 day = 29 else if( leap_year ) day_of_year = day_of_year - 1 ***** Now find month month = 2 do while ( day_of_year.gt. days_to_month(month) ) month = month + 1 end do month = month - 1 day = day_of_year - days_to_month(month) end if END IF ***** Now save the date date(1) = years_from_1600 + 1600 date(2) = month date(3) = day ***** Now convert the fraction of a day to hours, minutes and seconds date(4) = fraction*24.d0 date(5) = fraction*1440.d0 - date(4)*60.d0 seconds = 86400.d0*fraction - date(4)*3600.d0 . - date(5)*60.d0 ***** Thats all RETURN end subroutine XYZ_to_GEOD(rot_mat, site_pos, geod_pos ) * Routine to compute the geodetic latitudes, longitude and * elipsoidal height from the cartesian XYZ coordinates of * the site. The algorithm used is adopted from: * Heiskanen W A and H. Moritz, Physical Geodesy, W.H.Freeman * and company, San Francisco, 364, 1967. Chapter 5.3. * * The routine also returns the tranformation between XYZ and * local NEU system. * * RESTRICTION: User must supply the XYZ coordinates. If these * are not available the can be computed from geod_pos using: * * X = (N+h) cos(phi) cos(lambda) * Y = (N+h) cos(phi) sin(lambda) * Z = [(1-e**2)N+h] sin(phi) * * where e**2 = 2f-f**2; f is flattening * N**2 = a**/[1-(e*sin(phi))**2] * N is East-West Radius of curvature. * include 'const_param.h' * dlatdN - Derivative of latitudes with respect to the * - (horizontal) North component of position * - This derivative should not be confused with * - dNdlat (same symbol N but different meanings) * dNdlat - the derivative of the curvature (rad_curve) * - with respect to latitudes. * dXdlat - Derivative of the X coordinate wrt latitudes * dYdlat - Derivative of the Y coordinate wrt latitudes * dZdlat - Derivative of the Z coordinate wrt latitudes * eccsq - eccentricity squared computed from flattening * equ_rad - radial distance of site from rotation axis. * geod_pos(3) - the geodetic coordinates. The order of the * - values is: * - (1) - Geodetic co-latitudes (rads) * - (2) - Geodetic longitude (positive east) (rad) * - (3) - ellipsoidal height (m) * lat_i - approximate latitudes used in iteration * lat_p - approximate latitudes from previous iteration * long - longitude of the site. * h_i - approximate height used in iteration * h_p - approximate height from previous iteration * rad_lat - radius to be used in computing the geodetic * - latitudes * rad_curve - radius of curvature at the site (N above) * rot_mat(3,3) - Transformation matrix XYZ and NEU where * - N is -colatitudes, E is longitude and U is * - along ellipsoidal height * site_pos(3) - XYZ coordinates of the site (m) * tolerance - toleranace for convergence on geodetic * - coordinates. real*8 dlatdN, dNdlat, dXdlat, dYdlat, dZdlat, eccsq, equ_rad, . geod_pos(3), lat_i, lat_p, long, h_i, h_p, rad_lat, . rad_curve, rot_mat(3,3), site_pos(3), tolerance * niter - Number of iterations. EXITS if more than 50 * iterations are required. (Removes problem if * coordinates are NaN in IEEE floating point) integer*4 niter * converged - Indicate values have converged. logical converged data * ! Converge to 0.1 mm . tolerance / 0.0001d0 / ***** Start, get the latitudes and height by iteration. equ_rad = sqrt( site_pos(1)**2 + site_pos(2)**2 ) eccsq = 2.d0*earth_flat - earth_flat**2 * ! Set previous iteration values lat_p = atan2( site_pos(3), equ_rad) h_p = 0.d0 converged = .false. niter = 0 do while ( .not. converged ) * Get radius of curvature using previous latitudes estimate rad_curve = earth_rad / . sqrt(1.d0 - eccsq*sin(lat_p)**2 ) rad_lat = equ_rad * . ( 1.d0 - eccsq*rad_curve/(rad_curve+h_p) ) lat_i = atan2( site_pos(3), rad_lat) * ! Use cos lat formula if( abs(lat_i).lt. pi/4 ) then h_i = equ_rad/cos(lat_i) - rad_curve * ! Use sin lat formula else h_i = site_pos(3)/sin(lat_i) - (1.d0-eccsq)*rad_curve end if * Check for convergence if( abs(h_i-h_p) .lt. tolerance .and. * ! Converged . abs(lat_i-lat_p)*rad_curve.lt. tolerance ) then converged = .true. end if * Check for two many iterations niter = niter + 1 if( niter.gt.50 ) then write(*,'('' XYZ_to_GEOD ERROR: Failure to converge'')') converged = .true. end if * Save the latest values h_p = h_i lat_p = lat_i * ! iterating for latitudes and height end do ***** Save the final values long = atan2( site_pos(2),site_pos(1) ) * ! colatitudes geod_pos(1) = pi/2.d0 - lat_i * ! Add 2*pi if( long.lt.0 ) then geod_pos(2) = 2*pi + long else geod_pos(2) = long end if geod_pos(3) = h_i ***** Now do the transformation between XYZ and NEU * Compute derivative of the radius of curvature with to latitudes * MOD TAH 980216: Fixed error in dNdLat (derivative in error by factor * of two). Long standing bug pointed out by Dan Leibach at CfA. dNdlat = earth_rad*eccsq *sin(lat_i)*cos(lat_i) / . sqrt( (1.d0-eccsq* sin(lat_i)**2)**3 ) * Now do NORTH component (First do derivate wrt latitudes, then compute * the Northing derivative) dXdlat = (dNdlat * cos(lat_i) - . (rad_curve+h_i)* sin(lat_i) )*cos(long) dYdlat = (dNdlat * cos(lat_i) - . (rad_curve+h_i)* sin(lat_i) )*sin(long) dZdlat = (1.d0-eccsq)*dNdlat *sin(lat_i) + . ((1.d0-eccsq)*rad_curve + h_i)*cos(lat_i) dlatdN = sqrt(dXdlat**2 + dYdlat**2 + dZdlat**2 ) * Now do rotation matrix rot_mat(1,1) = dXdlat / dlatdN rot_mat(1,2) = dYdlat / dlatdN rot_mat(1,3) = dZdlat / dlatdN * Now do EAST component rot_mat(2,1) = -sin(long) rot_mat(2,2) = cos(long) rot_mat(2,3) = 0.d0 * Now do UP component rot_mat(3,1) = cos(lat_i) * cos(long) rot_mat(3,2) = cos(lat_i) * sin(long) rot_mat(3,3) = sin(lat_i) ***** Thats all return end CTITLE CHECK_ASCII subroutine check_ascii(string) * Routine to check a string for non-ascii characters. Any * non-ascii characters are replaced with the @ character * PASSED Variables * string - The string to checked. (ichar should be between 32-127 * for each character in the string. If non-ascii characters * found, all charcaters after 32 are set to blanks. character*(*) string * LOCAL variables * i - loop counter * ascii_ichar - Value of ascii character (nulls are set to blank) integer*4 i, ascii_ichar * non_ascii - Set true if non-ascii characters are found logical non_ascii ***** Set non_ascii false to start non_ascii = .false. do i = 1, len(string) ascii_ichar = ichar(string(i:i)) if( ascii_ichar.eq.0 ) then string(i:i) = ' ' ascii_ichar = 32 else if ( ascii_ichar.lt.32 .or. . ascii_ichar.ge.128 ) then string(i:i) = '@' non_ascii = .true. end if end do ***** If we had non-ascii characters then clear the end of the * string if( non_ascii .and. len(string).gt.32 ) then string(32:) = ' ' end if ***** Thats all return end CTITLE report_stat subroutine report_stat( typep, modulep, routinep, filep, . messagep, ierr ) * Routine to report on the status on gamit and globk processing. * The routine reports status messages which indicate where in the * processing we are; warnings which are are not fatal to the processing * but could be used as diagnostic of later problems; and fatal errors * which cause the program to stop running. * For all reports, the messsages are appended to files named * 'module' .status, .warn or .fatal depending on the type of error * These files are opened with append status and closed before the * routine ends (the later is to ensure that the UNIX file buffers * flushed so that the user can see the current status). * Standard examples of the use of this routine would be: * STATUS: * call report_stat('STATUS','MODEL','model',' ','Beginning process',0) * WARNING: * call report_stat('WARNING','MODEL','wrtout',' ', * . 'Site radius too large', 0) * call report_stat('WARNING','GLOBK','update_aprioris',line, * . 'IOSTAT error reading string ',ierr) * FATAL: * call report_stat('FATAL','SOLVE','invert',' ', * . 'Normal equations non-positive definite',0) * call report_stat('FATAL','GLOBK','read_glb_mar',glb_mar_file, * . 'IOSTAT error opening file ',ierr) * CLEAR: * call report_stat('CLEAR','GLOBK',' ',' ',' ',0) * ROUTINES used from HP1000 library: * trimlen - Length of string * casefold - Converts string to upper case * caseunfold - Converts string to lower case * systime - Returns system time (available on both HP and Sun) * * PASSED Variables * ierr - Error number associated with report. This value is * - printed for warnings and fatal errors when non-zero. integer*4 ierr * type - Type of message: Only three types are allowed, any * - other type will cause this routine to print a warning * - message. The types are (only first four characters * - checked: * - status -- General reporting of status for progress * - information * - warning -- Non-fatal error in the processing, could be * - diagonistic of later fatal errors * - fatal -- Fatal error, If ierr=0 then only the message * - is printed, otherwize the error number is * - printed with the message. * - clear -- Removes the status, warning and fatal files * for a given module name. * module - Name of module (main program) (e.g., model, solve) * subroutine - Name of the subroutine from which the report is * - made. * file - Name of the file or string being manipulated. This is * printed whenever a non-zero length string is passed. * It is included with (Name xxxxx) at the end of the line. * message - Message to be printed character*(*) typep, modulep, routinep, filep, messagep character*(256) type, module, routine, file, message * LOCAL VARIABLES * rep_date(5) - Yr, mon, day, hr, min that the report is made * trimlen - Library routine to returned the non-blank length * - of a string * len_mod - Length of the module string * len_sub - Length of the subroutine string * len_mes - Length of the message * len_fil - Length of the file name * len_typ - Length of the status type (variable type) * jerr - IOSTAT error writing to the status file * kerr - IOSTAT error writing to stdout * cerr - IOSTAT error closing the status file * rep_unit - Unit number assigned to the report file. If we have * - trouble with the opening the unit, rep_unit is set * - 6 so that output will goto to stdout. * i - Counter used to find / or start of string in the module * name integer*4 rep_date(5), trimlen, len_mod, len_sub, len_mes, . len_fil, len_typ, jerr, kerr, cerr, rep_unit, i * rep_sec - Seconds that the report is made real*8 rep_sec * rep_file - Name of the report file to which the messages are * - written. Generated from the module name with an an * - extent based on the type. (NOTE: Module is upper * - case when name generated). * mod_name - Names of the module for creating the file name * (Removes the full information if present) character*256 rep_file, mod_name character*6 prog_name ***** Make a copy of all passed variables. Stops us overwritting * a fixed string. type = typep module = modulep routine = routinep file = filep message = messagep **** Start: First get the system time so that we know when the report * was made call systime( rep_date, rep_sec ) * Remove the 19-- portiion of the date so that HP and Sun look the * same rep_date(1) = rep_date(1) - int(rep_date(1)/100)*100 * See what type of report we are to make (this saves * unnecessary opening of files). Casefold the type incase user * calls with the wrong case. Module and subroutine are set case for * "nice" looking output. call casefold( type ) call casefold( module ) call caseunfold( routine ) * Get the lengths of the strings. (1 is minumum length to stop fortran * error is string(1:0) is attempted.) if( ichar(module(1:1)) .eq. 0 ) module(1:1) = ' ' if( ichar(routine(1:1)).eq. 0 ) routine(1:1) = ' ' if( ichar(message(1:1)).eq. 0 ) message(1:1) = ' ' if( ichar(file(1:1)) .eq. 0 ) file(1:1) = ' ' if( ichar(type(1:1)) .eq. 0 ) type(1:1) = ' ' * Make sure all the strings contain ascii characters call check_ascii( module ) call check_ascii( routine ) call check_ascii( message ) call check_ascii( file ) call check_ascii( type ) len_mod = max(1,trimlen(module)) len_sub = max(1,trimlen(routine)) len_mes = max(1,trimlen(message)) len_fil = max(1,trimlen(file)) len_typ = max(1,trimlen(type)) * * Now strip out the module name from the end the module string. * (Look for /'s in the name in case full name passed) i = len_mod do while ( i.gt.1 .and. module(i:i).ne.'/' ) i = i - 1 end do if( module(i:i).eq.'/' ) i = i + 1 mod_name = module(i:len_mod) * The following will result in the report files being named 'GAMIT' * for the main GAMIT batch sequence, but individual names for the * non-batch modules (MAKEXP, MAKEX, FIXDRV, and utilities) and * for the /kf routines except for AUTCLN. Currently, the GAMIT * batch script erases the 'GAMIT' files but all the other modules * erase their files upon executation with a 'CLEAR' call of report_stat. if( mod_name(1:3).eq.'ARC'.or. . mod_name(1:6).eq.'YAWTAB'.or. . mod_name(1:5).eq.'MODEL'.or. . mod_name(1:6).eq.'AUTECL' .or. . mod_name(1:6).eq.'AUTCLN'.or. . mod_name(1:5).eq.'CFMRG'.or. . mod_name(1:5).eq.'SOLVE'.or. . mod_name(1:6).eq.'SCANDD' ) then prog_name = 'GAMIT' else prog_name = mod_name(1:6) endif * Set the error flags for errors in this routine to zero. At end we * check these are still zero. jerr = 0 kerr = 0 cerr = 0 ***** Report the appropriate type of messages. * STATUS report if( type(1:4).eq.'STAT' ) then * Generate the name the file to write the status do and get * an available unit number to open the file with status append. call rep_open( prog_name, '.status', rep_file, rep_unit) * Now write the status message if( len_fil.le.1 ) then write(rep_unit, 120, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) write(*, 120, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) 120 format('STATUS :',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a) else * if file name passed print this as well) write(rep_unit, 140, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) write(*, 140, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) 140 format('STATUS :',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a,' (Name ',a,')') end if **** Warning messages that are not fatal (can be for information or to * report problems with file access and/or string reading. else if ( type(1:4).eq.'WARN' ) then * Generate file name and open call rep_open( prog_name, '.warning', rep_file, rep_unit) * If the error flag is zero, assume the file name is not * needed if( ierr.eq. 0 ) then * If the file name is zero or 1 charcater long, just * print the message, otherwize output the file name * as well. if( len_fil.le.1 ) then write(rep_unit, 220, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) write(*, 220, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) 220 format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a) else * Print the file name as well. write(rep_unit, 230, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) write(*, 230, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) 230 format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a,' (Name ',a,')') end if else * Report the messsage and the file name and the error number write(rep_unit, 240, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr write(*, 240, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr 240 format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a,1x,a,' ERROR ',i5) end if ***** FATAL error messages: These are processed the same as warnings * but will cause the program to stop running. else if ( type(1:4).eq.'FATA' ) then * Generate file name and open call rep_open( prog_name, '.fatal', rep_file, rep_unit) * If the error flag is zero, assume the file name is not * needed if( ierr.eq. 0 ) then * Again check to see if the file name has been passed if( len_fil.le.1 ) then write(rep_unit, 320, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) write(*, 320, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) 320 format('FATAL :',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a) else * Print the file name as well write(rep_unit, 330, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) write(*, 330, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) 330 format('FATAL :',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a,' (Name ',a,')') end if else * Report the messsage and the file name and the error number write(rep_unit, 340, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr write(*, 340, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr 340 format('FATAL :',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': ',a,1x,a,' ERROR ',i5) end if * For the fatal errors we close the rep_file and stop the program * running close( rep_unit, iostat=cerr ) stop 'FATAL Error: Stop from report_stat' ***** See if we clearing the files else if( type(1:4).eq.'CLEA' ) then * Open each of the types of files and delete on closing if( len_mod.gt.1 ) then call rep_open( prog_name, '.status', rep_file, rep_unit) close(rep_unit, iostat=cerr, status='delete') call rep_open( prog_name, '.warning', rep_file, rep_unit) close(rep_unit, iostat=cerr, status='delete') call rep_open( prog_name, '.fatal', rep_file, rep_unit) close(rep_unit, iostat=cerr, status='delete') end if else * An unknown type of report has been called. First warn the user * that an unknown type has been attempted, and then print out the * message information as if it were a warning. * Generate file name and open call rep_open( prog_name, '.warning', rep_file, rep_unit) write(rep_unit, 420, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), routine(1:len_sub), . type(1:len_typ) write(*, 420, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), routine(1:len_sub), . type(1:len_typ) 420 format('WARNING:',3(i2.2),':',2(i2.2),':',f4.1,1x, . a,'/',a,': Unknown type of status. Error in call to ', . 'report_stat -- ',a) * Now print the standard warning format. if( ierr.eq. 0 ) then * See if the file name has been passed. if( len_fil.le.1 ) then write(rep_unit, 220, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) write(*, 220, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes) else * Print the file name as well. write(rep_unit, 230, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) write(*, 230, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil) end if else * Report the messsage and the file name and the error number write(rep_unit, 240, iostat=jerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr write(*, 240, iostat=kerr) rep_date, rep_sec, . module(1:len_mod), . routine(1:len_sub), . message(1:len_mes), . file(1:len_fil), ierr end if end if **** Now clean up before exiting. Close the report file so that it will * be updated and check if any errors occurred during the writing of * the messages close( rep_unit, iostat=cerr ) if( jerr.ne.0 ) then write(*,520) jerr, rep_file(1:max(1,trimlen(rep_file))) 520 format('WARNING: IOSTAT error ',i4,' occurred in ', . 'report_stat writing the status file ', a) end if if( kerr.ne.0 ) then write(*,540) kerr 540 format('WARNING: IOSTAT error ',i4,' occurred in ', . 'report_stat writing the status report to stdout') end if if( cerr.ne.0 ) then write(*,560) cerr, rep_file(1:max(1,trimlen(rep_file))) 560 format('WARNING: IOSTAT error ',i4,' occurred in ', . 'report_stat closing the status file ', a) end if **** Thats all return end CTITLE REP_OPEN subroutine rep_open( module, extent, rep_file, rep_unit ) * This routine will open the status/warning/fatal file with access * append. The unit number is searched to find an available unit * and if no unit can be found the unit is returned as 6 so that * stdout will be used. * PASSED VARIABLES * rep_unit - Unit number for the file. Returned by this * - routine integer*4 rep_unit * module - Name of the module calling report_stat * extent - extent to be given to file name (status, warning * - or fatal). * rep_file - Name given to output file (returned). character*(*) module, extent, rep_file * LOCAL VARIABLES * ierr - IOSTAT error when file opened. * trimlen - Length of non-blank portion of string integer*4 ierr, trimlen * unit_open - Set true by inquire if unit is already in use. logical unit_open **** First generate the file name rep_file = module(1:max(1,trimlen(module))) // extent * Now find an available unit. Units 60-100 are searched unit_open = .true. rep_unit = 60 do while ( unit_open .and. rep_unit.lt.100 ) inquire( unit=rep_unit, opened=unit_open) * If the unit is in use, try the next one if( unit_open ) rep_unit = rep_unit + 1 end do *** If the unit is still open, then we have a major problem. We can * find a unit number for output so warning user and set the unit * to 6 if( unit_open ) then write(*,120) extent(2:) 120 format('WARNING: Cannot find available unit number for ',a, . ' reports. Output to stdout') rep_unit = 6 else * Found unit number, so open file open(rep_unit, file=rep_file, status='unknown', . access='append', iostat=ierr) if( ierr.ne.0 ) then * Error opening the output file. Tell user of problem write(*,140) ierr, extent(2:), . rep_file(1:max(1,trimlen(rep_file))), . rep_unit 140 format('WARNING: IOSTAT error ',i4,' occurred opening ', . a,' file ',a,' on unit ',i3) rep_unit = 6 end if end if ***** Thats all return end CTITLE 'caseunfold' subroutine caseunfold( string ) * Routine to convert a string to lower case. * i - Loop counter * len_string - used length of string * trimlen - returns string length integer*4 i, len_string, trimlen * string - Used portion of string character*(*) string ***** Get length of string len_string = trimlen(string) * Now loop and convert do i = 1, len_string if( string(i:i).ge.'A' .and. string(i:i).le.'Z') then string(i:i) = CHAR( ICHAR(string(i:i))+32 ) end if end do **** Thats all return end ctitle systime subroutine systime( date, secs) * This routine will return the current system date and seconds. * The date is a five element array with year, month, day, hour, minute. * * Note: time is returned to nearest second at the moment. * Variables * date(5) - Year, month, day, hour, min, sec integer*4 date(5) * secs - Seconds part of the time. real*8 secs * Functions used: * IDATE ! Sun routine to return date as day, month, year * ITIME ! Returns time in hr, min, seconds * LOCAL VARIABLES * jdate(3) - Return from idate * jtime(3) - Return from itime. integer*4 jdate(3), jtime(3) external idate **** Get the date call idate(jdate) date(1) = jdate(3) date(2) = jdate(2) date(3) = jdate(1) **** Get the hours, min and seconds call itime(jtime) date(4) = jtime(1) date(5) = jtime(2) secs = jtime(3) **** Thats all return end LOGICAL FUNCTION KBIT(IARRAY,IBIT) C C KBIT C C 1. KBIT PROGRAM SPECIFICATION C C 1.1. KBIT is true if the IBIT-th bit in IARRAY is set (1), C false if it is not (0). C KBIT is designed to complement SBIT which sets or resets C bits identified in the same way. C Bits are numbered starting with 1 as the lowest-order bit C in the first word of IARRAY, 16 is the sign bit in the first C word of IARRAY, 17 is the low-order bit in the second word C of IARRAY, etc. C C 1.2. RESTRICTIONS - NONE C C 1.3. REFERENCES - FILE "SOLV2 = SOLV2 PURPOSE AND OVERALL STRUCTURE C C 2. KBIT INTERFACE C C 2.1. CALLING SEQUENCE: CALL KBIT(IARRAY,IBIT) C C INPUT VARIABLES: C integer*4 IARRAY(*) , ibit C - may or may not be an array in calling program. C - Variable in which the flag bits are located. C C IBIT = index of bit to test. Bits are numbered starting with C 1 as the lowest order bit in IARRAY(1), 16 is the sign C bit in IARRAY(1), 17 is the lowest order bit in IARRAY(2), C etc. C C C OUTPUT VARIABLES: C C KBIT = FUNCTION VALUE = TRUE if the indicated bit is 1 C = FALSE if the indicated bit is 0 C C C CALLING SEQUENCE VARIABLE SPECIFICATION STATEMENTS: C C C 2.2. COMMON BLOCKS USED C C integer*4 MASK(32) C C COMMON MSKCM CONTAINS THE STANDARD BIT MASK ARRAY C C TO C VARIABLES IN MSKCM TO WHICH KBIT GIVES A VALUE: NONE C C C FROM C VARIABLES IN MSKCM FROM WHICH KBIT GETS A VALUE: MASK C MASK(1)=1B, MASK(2)=2B ... MASK(16)=100000B (sign bit) C I.E. MASK(I) HAS I-TH BIT (SOLVE CONVENTION) SET C C C C C 2.3. DATA BASE ACCESSES C C C 2.4. EXTERNAL INPUT/OUTPUT C C 2.5. SUBROUTINE INTERFACE: C C CALLING SUBROUTINES: THIS IS A SOLV2 UTILITY C C CALLED SUBROUTINES: C UTILITIES IN SYSTEM: IAND(FUNCTION) C C 3. LOCAL VARIABLES C C C END EQUIVALENCE C C 4. CONSTANTS USED (description and references only. Data stmts below.) C C C 5. INITIALIZED VARIABLES C C C C C 6. LOCAL VARIABLES - SPECIFICATIONS AND DATA STATEMENTS integer*4 ia, ib C C 6.1. LOCAL VARIABLES - SPECIFICATION STATEMENTS C C 6.2. TEXT VARIABLES - SPECIFICATION,DATA STATEMENTS C C 6.3. INITIALIZED VARIABLES - DATA STATEMENTS data mask / 1, 2, 4, 8, 16, 32, . 64, 128, 256, 512, 1024, 2048, . 4096, 8192, 16384, 32768, 65536,131072, . 262144, 524288, 1048576, 2097152, 4194304, . 8388608, 16777216, 33554432, 67108864, . 134217728, 268435456, 536870912, . 1073741824,-2147483648 / C C C 7. PROGRAMMER: J.C.PIGG C LAST MODIFIED: C# LAST COMPC'ED 820423:19:54 # CMOD: JCP 1980 FEB 7: MOVED MASK FROM LOCAL TO COMMON MSKCM. CMOD: JCP 79 JUL 27: CREATED C C PROGRAM STRUCTURE C C 1. Decompose IBIT into an array index IA and a bit index IB. C IA = (IBIT+31)/32 IB = IBIT - (IA-1)*32 C C C 2. Test the appropriate bit C KBIT = IAND(IARRAY(IA),MASK(IB)) .NE. 0 C C C 3. Return to calling program. C END SUBROUTINE SBIT(IARRAY,IBIT,IVALUE) C C SBIT C C 1. SBIT PROGRAM SPECIFICATION C C 1.1. SBIT ... sets a bit in IARRAY indexed by IBIT to IVALUE (0 or 1). C Logical function KBIT can then test the bit (true if bit is 1). C The indexing scheme is described below under INPUT VARIABLES. C C 1.2. RESTRICTIONS - IF IVALUE IS NOT 0 IT IS ASSUMED TO BE 1. C - IF IVALUE IS OMMITTED IN CALLING SEQUENCE, 1 IS ASSUMED. C C 1.3. REFERENCES - FILE "SOLV2 = SOLV2 PURPOSE AND OVERALL STRUCTURE C C 2. SBIT INTERFACE C C 2.1. CALLING SEQUENCE: CALL SBIT(IARRAY,IBIT,IVALUE) C NOTE THAT IVALUE IS OPTIONAL ARG, IF OMITTED DEFAULTS TO 1. C C INPUT VARIABLES: C integer*4 IARRAY(*) , ibit, ivalue C -Array containing the bits to be set or reset. C - may or may not really be dimensioned in calling program. C C IBIT = index of bit to be set or reset. C Bit 1 is the lowest order bit of IARRAY(1), bit 16 is the sign C bit of IARRAY(1), bit 17 is the lowest order bit of IARRAY(2), C etc. This is the same indexing scheme employed in func. KBIT. C C IVALUE = 0 or 1 = value to which bit is to be set; C = 1 if IVALUE is ommitted from calling sequence. C IF NON-ZERO, 1 IS ASSUMED. C C C OUTPUT VARIABLES: C C C CALLING SEQUENCE VARIABLE SPECIFICATION STATEMENTS: C C C 2.2. COMMON BLOCKS USED C C integer*4 MASK(32) C C COMMON MSKCM CONTAINS THE SOLV2 STANDARD BIT MASK ARRAY C C TO: VARIABLES IN MSKCM TO WHICH SBIT GIVES A VALUE: NONE C C C FROM: VARIABLES IN MSKCM FROM WHICH SBIT GETS A VALUE: MASK C MASK(I), I=1,16, is a word in which only the I-th bit is set. C MASK(1)=1, MASK(16) = 100000B (SIGN BIT) C C C C C 2.3. DATA BASE ACCESSES C C C 2.4. EXTERNAL INPUT/OUTPUT C C 2.5. SUBROUTINE INTERFACE: C C CALLING SUBROUTINES: THIS IS A SOLV2 UTILITY C C CALLED SUBROUTINES: C UTILITIES IN SYSTEM: IAND(FUNCTION),IOR(FUNCTION),NOT(FUNCTION) C SOLVE-2 UTILITIES (%%S2US): IGNPS(INT.FN) C C 3. LOCAL VARIABLES C C IVAL = INTERNAL VALUE OF IVALUE C C C END EQUIVALENCE C C 4. CONSTANTS USED (description and references only. Data stmts below.) C C C 5. INITIALIZED VARIABLES C C C C C 6. LOCAL VARIABLES - SPECIFICATIONS AND DATA STATEMENTS C C 6.1. LOCAL VARIABLES - SPECIFICATION STATEMENTS integer*4 ia, ib, ival C C 6.2. TEXT VARIABLES - SPECIFICATION,DATA STATEMENTS C C 6.3. INITIALIZED VARIABLES - DATA STATEMENTS C data mask / 1, 2, 4, 8, 16, 32, . 64, 128, 256, 512, 1024, 2048, . 4096, 8192, 16384, 32768, 65536,131072, . 262144, 524288, 1048576, 2097152, 4194304, . 8388608, 16777216, 33554432, 67108864, . 134217728, 268435456, 536870912, . 1073741824,-2147483648 / C C 7. PROGRAMMER: J.C.PIGG C LAST MODIFIED: C# LAST COMPC'ED 820423:19:54 # C C TAH 87APR29: Replaced NARG with FTN77 routine PCOUNT C C JCP 81DEC16: INSTEAD OF NO-OP'ING WHEN IVAL IS NEITHER 0 NOR 1, C TREAT ANY NON-ZERO VALUE AS 1. I BELIEVE THIS MAKES C THE ROUTINE MORE CONVENIENT TO USE, AND THE RESULTS C BETTER DEFINED. C ALSO USE NARG INSTEAD OF IGNPS BECAUSE NARG IS IN C %%GSFC AND IGNPS IS NOT. THIS IS FOR THE CONVENIENCE C OF POTENTIAL NON-SOLV2 USERS. CC LAST COMPC'ED 811216:14:09 # C JCP 80 OCT 17: MADE IVALUE OPTIONAL AND DEFAULT TO 1. ONCE TOO MANY TIMES C I'VE HAD A CALLING PROGRAM SCREW UP BECAUSE I FORGOT TO C CODE IN THE 1 AS THIRD ARGUMENT! CMOD: JCP 79 Jul 27: Created. C C PROGRAM STRUCTURE C C 1. GET # PARAMETERS IN CALLING SEQUENCE AND SET INTERNAL IVALUE APPROPRIATELY C ival = ivalue IF(IVAL.NE.0) IVAL = 1 C C C 1. Decompose IBIT into an array index IA and a bit index IB. C IA = (IBIT+31)/32 IB = IBIT - (IA-1)*32 C C C 2. IF IVALUE is 1, force the appropriate bit to 1. C if( ival.eq.1) then iarray(ia) = ior(iarray(ia),mask(ib)) else iarray(ia) = iand(iarray(ia),not(mask(ib))) end if C C C 3. IF IVALUE is 0, force the appropriate bit to 0. C C IF$(IVAL.EQ.0) - SECOND TEST UNECESSARY 81DEC16. C ELSE IARRAY(IA) = IAND(IARRAY(IA),NOT(MASK(IB))) C C C 4. Return to calling program. C END CTITLE 'YMDHMS_TO_MJD' SUBROUTINE YMDHMS_to_MJD ( date, seconds, epoch ) * ------------------------- * ., 26JUL85 <871210.1317> * Author: T.Herring 3:53 PM THU., 10 JULY, 1986 * *$S "YMDHMS_to_MJD" *----------------------------------------------------------------------- * Routine to convert a calender date with hours, minutes and * seconds to a Modified Julian date. The calender date is * ordered as year, month, day of month, hours, and * minutes. These values are stored in a single I*2 array. * * If the year is greater than 1000 then the it is assumed to * contain the centuries. * * This routine is only valid for dates after 1600 Jan 0. * * CALLING SEQUENCE: * ================= * CALL ymdhms_to_JD ( date, seconds, epoch) * * WHERE: * date is a 5 element array containing year, month, day of month, * hour, and minutes. * (I*4 5 element array INPUT) * seconds is the seconds part of the epoch * (REAL*8 INPUT) * epoch is the JD with fractional days corresponding to date and * seconds. * (REAL*8 OUTPUT) * *---------------------------------------------------------------------- *$E * day - day of month * date(5) - 5 element date array with year, month, day, * - hours and minutes. * days_to_month(12) - number of days from start of year to * - each month * leap_days - the number of leap days we need to include * month - month of year * year - the years since 1900.0 * years_from_1600 - number of years since 1600 Jan 0. * days_from_1600 - number of days since 1600 Jan 0. integer*4 day, date(5), days_to_month(12), leap_days, month, . year, years_from_1600, days_from_1600 * epoch - The JD corresponding to date and seconds * fraction - Fraction of a day. We compute this separately * - to avoiod some rounding error when added to MJD * mjd - the computes MJD with fractional days * seconds - the seconds part of the epoch. real*8 epoch, fraction, mjd, seconds * leap_year - Indicates that this is a leap year logical leap_year data days_to_month / 0, 31, 59, 90, 120, 151, . 181, 212, 243, 273, 304, 334 / ***** START, Make sure year is from 1900 year = date(1) month = date(2) day = date(3) if( year.lt.1000 ) year = year + 1900 ***** Compute number of years from 1600 years_from_1600 = year - 1600 ***** Now compute number of leap days upto the start of the year * we are in (i.e., if the year we are in is a leap year do not * add the leap day yet) leap_days = (years_from_1600 - 1)/ 4 . - (years_from_1600 + 99)/100 . + (years_from_1600 + 399)/400 + 1 if( years_from_1600.eq.0 ) leap_days = leap_days - 1 ***** Now see if we are in leap year leap_year = .false. if( mod(years_from_1600, 4).eq.0 .and. . (mod(years_from_1600,100).ne.0.or. . mod(years_from_1600,400).eq.0) ) leap_year = .true. ***** Now compute number of days sinec 1600 days_from_1600 = years_from_1600*365 + leap_days + . days_to_month(month) + day ***** Add extra day if we are after Februrary and this is a leap year if( month.gt.2 .and. leap_year ) then days_from_1600 = days_from_1600 + 1 end if ***** Compute the mjd and add in the fraction of a day part * The MJD of 1600 Jan 0 is -94554 fraction = seconds/86400.d0 + date(5)/1440.d0 + date(4)/24.d0 mjd = -94554.d0 + days_from_1600 + fraction epoch = mjd ***** THATS ALL RETURN end CTITLE 'MJD_TO_YMDHMS' SUBROUTINE MJD_to_YMDHMS ( epoch, date, seconds ) * ------------------------ * ., 10JUL86 <871210.1317> * Author: T.Herring 11:28 AM THU., 10 JULY, 1986 * *$S "MJD_to_YMDHMS" *---------------------------------------------------------------------- * Routine to convert a modified Julian date (with fractions of a day) * or a full Julian date, to a calender data with hours and minutes, * and floating-point seconds (Real*8). * * NOTE: If a full Julian date is used the resolution of the seconds * will only be about 10 microseconds, a MJD should yield a resolution * of about 0.1 microseconds in the seconds. * * This routine is only valid for dates after 1600 Jan 0 (MJD -94554) * * CALLING SEQUENCE: * ================= * CALL MJD_to_YMDHMS ( epoch, date, seconds ) * * WHERE: * epoch is a modified or full julian date with possibly a fractional * part of a day. The two type of input are destinguished by * the full Julian date being greater than 2,000,000. * (REAL*8 INPUT) * date is an array containing the calender date with (full) year, * month of year, day of month, hours, minutes. Note the year * will be returned with the centuries added. * (I*4 5 element array OUTPUT) * seconds is the floating point seconds part of the MJD * (REAL*8 OUTPUT) * *---------------------------------------------------------------------- * *$E * century - Number of century from 1600 Jan 0. * date(5) - the calender date corresponding to * - 'epoch' resolution to the minute * day - day of month * day_of_year - Number of days from start of year * days_to_month(13) - number of days to start of each month * - in a non-leap-year * month - month of year * year - years since start of century * years_from_1600 - Number of years since 1600 Jan 0. * days_from_1600 - Number of days elapsed since 1600 Jan 0. * - (MJD -94554.0 Julian date 2305447.0) integer*4 century, date(5), day, day_of_year, days_to_month(13), . month, year, years_from_1600, days_from_1600 * epoch - the julian date or modified julian date * - to be converted to calender date * fraction - the fraction of a day part of MJD * mjd - epoch converted to a MJD * mjd_day - the whole number days in the mjd * seconds - the seconds part of the MJD (<60) real*8 epoch, fraction, mjd, mjd_day, seconds * leap_year - Indicates that this days is a leap year logical leap_year data days_to_month / 0, 31, 59, 90, 120, 151, 181, . 212, 243, 273, 304, 334, 365 / ***** START, convert full jd to mjd if we have too. * ! epoch must be Full julian date if( epoch.gt.2000000.d0 ) then * ! convert to MJD. mjd = epoch - 2 400 000.5d0 else mjd = epoch end if ***** Remove the fractional part of the mjd mjd_day = aint ( mjd ) fraction = mod ( mjd,1.d0 ) if( mjd.lt.0 .and. fraction.ne.0.d0 ) then mjd_day = mjd_day - 1 fraction = fraction + 1.d0 end if ***** Now convert MJD (even day) to date (year, month, day ) * Get number of days since 1600. days_from_1600 = mjd_day - (-94554.0d0) years_from_1600 = days_from_1600/365.d0 ***** Now compute day_of_year and years_from_1600 accounting for * leap years * ! Just to get us into loop day_of_year = 0 do while ( day_of_year.le.0 ) century = years_from_1600/100 day_of_year = days_from_1600 - years_from_1600*365 . - (years_from_1600 - 1)/ 4 . + (years_from_1600 + 99)/100 . - (years_from_1600 + 399)/400 - 1 * If we are 1600 then add one day if( years_from_1600.eq.0 ) then day_of_year = day_of_year + 1 end if * See if the leap days have taken us to a earlier year if( day_of_year.le.0 ) then years_from_1600 = years_from_1600 - 1 end if end do ***** We now have number of days from start of year and the year * Convert years back to start of century year = mod( years_from_1600, 100) ***** See if this is a leap year leap_year = .false. * ! we are at beginning of century if( year.eq.0 ) then if( mod(century,4).eq.0 ) leap_year = .true. else if( mod(year,4) .eq.0 ) leap_year = .true. end if ***** If the day of year is less than 60 then the leap years do no not * matter, if we are greater than or equal to 60, need to account * for the leap years * ! Dont worry about leap years IF( day_of_year.lt.60 ) THEN if( day_of_year.le.31 ) then month = 1 day = day_of_year * ! we are in February else month = 2 day = day_of_year - 31 end if * ! Need to account for leap years ELSE if( leap_year .and. day_of_year.eq.60 ) then month = 2 day = 29 else if( leap_year ) day_of_year = day_of_year - 1 ***** Now find month month = 2 do while ( day_of_year.gt. days_to_month(month) ) month = month + 1 end do month = month - 1 day = day_of_year - days_to_month(month) end if END IF ***** Now save the date date(1) = years_from_1600 + 1600 date(2) = month date(3) = day ***** Now convert the fraction of a day to hours, minutes and seconds date(4) = fraction*24.d0 date(5) = fraction*1440.d0 - date(4)*60.d0 seconds = 86400.d0*fraction - date(4)*3600.d0 . - date(5)*60.d0 ***** Thats all RETURN end CTITLE DECODE_OPTION subroutine decode_option( buffer, cont_types, num_types, conts, . default_cont ) * Routine to read a buffer line from the markov file and decode * the contribution types and set the corresponding bits in the * 'conts' array. * If the 'RESET' command is given then conts is set equal to * the default_cont. * The type may be proceeded with a minus sign to force a * contribution to be turned off. (An optional plus sign can * also be used, but no sign defaults to plus). * MOD TAH 961206: Added feature to ignore anything after ! or # in line. * * 08:47 PM WED., 18 Feb., 1987 * * bit_state - state of the bit to be set (set to zero if * - minus sign preceeds the type) * conts(1) - The contributions bit map. A bit is set for * - each contribution to be set * default_cont(1) - The default contribution pattern. Set when * - the RESET command is given. * dummy - dummy value used as place holder in READ_LINE * i - Loop counter * iel - index in cont_types which matches the next * - command in buffer. * ierr - Error flag from READ_LINE. Error generated * - by either by IOSTAT error on read, or -1 * - generated by reaching the end of the string. * indx - Index used to keep track of where we are * - the buffer (used by READ_LINE) * next_cont_len - length of the next_cont string found * num_types - number of contributions allowed. Gives the * - length cont_types array. This value is also * - used to compute number of words in the bit * - mapped conts when the default value needs to * - be assigned. * num_words - number of words in Default_cont (computed from * - num_types) * trimlen - HP function to return length of a string integer*4 bit_state, conts(*), default_cont(*), dummy, i, iel, . ierr, indx, next_cont_len, num_types, num_words, trimlen * buffer - line read from the Markov file. Should contain * - the list of contributions to be applied or * - reset (if minus sign used) * cont_types(1) - the array of types which are set. This * - array must have the entries in the same order * - which the bits will be set. character*(*) buffer, cont_types(*) * next_cont - The next contribution string from buffer. One * - extra character is added to allow for plus or * - minus sign. character*25 next_cont ***** Loop through the buffer getting all the contrubution names * until we run out of entries or an error (ierr.ne.0) ierr = 0 * ! Start at the beginning of the line indx = 1 do while ( ierr.eq.0 ) * Get the next string in buffer call read_line( buffer, indx, 'CH', ierr, dummy, next_cont) call casefold(next_cont) if( next_cont(1:1).eq.'#' .or. . next_cont(1:1).eq.'!' ) ierr = -1 * Only process if ierr is zero (no error) if( ierr.eq.0 ) then * strip + or - sign from start of next_cont and set * bit_state acordingly call get_bit_state( next_cont, bit_state) * Now see if we can find this contribution, if we cant * ignore entry unless it is RESET. * ! Set not found value iel = -1 i = 0 next_cont_len = min( len(cont_types(1)), * ! We use this . trimlen( next_cont ) ) * ! quantity so that shortened version * ! of the contributions can be used. * ! The first match will be used. do while ( i.lt.num_types .and. iel.lt.0 ) i = i + 1 if( cont_types(i)(1:next_cont_len) .eq. . next_cont(1:next_cont_len) ) iel = i end do * If we found next_cont then set bit in cont according to * bit state * ! Found if( iel.gt.0 ) then call sbit(conts,iel,bit_state) * ! Not found, check for reset else * ! or clear * ! Reset to num_words = (num_types+31)/32 if( next_cont(1:6).eq.'RESET ' ) then * ! default value. do i = 1,num_words conts(i) = default_cont(i) end do * ! RESET cont end if if( next_cont(1:4).eq.'ALL ' ) then do i = 1, num_words conts(i) = -1 end do end if * ! Clear values if( next_cont(1:6).eq.'CLEAR ' .or. . next_cont(1:5).eq.'NONE ' ) then do i = 1,num_words conts(i) = 0 end do * ! CLEAR conts end if * ! contribution found end if * ! no error reading buffer end if * ! looping until end of string end do ***** Thats all return end subroutine force_parms( cov_parm, sol_parm, num_parm, . cov_col, sol_col, nc, num_force, force_values, . dchi, avat, equ_gn, ipivot, scale ) * This routine will force parameters in a covariance matrix to have * specific values. The parameters to be forced are given in the * NC array. There are num_force of them and the values to be forced * to are given in the force_values array. If local array mapping * is to be then map_force should be called before this routine * setup the mapping of the arrays needed. * num_parm - Number of parameters in cov_parm. Cov_parm * - is assumed to be dimensioned to this size. * num_force - Number of parameters to be focesd to specific * - values * nc(num_force) - The list of parameters to be made equal. * ipivot(num_force) - The pivot elements for the matrix * - inversion integer*4 num_parm, num_force, nc(num_force), ipivot(num_force) * cov_parm(num_parm,num_parm) - Full covariance matrix to * - be modified. * sol_parm(num_parm) - The solution vector. The nc * - elements of this vector will be set equal. * cov_col(num_parm, num_force) - The columns of the * - covariance matrix for the parameters * - being forced. * sol_col(num_force) - The change in the parameter * - estimates needed to get the forced values. * force_values(num_force) - The values the forced parameters * - should take on. * avat(num_force,num_force) - The multiplication of * - the partials matrix of (00001000,000001000) * - and the covaraiance matrix. * equ_gn(num_parm,num_force) - Kalman gain matrix for * - getting solution. * scale(num_force) - Scaling vector for solution. (passed * - to invert_vis) * dchi - Change in Chi**2 real*8 cov_parm(num_parm,num_parm), sol_parm(num_parm), . cov_col(num_parm, num_force), sol_col(num_force), . force_values(num_force), avat(num_force,num_force), . equ_gn(num_parm,num_force), scale(num_force), dchi * LOCAL PARAMETERS * i,j,k - Loop counters integer*4 i,j,k * dsol, dcov - Summation variables for computing corrections * - to the solution vector and covariance matrix. real*8 dsol, dcov * T ***** Start, form the AVA matrix, where A is of the form: * y = Ax (y is vector of zero observations) * and * A = 0 0 0 1 0 0 0 0 0 0....... to num_parm * 0 0 0 0 0 0 0 1 0 0 * 0 0 0 0 0 0 0 0 1 0 * where the above form would set parameters 4,8, and 9. * For num_force parameters being equated, there are num_force * rows in A. * * The above form is much simpler to compute is we save the * columns of the covariance matrix for thhose parameters to * be forced first. do i = 1, num_parm do j = 1, num_force cov_col(i,j) = cov_parm(i,nc(j)) end do end do * Save the forced parts of the solution vector as well do j = 1, num_force sol_col(j) = sol_parm(nc(j)) end do * T **** Now do AVA * do i = 1, num_force do j = 1, num_force avat(i,j) = cov_col(nc(i),j) end do end do **** Now invert this matrix. (If "sort of equal" was desired we could * add value to diagonal now representing variance of y above). Pass * zero as number in solution vector, we dont want to multiply. * kgain below is dummy argument. call invert_vis(avat, equ_gn, scale, ipivot, num_force, . num_force, 0 ) **** Before continuing compute the change in Chi**2 due to condition dchi = 0.d0 do i = 1, num_force do j = 1, num_force dchi = dchi + sol_col(i)*avat(i,j)*sol_col(j) end do end do dchi = dchi/num_force * Now form the Kalman gain, equ_gn given by * T T -1 * equ_gn = VA (AVA ) * do i = 1, num_parm do j = 1, num_force * Do the multiply (could use VIS but stick to straight) * call dwmul(equ_gn(i,j), col_col(i,1), num_force, * . avat(1,j),1, num_force) equ_gn(i,j) = 0.d0 do k = 1, num_force equ_gn(i,j) = equ_gn(i,j) + cov_col(i,k)* . avat(k,j) end do end do end do **** Now get the change to the solution vector * * x = x - equ_gn*(force_values-sol_col) * n o * do i = 1,num_parm dsol = 0.d0 do j = 1, num_force dsol = dsol + equ_gn(i,j)*(force_values(j)-sol_col(j)) end do sol_parm(i) = sol_parm(i) + dsol end do * Now update the covariance matrix * * V = V - equ_gn*cov_col * n o do i = 1, num_parm do j = 1, num_parm * Do summation loop dcov = 0.d0 do k = 1, num_force dcov = dcov + equ_gn(i,k)*cov_col(j,k) end do cov_parm(i,j) = cov_parm(i,j) - dcov end do end do **** Thats all return end CTITLE MAP_FORCE subroutine map_force(first, icov_col, isol_col, iavat, . iequ_gn, iscale, iipivot, max_space, num_parm, . num_force) * This routine sets up the mapping for an Interger*4 array * to contain the variables we need for equating parameters. * * first - location in array of first available word * icov_col - Start of the cov_col matrix * isol_col - Start of the sol_col vector * iavat - Start of the avat matrix * iequ_gn - Start of the equ_gn matrix * iscale - Start of the scale vector * iipivot - Start of the ipivot vector * max_space - maximum number of words allowed to be used. * num_parm - NUmber of parameters in solution * num_force - number of conditions integer*4 first, icov_col, isol_col, iavat, iequ_gn, iscale, . iipivot, max_space, num_parm, num_force * LOCAL VARIABLES * final - Last word used. integer*4 final ***** Start the mapping. icov_col = first isol_col = icov_col + 2*num_parm*num_parm iavat = isol_col + 2*num_parm iequ_gn = iavat + 2*num_force*num_force iscale = iequ_gn + 2*num_parm*num_force iipivot = iscale + 2*num_force final = iipivot + num_force ***** See if have enough memory * ! Not enough memory. if( final.gt.max_space ) then write(*,100) final, max_space 100 format(' MAP_EQUATE ERROR: Not enough memory: ',i9, . ' Integer*4 words required,',/, . 19x,'and only ',i8,' words available.') stop ' MAP_EQUATE ERROR: Not enough memory' end if **** Thats all return end CTITLE 'INVERT_VIS' subroutine invert_vis( norm_eq, sol_vec, scale, ipivot, nsize, . ndim, nsol) c c c Matrix inversion taken from Kalman_lib.ftn on HP1000. * MOD TAH 980611: Changed to not scale the matrix if scale is passed * with -1 for all values from 1 to nsize. This is to handle * non-symetric, non-postive definate matrices. c * i - Loop counters * ndim - dimension of matrix * ipivot(ndim) - Pivot elements for inversion * nsol - Number of solutions to be computed * nsize - Size of matrix to be inverted integer*4 i, ndim, ipivot(ndim), nsol , nsize * * norm_eq(ndim,ndim) - Symetric matrix to be inverted * sol_vec(ndim,nsol) - Solution vector * scale(ndim) - Scale needed to normalize the matrix real*8 norm_eq(ndim,ndim), sol_vec(ndim,nsol), scale(ndim) * * no_scale -- Logical set true if we should not scale the matrix * (indicated by the scale factors being preset to -1) logical no_scale c * Check to see if we need to scale no_scale = .true. do i = 1, nsize if( scale(i).ne. -1.d0 ) no_scale = .false. end do c ***** Scale the matrx before inversion if( .not.no_scale ) then call dwvmv( norm_eq(1,1), ndim+1, scale, 1, nsize) c * Check the scale factors do i = 1, nsize if( scale(i).ne.0 ) then scale(i) = 1.d0/sqrt(abs(scale(i))) else scale(i) = 1.d0 norm_eq(i,i) = 1.d0 end if end do c ***** Now scale the matric call scale_matrix( norm_eq, sol_vec, scale, nsize, nsol, ndim) end if c ***** Now do the Gauss Elimination call Gauss_elim( norm_eq, sol_vec, ipivot, nsize, nsol, ndim) c ***** Descale the matrix if( .not. no_scale ) then call scale_matrix( norm_eq, sol_vec, scale, nsize, nsol, ndim) end if c ***** Thats all return c END CTITLE 'GAUSS_ELIM' subroutine gauss_elim( norm_eq, sol_vec, ipivot, nsize, nsol, . ndim) c c * Gauss ellimination routine taken from Kalman_lib.ftn on HP1000 c * i - Loop counters * ndim - dimension of matrix * nsize - Size of matrix to be inverted * icol, irow - Row and column counters * jcol,jrow - other row and solumns counters * ipivot(ndim) - Pivot elements for inversion * imab - Index to max pivot * imax, jmax - Row and colmum with max values * nsol - Number of solutions to be computed integer*4 i, ndim, nsize, icol, irow, jrow, ipivot(ndim), . imab, imax, jmax, nsol * c * norm_eq(ndim,ndim) - Symetric matrix to be inverted * pivot - Value of pivot element * sol_vec(ndim,nsol) - Solution vector * temp - Temporary value real*8 norm_eq(ndim,ndim), pivot, sol_vec(ndim,nsol), temp * c c **** Initialize the pivot counters do i = 1, nsize ipivot(i) = i end do c ***** Start looping on diagonal for searching for pivot do irow = 1, nsize icol = ipivot(irow) c * Find maximum value in this row call dwmab( imab, norm_eq(icol, irow), ndim, nsize-irow+1) imax = imab + irow - 1 c * Get the pivot pivot = norm_eq(icol, imax) c * See if pivot is too small if( abs(pivot).lt.1.d-12) then write(*,100) irow, pivot, 1.d-12 100 format(/' Pivot element too small for row ',i3,/, . ' Value is ',d15.6,', tolerancc is ',d15.6) end if c ***** Interchange the columns for the pivot jmax = ipivot(imax) * ! Interchange if( imax.ne.irow ) then ipivot(irow) = ipivot(imax) ipivot(imax) = icol c * Swap matrix elements call dwswp( norm_eq(1,irow),1, norm_eq(1,imax), 1, nsize) end if c ***** Zero out the elements of this column and form col of inverse norm_eq(icol,irow) = norm_eq(jmax,irow) norm_eq(jmax,irow) = 1.d0 c * Form the column call dwsmy(1.d0/pivot, norm_eq(1,irow),1, norm_eq(1,irow),1, . nsize) c do jrow = 1, nsize if( jrow.ne.irow) then temp = - norm_eq(icol,jrow) norm_eq(icol,jrow) = norm_eq(jmax, jrow) call dwpiv(temp, norm_eq(1,irow),1, norm_eq(1,jrow),1, . norm_eq(1,jrow),1, nsize) norm_eq(jmax,jrow) = temp/pivot end if end do c * Do solution do jrow = 1, nsol temp = - sol_vec(icol,jrow) sol_vec(icol,jrow) = sol_vec(jmax, jrow) call dwpiv( temp, norm_eq(1,irow),1, sol_vec(1,jrow),1, . sol_vec(1,jrow),1, nsize) sol_vec(jmax,jrow) = temp/pivot end do c * end do loop over rows end do c **** Now change sign of solution vectoe if( nsol.gt.0 ) then * MOD TAH 980609: Only change the sign of the used portion of the * solution vector (used to be single call for ndim*nsol) do jrow = 1, nsol call dwsmy( -1.d0, sol_vec(1,jrow),1, sol_vec,1, nsize) end do end if c ***** Thats all Return END CTITLE 'SCALE_MATRIX' subroutine scale_matrix( norm_eq, sol_vec, scale, nsize, nsol, . ndim) c c * Routine to scale the matrix before inversion c * i,j - Loop counters * ndim - dimension of matrix * nsol - Number of solutions to be computed * nsize - Size of matrix to be inverted integer*4 i,j, ndim, nsol , nsize * * norm_eq(ndim,ndim) - Symetric matrix to be inverted * sol_vec(ndim,nsol) - Solution vector * scale(ndim) - Scale needed to normalize the matrix real*8 norm_eq(ndim,ndim), sol_vec(ndim,nsol), scale(ndim) * c c **** Start by scaling the rowa do i = 1, nsize call dwsmy( scale(i), norm_eq(i,1), ndim, norm_eq(i,1), ndim, . nsize) do j = 1, nsol sol_vec(i,j) = sol_vec(i,j)*scale(i) end do end do c ***** Scale the columns do j = 1, nsize call dwsmy( scale(j), norm_eq(1,j),1, norm_eq(1,j),1, nsize) end do c ***** Thats all return END CTITLE 'DWVMV' subroutine dwvmv( from, incf, to, inct, iter) c c * Routine to copy one vector to another one c * if, it - Indexes in from and to vectors * i - Loop counte * incf - Increment for from * inct - Increment to * iter - Number of iterations integer*4 if, it, i, incf, inct, iter * c C real*8 c . from(incf,1) ! the from vector c ., to(inct,1) ! the two vector real*8 from(1), to(1) * c c ***** Copy if = 1 it = 1 do i = 1, iter to(it) = from(if) c to(1,i) = from(1,i) it = it + inct if = if + incf end do c c write(*,*) if,it, iter, incf,inct, to(1), from(1) c ***** Thats all return END CTITLE 'READ_LINE' subroutine read_line(line,indx,type,err,value,cvalue) c c Routine to read the next 'thing' in string LINE, starting c with character indx. The 'thing' is returned in VALUE, c and the type of VALUE is given by the user in TYPE: c c TYPE type c ---- ---- c 'CH' character c 'I2' integer*2 c 'I4' integer*4 c 'R4' real*4 c 'R8' real*8 c 'R16' real*16 c c The IOSTAT error is returned in ERR. indx is updated to point c at the next blank character. If type is 'CH', the result is c returned in CVALUE, otherwise in VALUE. c c c J.L. Davis 870729 Transported to Vax from HP. Added REAL*16 c feature. c T.A. Herring 890301 Removed I*2 references (except for value) c so that routine is more generic Fortran. c R. W. King 970111 Dimensioned value to 8 rather than 1 to allow c bounds checking character*(*) line, type, cvalue c * value(1) - Passed as shortest length so that all * - will match. integer*2 value(8) integer*4 indx, err, length, indx_end, num_words integer*4 i c integer*2 i2_value(8) integer*4 i4_value real*4 r4_value real*8 r8_value c real*16 r16_value ! Real*16 not supported on the convex. Subs c ! Real*8 as fix. real*8 r16_value * option_found - Indicates that the type passed was valid * - Input value is left unchanged. logical option_found c equivalence (i2_value, i4_value) equivalence (i2_value, r4_value) equivalence (i2_value, r16_value) equivalence (i2_value, r8_value) c c.... Initialize results num_words = 0 * ! Default no error err = 0 c c.... Determine length of string length = len(line) if( indx.le.0 ) indx = 1 * MOD TAH 951026: Check to make sure that indx does not get too long if( indx.gt.length ) indx = length c c.... Find next nonblank character do while (line(indx:indx) .eq. ' ' .and. indx .lt. length) c c.... Increment indx indx = indx + 1 c end do c c.... Did we reach the end of the line? At this point we do not know if c we exited the above line because we were are the end of the line c or because the next non-blank character was the last character in c the line. Check now by seeing if line(indx:indx) is blank. If it c we have reached end of line so set ERR=-1 (EOF) * ! We did reach end of line if ( line(indx:indx).eq.' ' .and. indx.eq.length ) then err = -1 * MOD TAH 921102: We are not going to do anything so say we found the * option, so that an error message is not generated. * MOD TAH 951026: Set the end index value and added check on indx itself. option_found = .true. indx_end = indx * ! We still have characters to else * ! read c c.... Find the next blank character indx_end = indx do while (line(indx_end:indx_end) .ne. ' ' . .and. indx_end .lt. length) c c.... Increment end indx indx_end = indx_end + 1 c end do option_found = .false. call casefold( type ) c c.... Do we want character data? if (type .eq. 'CH') then c c.... Simply assign character substrings cvalue = line(indx:indx_end) option_found = .true. c end if c c.... Integer*2? if (type .eq. 'I2') then c c.... Read into I2 variable i2_value(1) = 0 read (line(indx:indx_end),*,iostat=err) i2_value(1) if (err .eq. 0) num_words = 1 option_found = .true. c end if c c.... Integer*4? if (type .eq. 'I4') then c c.... Read into I4 variable i4_value = 0 read(line(indx:indx_end),*,iostat=err) i4_value if (err .eq. 0) num_words = 2 option_found = .true. c end if c c.... Real*4? if (type .eq. 'R4') then c c.... Read into R4 variable r4_value = 0 read(line(indx:indx_end),*,iostat=err) r4_value if (err .eq. 0) num_words = 2 option_found = .true. c end if c c.... Real*16? if (type .eq. 'R6') then c c.... Read into R16 variable r16_value = 0 read(line(indx:indx_end),*,iostat=err) r16_value if (err .eq. 0) num_words = 8 option_found = .true. c end if c c.... Real*8? if (type .eq. 'R8') then c c.... Read into R8 variable r8_value = 0 read(line(indx:indx_end),*,iostat=err) r8_value if (err .eq. 0) num_words = 4 option_found = .true. c end if end if c.... Transfer to VALUE, if necessary if (num_words .gt. 0) then do i = 1, num_words value(i) = i2_value(i) end do end if c c.... Update indx indx = indx_end c.... Make sure we founf option if( .not.option_found ) then call bad_option( type, 'READ_LINE') err = 1 end if c END CTITLE 'DWPIV' subroutine dwpiv(scalar, v1,inc1, v2,inc2, v3,inc3, iter) c c * Does v3 = scalar*v1 + v2 -- Pivot operationa c * i - Loop counter * inc1,inc2,inc3 - increments fpor each vector * i1,i2,i3 - pointers for each vector * iter - number of iterations integer*4 i, inc1,inc2,inc3, i1,i2,i3, iter * C real*8 c . v1(inc1,1), v2(inc2,1), v3(inc3,1) ! Vectors to be operated on * scalar - Scaler multiplier real*8 v1(*), v2(*), v3(*), scalar * c c ***** Set up pointers and start i1 = 1 i2 = 1 i3 = 1 do i = 1 ,iter v3(i3) = scalar*v1(i1) + v2(i2) c v3(1,i) = scalar*v1(1,i) + v2(1,i) i1 = i1 + inc1 i2 = i2 + inc2 i3 = i3 + inc3 end do **** Thats all return end CTITLE 'DWSWP' subroutine dwswp(v1, inc1, v2,inc2, iter) c c * Routine to swap elements in V1 and V2 c * inc1,inc2 - increment on V1, and v2 * i - loop counter * i1,i2 - position in V1 and v2 * iter - number of values to be searched integer*4 inc1,inc2, i, i1,i2, iter * c C real*8 c . v1(inc1,1),v2(inc2,1) ! vector to be swapped * temp - Copy for swapping real*8 v1(1), v2(1), temp * c c **** Loop swapping values i1 = 1 i2 = 1 do i = 1, iter temp = v1(i1) v1(i1) = v2(i2) v2(i2) = temp i1 = i1 + inc1 i2 = i2 + inc2 c temp = v1(1,i) c v1(1,i) = v2(1,i) c v2(1,i) = temp end do c ***** Thats all return END CTITLE 'DWSMY' subroutine dwsmy(scalar, v1, inc1, v2, inc2, iter) c c * Scalar multiply routine v2 = scalar*v1 c * inc1, inc2 - Increments for vectors * iter - number of elements * i1, i2 - Indices for vector one and two * i - Loop counter integer*4 inc1, inc2, iter, i1, i2, i * c C real*8 c . v1(inc1,1), v2(inc2,1) ! Input vectors * scalar - multipling scalar real*8 v1(1), v2(1), scalar * c c **** Do the multipy i1 = 1 i2 = 1 do i = 1, iter v2(i2) = scalar*v1(i1) c v2(1,i) = scalar*v2(1,i) i1 = i1 + inc1 i2 = i2 + inc2 end do c ***** Thats all return END CTITLE 'DWMAB' subroutine dwmab(scalar, v1, inc1, iter) c c * Returns index of largest absolute value in v1 c * inc1 - increment on V1 * i - loop counter * i1 - position in V1 * iter - number of values to be searched * scalar - index of maximum value integer*4 inc1, i, i1, iter, scalar * c C real*8 c . v1(inc1,1) ! vector to be checked * valmax - maxvalue found real*8 v1(*), valmax * c c ***** Find maximum value * MOD TAH 980609: Initialized valmax to -1 to ensure that a * value is found. valmax = -1 i1 = 1 do i = 1, iter if( abs(v1(i1)).gt.valmax ) then valmax = abs(v1(i1)) i1 = i1 + inc1 scalar = i end if c if( abs(v1(1,i)).gt.valmax ) then c valmax = abs(v1(1,i)) c scalar = i c end if end do c ***** Thats all return END CTITLE BAD_OPTION subroutine bad_option(option, name) c Writes a message to the CRT that we could not understand c the option passed to a "RW" subroutine. * option - OPtion or block name used character*(*) name, option ***** Write the message write(*,100) option, name 100 format(' *** Bad option ',a,' passed to routine ',A) end c...................................................................... c c used to be FTN4X,X SUBROUTINE FOURG(DATA,N,ISIGN,WORK,ICRT) C ONE DIMENSIONAL COMPLEX FAST FOURIER TRANSFORM FOR ANY NUMBER OF C POINTS. DATA IS COMPLEX OF LENGTH N, AN ARBITRARY INTEGER. C THE FOLLOWING TRANFORM REPLACES DATA IN STORAGE: C TRAN(K) = 1/SQRT(N) * SUM(DATA(J)*EXP(2*PI*i*ISIGN*(J-1)*(K-1)/N)) C SUMMED OVER ALL J FROM 1 TO N AND ALL K FROM 1 TO N. C ISIGN = +1 OR -1, THE DIRECTION OF THE TRANSFORM. C WORK IS A COMPLEX WORKING STORAGE ARRAY OF LENGTH N. C RUNNING TIME IS PROPORTIONAL TO N*(SUM OF PRIME FACTORS OF N). C ORIGINAL SUBROUTINE BY NORMAN BRENNER, MIT, NOVEMBER 1971. C C# LAST COMPC'ED 850520:1 C implicit none integer*4 ip0,n,isign,icrt,nrem,in,idiv,igo,ndone,i3,i2,i23 .,istep,i1,i231,i321,j2,imax,i * MOD TAH 971121: Changed the declaration from real*8 to real*4 * which is what it is suppposed to be. real*4 pi2,data(*),work(*),theta,wr,wi,temp,s c DIMENSION DATA(*),WORK(*) DATA PI2/6.283185308/ IP0=2 IF(N.GT.0 .OR. IABS(ISIGN).EQ.1) GO TO 20 WRITE(ICRT,30) N,ISIGN 30 FORMAT(/"ERROR IN FOURG. N = ",I10,", ISIGN = ",I4) RETURN 20 NREM=N IN=-1 IDIV=2 IGO=-1 GO TO 80 50 IDIV=1 IGO=0 60 IDIV=IDIV+2 IF(IDIV**2.LE.NREM) GO TO 80 IDIV=NREM IGO=1 IF(IDIV.LE.1) GO TO 160 80 IF(NREM.NE.IDIV*(NREM/IDIV)) GO TO 150 NDONE=N/NREM NREM=NREM/IDIV DO 140 I3=1,NDONE DO 145 I2=1,IDIV I23=((I2-1)*NDONE+I3-1)*NREM THETA=PI2*FLOAT(I23)/FLOAT(ISIGN*N) WR=COS(THETA) WI=SIN(THETA) ISTEP=NREM*IP0 DO 143 I1=1,NREM I231=(I23+I1-1)*IP0+1 I321=((I3*IDIV-1)*NREM+I1-1)*IP0+1 IF(IN.GT.0) GO TO 120 WORK(I231)=DATA(I321) WORK(I231+1)=DATA(I321+1) DO 110 J2=2,IDIV I321=I321-ISTEP TEMP=WORK(I231) WORK(I231)=WR*TEMP-WI*WORK(I231+1)+DATA(I321) WORK(I231+1)=WR*WORK(I231+1)+WI*TEMP+DATA(I321+1) 110 CONTINUE GO TO 143 120 DATA(I231)=WORK(I321) DATA(I231+1)=WORK(I321+1) DO 130 J2=2,IDIV I321=I321-ISTEP TEMP=DATA(I231) DATA(I231)=WR*TEMP-WI*DATA(I231+1)+WORK(I321) DATA(I231+1)=WR*DATA(I231+1)+WI*TEMP+WORK(I321+1) 130 CONTINUE 143 CONTINUE 145 CONTINUE 140 CONTINUE IN=-IN GO TO 80 150 IF(IGO.LT.0) GO TO 50 IF(IGO.EQ.0) GO TO 60 160 S=SQRT(FLOAT(N)) IMAX=IP0*N IF(IN.GT.0) GO TO 190 DO 180 I=1,IMAX,IP0 DATA(I)=DATA(I)/S DATA(I+1)=DATA(I+1)/S 180 CONTINUE RETURN 190 DO 200 I=1,IMAX,IP0 DATA(I)=WORK(I)/S DATA(I+1)=WORK(I+1)/S 200 CONTINUE RETURN END CTITLE GET_BIT_STATE subroutine get_bit_state( next_cont, bit_state ) * Routine to check the first character in next_cont and set * bit_state accordingly. If the first character is - then * bit_state is set to zero, if + or any other character then * bit_state is set to one. The + or - sign is stripped from * next_cont if either value is present * bit_state - Bit state to be set depending on + or - * - at start of string * i - Loop counter * string_len - length of the next_cont string integer*4 bit_state, i, string_len character*(*) next_cont ***** Initialize bit_state, and then see if first character * is minus. Change bit state if it is bit_state = 1 if( next_cont(1:1).eq.'-' ) then bit_state = 0 end if * Strip + or - from string if( next_cont(1:1).eq.'-' .or. next_cont(1:1).eq.'+' ) then string_len = len(next_cont) do i = 2, string_len next_cont(i-1:i-1) = next_cont(i:i) end do end if ***** Thats all return end