c Last Update to Flux Routines: 28 Apr 1994 c SUBROUTINE albedo(iy,m,id,ihr,min,poslat,poslong,solar_in, * albedo_payne,albedo_fresnel,ierr) c Peter A. Coppin c CSIRO Centre for Environmental Mechanics c email: coppin@enmech.csiro.au c version 1.0, 28-Apr-94 c this subroutine calculates the albedo of sea water using two methods; c c the first is via the model of PAYNE (1972) c Payne, R.E. (1972) "Albedo of the sea Surface" c J. Atm. Sci. 29: 959-970 c which is essentially a fit to data obtained over open water. c It uses only date, time, position, the sun's elevation and c the total short wave solar radiation as inputs. The latter is compared c to the maximum solar input at that sun angle with no atmosphere. This is c a simple parameterisation for cloudiness. The values from the model fit c well with observations in a statistical sense, but due to the simplicity c of the model individual values may be incorrect. c c the second method simply calculates the theoretical albedo for sea water c with properties encountered in TOGA-COARE using Fresnel theory c i.e. point radiation source, perfectly flat surface and c perfectly clear water. This provides an interesting theoretical c limit to albedo values. Real values will be higher in the middle of the day c and lower at the extremes of the day due to the presence of diffuse c radiation and waves on the sea surface. c NOTE: c the time used is important. It should be the centre time of any averaging c period. The routine gives an error if the sun is below the horizon, so c adjustments to the nominated time may be necessary. c INPUTS c iy - (integer) year between 50 and 99 (1950 and 1999) c m - (integer) month between 1 and 12 c id - (integer) day between 1 and end of month (<31!) c ihr - (integer) hour between 1 and 24 (this must be UTC!!) c min - (integer) minute between 1 and 60 (also UTC) c poslat - (real) latitude in decimal degrees, S is negative c poslong - (real) longitude in decimal degrees, E is positive c solar_in - (real) incoming solar short wave radiation in W/m^2 c c OUTPUTS c albedo_payne - (real) albedo from PAYNE model ( value between 0.0 and 1.0) c sometimes returned as nonsense value when error c albedo_fresnel - (real) albedo from Fresnel theory ( value between 0.0 and 1.0) c sometimes returned as nonsense value when error c ierr - error flag, see following comments c c users MUST check that the error code is clear (ierr=0) as there is a multitude of c ways for the routine to fail c c ERROR CODES c ierr=1 (from SUNVEC) year not between 50 and 99 c ierr=2 (from SUNVEC) month > 12 or day >31 c ierr=3 (from SUNVEC) month or day negative c ierr=4 solar zenith distance above 90 deg i.e. NIGHT TIME!! c ierr=5 atmospheric transmission negative of exceeds 1, check solar radiation c value or time of day used is not appropriate. c ierr=6 error in PAYNE model lookup table (routine should not get this far) c there is some double precision arithmetic in the SUNVEC subroutine which c should be OK on most FORTRAN77 compilers c pi = 4. * ATAN(1.) drad = pi / 180. ierr=0 albedo_payne=999.99 albedo_fresnel=999.99 c call solar angle calculator call SUNVEC(poslat,poslong,iy,m,id,ihr,min,SL,SM,SN,ierrsv) IF (ierrsv.ne.0) THEN ierr=ierrsv return END IF azimuth = ATAN2(-sm , sl) / drad IF( azimuth .lt. 0) THEN azimuth = 360 + azimuth END IF zenangle = ATAN2(SQRT(1 - sn ** 2) , sn) / drad if (zenangle .gt. 90.0)then ierr=4 return else elevation=90.0-zenangle c calculate fresnel limit albedo call fresnel(zenangle,albedo_fresnel) c calculate atmospheric transmissivity trans=solar_in/(1353.0*cos(zenangle*drad)) if(trans.ge.0.0.and.trans.le.1.0)then c call PAYNE model call PAYNE(trans,elevation,albedo_payne,ierrp) if(ierrp.ne.0)then ierr=6 return end if else ierr=5 return end if end if RETURN END c----------------------------------------------------------------------- SUBROUTINE PAYNE(trans,elevation,albedo,ierr) DIMENSION payneobs(21,46) DATA (payneobs(i,1),i=1,21)/0.061,0.062,0.072,0.087,0.115,0.163, *0.235,0.318,0.395,0.472,0.542,0.604,0.655,0.693,0.719,0.732, *0.73, 0.681,0.581,0.453,0.425/ DATA (payneobs(i,2),i=1,21)/0.061,0.062,0.07, 0.083,0.108,0.145, *0.198,0.263,0.336,0.415,0.487,0.547,0.595,0.631,0.656,0.67, *0.652,0.602,0.494,0.398,0.37/ DATA (payneobs(i,3),i=1,21)/0.061,0.061,0.068,0.079,0.098,0.13, *0.174,0.228,0.29, 0.357,0.424,0.498,0.556,0.588,0.603,0.592, *0.556,0.488,0.393,0.342,0.325/ DATA (payneobs(i,4),i=1,21)/0.061,0.061,0.065,0.073,0.086,0.11, * 0.15, 0.192,0.248,0.306,0.36, 0.407,0.444,0.469,0.48, 0.474, *0.444,0.386,0.333,0.301,0.29/ DATA (payneobs(i,5),i=1,21)/0.061,0.061,0.065,0.07, 0.082,0.101, *0.131,0.168,0.208,0.252,0.295,0.331,0.358,0.375,0.385,0.377, *0.356,0.32, 0.288,0.266,0.255/ DATA (payneobs(i,6),i=1,21)/0.061,0.061,0.063,0.068,0.077,0.092, *0.114,0.143,0.176,0.21, 0.242,0.272,0.288,0.296,0.3, 0.291, *0.273,0.252,0.237,0.266,0.22/ DATA (payneobs(i,7),i=1,21)/0.061,0.061,0.062,0.066,0.072,0.084, *0.103,0.127,0.151,0.176,0.198,0.219,0.236,0.245,0.25, 0.246, *0.235,0.222,0.211,0.205,0.2/ DATA (payneobs(i,8),i=1,21)/0.061,0.061,0.061,0.065,0.071,0.079, *0.094,0.113,0.134,0.154,0.173,0.185,0.19, 0.193,0.193,0.19, *0.188,0.185,0.182,0.18, 0.178/ DATA (payneobs(i,9),i=1,21)/0.061,0.061,0.061,0.064,0.067,0.072, *0.083,0.099,0.117,0.135,0.15, 0.16, 0.164,0.165,0.164,0.162, *0.16, 0.159,0.158,0.157,0.157/ DATA (payneobs(i,10),i=1,21)/0.061,0.061,0.061,0.063,0.067,0.072, *0.08, 0.092,0.107,0.125,0.136,0.141,0.145,0.145,0.145,0.144, *0.143,0.142,0.141,0.14, 0.14/ DATA (payneobs(i,11),i=1,21)/0.061,0.061,0.061,0.062,0.065,0.068, *0.074,0.084,0.097,0.111,0.121,0.127,0.13, 0.131,0.131,0.13, *0.129,0.127,0.126,0.125,0.122/ DATA (payneobs(i,12),i=1,21)/0.061,0.061,0.061,0.061,0.063,0.067, *0.074,0.082,0.091,0.102,0.11, 0.116,0.119,0.118,0.116,0.114, *0.113,0.111,0.11, 0.109,0.108/ DATA (payneobs(i,13),i=1,21)/0.061,0.061,0.061,0.061,0.062,0.064, *0.07, 0.076,0.085,0.094,0.101,0.105,0.107,0.106,0.103,0.1, *0.097,0.096,0.095,0.095,0.095/ DATA (payneobs(i,14),i=1,21)/0.061,0.061,0.061,0.06, 0.061,0.063, *0.067,0.072,0.079,0.086,0.093,0.097,0.098,0.097,0.092,0.088, *0.086,0.084,0.083,0.083,0.083/ DATA (payneobs(i,15),i=1,21)/0.061,0.061,0.061,0.06, 0.061,0.062, *0.065,0.07, 0.075,0.081,0.086,0.089,0.09, 0.088,0.084,0.08, *0.077,0.075,0.074,0.074,0.074/ DATA (payneobs(i,16),i=1,21)/0.061,0.061,0.061,0.06, 0.06, 0.061, *0.064,0.067,0.071,0.076,0.081,0.083,0.084,0.081,0.076,0.072, *0.069,0.067,0.066,0.065,0.065/ DATA (payneobs(i,17),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.061, *0.063,0.065,0.068,0.072,0.076,0.077,0.076,0.074,0.071,0.067, *0.064,0.062,0.061,0.061,0.061/ DATA (payneobs(i,18),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.061, *0.062,0.064,0.067,0.071,0.073,0.074,0.073,0.069,0.065,0.062, *0.06, 0.058,0.057,0.057,0.056/ DATA (payneobs(i,19),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, *0.061,0.062,0.065,0.068,0.069,0.069,0.068,0.065,0.061,0.058, *0.055,0.054,0.053,0.052,0.052/ DATA (payneobs(i,20),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, *0.06, 0.062,0.063,0.066,0.067,0.066,0.064,0.061,0.057,0.054, *0.051,0.05, 0.049,0.048,0.048/ DATA (payneobs(i,21),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, *0.06, 0.06, 0.062,0.065,0.065,0.063,0.06, 0.057,0.054,0.05, *0.047,0.046,0.045,0.044,0.044/ DATA (payneobs(i,22),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, *0.06, 0.06, 0.061,0.063,0.064,0.061,0.058,0.055,0.051,0.047, *0.044,0.042,0.041,0.04, 0.04/ DATA (payneobs(i,23),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.06, *0.059,0.06, 0.06, 0.062,0.062,0.059,0.056,0.052,0.049,0.045, *0.042,0.04, 0.039,0.038,0.038/ DATA (payneobs(i,24),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.06, 0.061,0.06, 0.057,0.054,0.05, 0.047,0.043, *0.039,0.036,0.034,0.033,0.033/ DATA (payneobs(i,25),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059, *0.059,0.059,0.06, 0.06, 0.059,0.056,0.053,0.049,0.045,0.041, *0.037,0.035,0.033,0.032,0.032/ DATA (payneobs(i,26),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.059,0.059,0.058,0.055,0.051,0.047,0.043,0.039, *0.035,0.033,0.032,0.031,0.031/ DATA (payneobs(i,27),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.059,0.058,0.057,0.054,0.05, 0.046,0.043,0.039, *0.035,0.032,0.031,0.03, 0.03/ DATA (payneobs(i,28),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.058,0.057,0.056,0.053,0.049,0.046,0.042,0.038, *0.035,0.032,0.03, 0.029,0.029/ DATA (payneobs(i,29),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.058,0.057,0.055,0.053,0.048,0.044,0.041,0.037, *0.034,0.031,0.029,0.028,0.028/ DATA (payneobs(i,30),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.059,0.058,0.057,0.055,0.052,0.048,0.044,0.04, 0.036, *0.033,0.03, 0.028,0.027,0.027/ DATA (payneobs(i,31),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.056,0.054,0.051,0.047,0.043,0.039,0.036, *0.033,0.03, 0.028,0.027,0.026/ DATA (payneobs(i,32),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.039,0.035, *0.032,0.03, 0.028,0.026,0.026/ DATA (payneobs(i,33),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.038,0.035, *0.032,0.03, 0.028,0.026,0.026/ DATA (payneobs(i,34),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.041,0.038,0.034, *0.032,0.029,0.027,0.026,0.026/ DATA (payneobs(i,35),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.041,0.037,0.033, *0.029,0.027,0.026,0.025,0.025/ DATA (payneobs(i,36),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.052,0.049,0.045,0.04, 0.036,0.032, *0.029,0.027,0.026,0.025,0.025/ DATA (payneobs(i,37),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.057,0.055,0.051,0.048,0.044,0.04, 0.036,0.032, *0.029,0.027,0.026,0.025,0.025/ DATA (payneobs(i,38),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.056,0.054,0.051,0.047,0.043,0.039,0.035,0.032, *0.029,0.027,0.025,0.025,0.025/ DATA (payneobs(i,39),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.057,0.056,0.053,0.05, 0.047,0.043,0.039,0.035,0.031, *0.028,0.026,0.025,0.025,0.025/ DATA (payneobs(i,40),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.056,0.054,0.05, 0.047,0.043,0.039,0.034,0.031, *0.028,0.026,0.025,0.025,0.025/ DATA (payneobs(i,41),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.031, *0.028,0.026,0.025,0.025,0.025/ DATA (payneobs(i,42),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, *0.027,0.026,0.025,0.025,0.025/ DATA (payneobs(i,43),i=1,21)/0.061,0.061,0.06, 0.06, 0.06, 0.059, *0.059,0.058,0.056,0.054,0.051,0.047,0.043,0.038,0.034,0.03, *0.027,0.026,0.025,0.025,0.025/ DATA (payneobs(i,44),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059, *0.058,0.057,0.056,0.054,0.05, 0.047,0.042,0.038,0.034,0.03, *0.027,0.026,0.025,0.025,0.025/ DATA (payneobs(i,45),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059, *0.058,0.057,0.056,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, *0.028,0.026,0.025,0.025,0.025/ DATA (payneobs(i,46),i=1,21)/0.061,0.061,0.06, 0.06, 0.059,0.059, *0.058,0.057,0.055,0.053,0.05, 0.046,0.042,0.038,0.034,0.03, *0.028,0.026,0.025,0.025,0.025/ ierr=0 if(elevation.ge.0.0.and.elevation.le.90.0)then ielev=nint(elevation/2.)+1 else ierr=1 end if if(trans.ge.0.0.and.trans.le.1.0)then itrans=nint(trans*20.)+1 else ierr=1 end if if(ierr.eq.0)then albedo=payneobs(itrans,ielev) else albedo=999.99 end if RETURN END c--------------------------------------------------------------------- subroutine fresnel(zenangle,albedo_fresnel) c this routine calculates the Fresnel limit albedo for sea water c as a funstion of solar zenith distance in degrees c for water at 35 ppt salinity and 30 deg C i.e.ref index of 1.339675 pi = 4. * ATAN(1.) drad = pi / 180. refindex=1.339675 alpha=zenangle*drad beta=asin(sin(alpha)/refindex) albedo_fresnel=0.5*((tan(alpha-beta)**2/tan(alpha+beta)**2)+ * (sin(alpha-beta)**2/sin(alpha+beta)**2)) return end c--------------------------------------------------------------------------------- SUBROUTINE SUNVEC(poslat,poslong,iy,m,id,ihr,min,SL,SM,SN,ierr) C SHELL 14.8.73 FOLLOWS GOODSPEED. C SUBROUTINE VERSION FOR CALCULATING VECTORS AS REQUIRED. C CALCULATES SUN'S VECTOR DIRECTION COSINES(SL,SM,SN) C TAKES IN TIME,DAY,MON,YEAR. CONV. TIME TO HR AND MIN C C RESULTS HAVE BEEN CHECKED AGAINST GOODSPEED'S PROGRAM C ON THE 3600, AND THE ANGLES GIVEN BY THIS PROGRAM C ON THE PDP-9 ARE CORRECT TO WITHIN 0.01 DEGREES. C MOD DE 14.8.73 C C DOUBLE PRECISION DAY,DAYT,DELTDY,TWOPI,T0,DAY2,A,ARG REAL LONG,LAT,NORTH C DIMENSION B(3),OM(3),EAST(3),NORTH(3),GL(3),GN(3),OVERH(3) DIMENSION SINTOM(3),RADVEC(3),MON(12) DATA GL(1),GL(2),GL(3)/-0.226715,-0.973961,0.0/ DATA GN(1),GN(2),GN(3)/0.89356453,-0.2080006,-0.3978419/ DATA OM(1),OM(2),OM(3)/0.3874825,-0.090196746,0.917454/ TWOPI=.6283185308D01 T0=60.0 * 24.0 * 1.90123508D - 06 PI=3.141593 ierr=0 LONG=poslong*TWOPI/360.0 LAT=poslat*TWOPI/360.0 SINLAT=SIN(LAT) COSLAT=COS(LAT) C C OLD SUBR GSDATE C SHELL 21.6.73 CALCULATES NUMBER OF DAYS FROM 4.1.50 FOR GSSUN C DO 3 I=1,12 3 MON(I)=31 HR=IHR RMIN=MIN C MON(4)=30 MON(6)=30 MON(9)=30 MON(11)=30 MON(2)=28 IF(IY.LT.50.OR.IY.GT.100) then ierr=1 return END IF Y=IY IF(ID.GT.31.OR.M.GT.12) then ierr=2 return END IF IF(ID.LE.0.OR.M.LE.0) then ierr=3 return END IF C NOW CALCULATE NUMBER OF DAYS SINCE 1.1.IY C IDAY=0 DO 10 K=1,12 IF(K.EQ.M) GO TO 20 10 IDAY=IDAY+MON(K) 20 IDAY=IDAY+ID C NOW ADD NUMBER OF DAYS FROM 4.1.50 TO 1.1.IY IIY=IY-50 IDAY=IDAY+365*IIY C ADD 1 DAY FOR EACH LEAP YEAR BETWEEN 1950 AND 19IY IF(IY.LT.52) GO TO 31 DO 35 I=52,IY,4 IF(I.EQ.IY.AND.M.LT.3) GO TO 31 IDAY=IDAY+1 35 CONTINUE C 31 IDAY=IDAY-4 DAY=IDAY HOUR=HR+RMIN/60.0 c change to zero hour zone offset to suit UTC DELTDY=(HOUR-(1.0+40.0/60.0))/24.0 c DELTDY=(HOUR-(1.0+10.0+40.0/60.0))/24.0 C ALLOWS FOR PERIHELION AT 0140HRS, 1950, AND HOURZONE. DAY=DAY+DELTDY DAYT=DAY C C OLD SUBR GSNUSS C SHELL 6.7.73 FOLLOWS GOODSPEED C CALCULATES VECTOR FROM SUN TO EARTH AT GIVEN TIME. C EC IS ECCENTRICITY OF EARTH'S ORBIT AROUND SUN. C C EC=0.016751 DAY2=DAY*T0 DAY2=DMOD(DAY2,1.0D0) C THIS TAKES THE REMAINDER ONLY !!! A=TWOPI*DAY2 C DAY2 IS EXCESS NUMBER OF DAYS IN CURRENT "YEAR". C NEWTON RAPHSON SOLUTION FOR E, THE ECCENTRIC ANOMOLY. N=0 C FIRST APPROXIMATION E=A SINEN=SIN(E) HOLD=SINEN COSEN=COS(E) 51 N=N+1 F=E-EC*SINEN-A F1=1.-EC*COSEN RF=F/F1 SINEN=-RF*COSEN+SINEN COSEN=RF*HOLD+COSEN E=E-RF RE=ABS(RF/E) IF(RE.GT.1.E-05) GO TO 51 C SINTH=SQRT(1.0-EC*EC)*SINEN SINTH=0.99985886*SINEN C THIS IS THE "TRUE VALUE" C* SINTH=.9997194*SINEN C THIS IS GOODSPEED'S VALUE--INCORRECT! C HOWEVER THE RESULTANT ERROR IN SUN'S ANGLES IS LESS THAN 0.01 DEG. COSTH=COSEN-EC FAC=1.0/SQRT(SINTH*SINTH+COSTH*COSTH) RADVEC(1)=FAC*COSTH C X-COMPONENT OF RADVEC, E.T.C...... RADVEC(2)=FAC*SINTH RADVEC(3)=0 C C OLD SUBR GSOLVC C SHELL 6.7.73 FOLLOWS GOODSPEED 5.4.72 CALCULATES SUN'S DIRECTIONAL COSINES IN LOCAL COORDINATE FRAME C CODING SLIGHTLY ALTERED TO GIVE CORRECT(!!!) RESULTS. C FOR GRIFFITH (NSW, AUSTRALIA) c but good enough for TOGA-COARE water albedos (PAC) C ARG=60.0*24.0*.4375269285D-02*DAYT+LONG+2.2312617D0 ARG=ARG/TWOPI ARG=DMOD(ARG,1.0D0) ARG=ARG*TWOPI COSA=COS(ARG) SINA=SIN(ARG) C C DO 70 I=1,3 OVERH(I)=COSA*GL(I)+SINA*GN(I) SINTOM(I)=SINLAT*OM(I) 70 OVERH(I)=OVERH(I)*COSLAT+SINTOM(I) CALL GSVEC(OM,OVERH,EAST,4) CALL GSVEC(EAST,B,EAST,8) CALL GSVEC(OVERH,EAST,NORTH,4) CALL GSVEC(RADVEC,NORTH,B,3) SL=-B(1) CALL GSVEC(RADVEC,EAST,B,3) SM=B(1) CALL GSVEC(RADVEC,OVERH,B,3) SN=-B(1) C C RETURN END FUNCTION ACOS(A) C COS -VE IN 2,3 QUAD.RETURN PI + ATAN(ANGLE) B=1.0-A*A IF (B) 5,6,6 6 IF (A) 1,2,3 1 ACOS= 3.14159 + ATAN(SQRT(1.0- A*A)/A) RETURN 2 ACOS = 3.14159/2.0 RETURN 3 ACOS = ATAN(SQRT(1.0-A*A)/A) RETURN 5 ACOS=999. C THIS INTENDED AS DIAGNOSTIC RETURN END C------------------------------------------------------ SUBROUTINE GSVEC(A,B,C,IOP) C SUBROUTINE TO PERFORM VECTOR OPERATIONS C SHELL 18.5.73 MODIFIED 8.6.73 CHECKED 12.6.73 C A,B,C ARE VECTORS, IOP IS VECTOR OPERATION C 1=ADDITION C 2=SUBTRACTION C 3=DOT PRODUCT C 4=CROSS PRODUCT C 5=MAGNITUDE C 6=MULTIPLICATION OF VECTOR BY A SCALAR C 7=COS(ANGLE)BETWEEN 2 VECTORS, RETURNING C(1) AS ANGLE, C AND C(2) AS COS OF THE ANGLE. C 8=NORMALIZATION OF VECTOR A, RETURNING UNDER SAME NAME C C C DIMENSION A(3),B(3),C(3) GO TO (1,2,3,4,5,6,7,8), IOP 2 DO 12 I=1,3 C(I)=A(I)-B(I) 12 CONTINUE RETURN 6 DO 16 I=1,3 C B(1) IS THE SCALAR 16 C(I)=B(1)*A(I) RETURN C C 7 C(2)=A(1)*B(1)+A(2)*B(2)+A(3)*B(3) C C(2) IS DOT PRODUCT C RMA AND RMB ARE MAGNITUDES OF A AND B RESPECTIVELY RMA=SQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3)) RMB=SQRT(B(1)*B(1)+B(2)*B(2)+B(3)*B(3)) C(2)=C(2)/(RMA*RMB) C C(2) IS COS(ANGLE) BETWEEN VECTORS A AND B C(3)=1.0E+38 C INDEX TO 38 FROM 70 FOR PDP 11/10 C THIS IS A DIAGNOSTIC CTH=C(2) C(1)=ACOS(CTH) C C(1) LIES BETWEEN ZERO AND PI (3.141593) RETURN C C 1 DO 11 I=1,3 C(I)=A(I)+B(I) 11 CONTINUE RETURN 3 C(1)=A(1)*B(1)+A(2)*B(2)+A(3)*B(3) C(3)=1.0E+38 C(2)=C(3) C THIS IS A DIAGNOSTIC WHICH GIVES NONSENSE VALUES C IF THE RESULT OF THE DOT PRODUCT IS USED AS A VECTOR. RETURN 4 C(1)=A(2)*B(3)-A(3)*B(2) C(2)=A(3)*B(1)-A(1)*B(3) C(3)=A(1)*B(2)-A(2)*B(1) RETURN 5 C(1)=SQRT(A(1)*A(1)+A(2)*A(2)+A(3)*A(3)) C(3)=1.0E+38 C(2)=C(3) C THIS IS A DIAGNOSTIC. RETURN 8 D=0.0 DO 18 I=1,3 D=D+A(I)*A(I) 18 CONTINUE IF(D.EQ.0.0) GO TO 80 D=SQRT(D) DO 19 I=1,3 A(I)=A(I)/D C(I)=A(I) 19 CONTINUE 80 RETURN END c----------------------------------------------------------------------