! Program 2.1 (WINDCHIL.F90) PROGRAM WindChill ! ! Program file name WINDCHIL.F90. ! Given a temperature (deg F) and wind speed (mi/hr), calculate ! the equivalent windchill temperature (deg F). ! ! Variable declarations... IMPLICIT NONE REAL :: temperature,wind_speed,wind_chill_temperature ! ! Get input... ! WRITE(*,*)'Give temperature in deg F and wind speed in mi/hr...' READ(*,*)temperature,wind_speed ! ! Do calculation... ! wind_chill_temperature=0.0817*(3.71*SQRT(wind_speed)+ & 5.81 - 0.25*wind_speed)*(temperature-91.4)+91.4 ! ! Display output... ! WRITE(*,*)wind_chill_temperature ! STOP END PROGRAM WindChill ! Program 3.1 (CIRCLE.F90) PROGRAM Circle ! ! Purpose: Calculate the circumference and area of a circle. ! ! variable declarations... ! IMPLICIT NONE !Forces explicit typing of all variables. REAL radius, circumference, area REAL, PARAMETER :: pi=3.1415927 ! ! get input... ! PRINT*,' Give radius of a circle: ' READ*,radius ! ! do calculations... ! circumference=2.0*pi*radius area=pi*radius*radius ! ! display output... ! PRINT*,' circumference=',circumference,' area=',area ! ! terminate program... ! END ! Program 3.1(a) (CIRCLE1A.F90) parameter (pi=3.1415927) print*,'Give radius of a circle: ';read*,r;c=2.0*pi*r;a=pi*r*r print*,'circumference=',c,' area=',a;end ! Program 3.2 (BOX.F90) PROGRAM box ! IMPLICIT NONE REAL height, width, length, SurfaceArea, volume ! PRINT *,'Give the height, width, and length of a box: ' READ *,height,width,length ! SurfaceArea=2.0*(height*width+height*length+length*width) volume=height*width*length ! PRINT *,' surface area = ',SurfaceArea,' volume = ',volume ! END ! Program 3.3 (NAMES.F90) PROGRAM names ! CHARACTER*20 name INTEGER age ! PRINT *,' What is your first name? ' READ *, name PRINT *,' How old are you? ' READ *,age PRINT *,name,age ! END ! Program 3.6 (MIXED.F90) PROGRAM mixed ! IMPLICIT NONE INTEGER i,j REAL x,y ! i=5 j=2 x=i/j ! PRINT *,i,j,x x=3.3 y=1.5 i=x/y PRINT *,i x=2/3*4. PRINT *,x ! END ! Program 3.7 (AVERAGE.F90) PROGRAM average ! IMPLICIT NONE REAL x,y,z,avg INTEGER n ! PRINT *,'Give three numbers: ' READ *,x,y,z n=3 ! avg=(x+y+z)/n ! PRINT *,'The average of these three numbers is ',avg ! END ! Program 3.11 (BEAM.F90) PROGRAM beam ! ! Calculate maximum deflection of a beam supported at both ends, ! with the load concentrated at the middle of the beam. ! IMPLICIT NONE REAL length, force, elasticity, moment_of_inertia,deflection ! PRINT *,' Give length (ft), force (lb): ' READ *,length,force PRINT *,' Give elasticity (lb/in^2), moment of inertia (in^4): ' READ *,elasticity,moment_of_inertia ! length=length*12.0 deflection=-force*length**3/(48.0*elasticity*moment_of_inertia) ! PRINT *,'The deflection (in) is: ',deflection ! END ! Program 3.12 (REL_MASS.F90) PROGRAM Electron ! ! Calculates relativistic mass and speed of an electron. ! See Schaum's Outline of Theory and Problems of College Physics. ! REAL rest_mass,relativistic_mass ! kg REAL voltage ! volts REAL speed ! m/sec REAL e ! electronic charge, Coulomb REAL c ! speed of light, m/sec PARAMETER (e=1.602e-19, c=2.9979e8, rest_mass=9.109e-31) ! PRINT *,' Give electron gun voltage in volts: ' READ *,voltage ! relativistic_mass=(voltage*e+rest_mass*c**2)/c**2 speed=c*(1.-(rest_mass/relativistic_mass)**2)**0.5 ! PRINT *,' relativistic mass (kg) and speed (m/s): ', & relativistic_mass,speed END ! Program 4.2 (STRING.F90) PROGRAM string ! ! Demo program for string operations. ! IMPLICIT NONE CHARACTER *10 first_name,last_name,name*15 ! first_name='Laura' last_name='Brooks' name=first_name//last_name PRINT *,name PRINT *,'untrimmed length: ',LEN(name) name=TRIM(first_name)//' '//TRIM(last_name) PRINT *,name PRINT *,'untrimmed length: ',LEN(name) PRINT *,'trimmed length: ',LEN_TRIM(name) ! END ! Program 4.3 (POLAR.F90) PROGRAM polar ! ! Convert polar coordinates to Cartesian and check the ! results by converting them back to polar coordinates. ! IMPLICIT NONE REAL X,Y,r,theta,pi,DegToRad ! pi=4.0*ATAN(1.0); DegToRad=pi/180.0 PRINT*,' Give polar coordinates r, theta (deg): ' READ*,r,theta ! X=r*COS(theta*DegToRad) Y=r*SIN(theta*DegToRad) PRINT*,'x and y: ',X,Y ! Recalculate values of r and theta... r=SQRT(X*X+Y*Y) theta=ATAN2(Y,X)/DegToRad ! PRINT*,'recalculated r and theta: ',r,theta ! END ! Program 4.4 (RING.F90) PROGRAM Ring ! ! Calculate the area of a circular ring. ! IMPLICIT NONE REAL Inner_Radius,Outer_Radius,pi,area,radius PARAMETER (pi=3.1415927) area(radius)=pi*radius*radius ! PRINT*,' Give outer radius, then inner radius: ' READ*,Outer_Radius,Inner_Radius ! PRINT*,' The area of the ring is: ', & area(Outer_Radius)-area(Inner_Radius) ! END ! Program 4.5 (REFRACT.F90) PROGRAM Refract ! ! Calculate angle of refraction for an incident ray, ! using Snell's Law. ! IMPLICIT NONE REAL ni,nr ! indices of refraction (dimensionless) REAL incident,refracted ! angles from perpendicular (deg) REAL pi,DegToRad PARAMETER (pi=3.1415927) ! DegToRad=pi/180. PRINT*,& ' Give indices of refraction for incident and refracting medium:' READ *,ni,nr PRINT*,' What is the angle of incidence? ' READ *,incident ! ! Convert refracted angle to degrees before displaying its value. ! refracted=ASIN(ni/nr*SIN(incident*DegToRad)) PRINT *,' refracted angle = ',refracted/DegToRad END ! Program 4.6 (HYPERBOL.F90) PROGRAM hyperbol ! ! Shows how to calculate inverse hyperbolic functions. ! IMPLICIT NONE REAL x,hyperbolic_sin,hyperbolic_cos,hyperbolic_tan REAL z,InvSinH,InvCosH,InvTanH InvSinH(z)=LOG(z+SQRT(z*z+1.0)) InvCosH(z)=LOG(z+SQRT(z*z-1.0)) InvTanH(z)=LOG((1.0+z)/(1.0-z))/2.0 ! PRINT*,'Give any real number: ' READ*,x hyperbolic_sin=SINH(x) hyperbolic_cos=COSH(x) hyperbolic_tan=hyperbolic_sin/hyperbolic_cos PRINT*,'Hyperbolic sin, cos, tan: ',hyperbolic_sin,hyperbolic_cos,& hyperbolic_tan PRINT*,'Inverse hyperbolic sin, cos, tan: ', & InvSinH(hyperbolic_sin), & SIGN(InvCosH(hyperbolic_cos),x),InvTanH(hyperbolic_tan) ! END ! Program 4.7 (NUMBERS.F90) PROGRAM Numbers ! Performs some tests on default REAL and INTEGER values. IMPLICIT NONE REAL x INTEGER i ! x=1. !can be any value of type REAL i=1 !can be any value of type INTEGER PRINT *,' For real numbers...' PRINT *,' Digits(x) ',Digits(x) PRINT *,' Huge(x) ',Huge(x) PRINT *,' Tiny(x) ',Tiny(x) PRINT *,' MaxExponent(x) ',MaxExponent(x) PRINT *,' MinExponent(x) ',MinExponent(x) PRINT *,' Precision(x) ',Precision(x) PRINT *,' Epsilon(x) ',Epsilon(x) PRINT *,' For integers...' PRINT *,' Digits(i) ',Digits(i) PRINT *,' Huge(i) ',Huge(i) END ! Program 5.1 (POLAR2.F90) PROGRAM polar2 ! ! Convert polar coordinates to Cartesian and check the ! result by converting them back to polar coordinates. ! Demonstrates formatted output. ! IMPLICIT NONE REAL X,Y,r,theta,pi,DegToRad ! pi=4.0*ATAN(1.0) DegToRad=pi/180.0 PRINT*,' Give polar coordinates r, theta (deg): ' READ*,r,theta ! X=r*COS(theta*DegToRad) Y=r*SIN(theta*DegToRad) ! PRINT 1000,X,Y ! ! Recalculate values of r and theta... ! r=SQRT(X*X+Y*Y) theta=ATAN2(Y,X)/DegToRad ! PRINT 1001,r,theta ! ! FORMAT statements... ! 1000 FORMAT(1x,'x and y: ',2f6.2) 1001 FORMAT(1x,'recalculated r and theta: ',2f7.3) ! END ! See Table 5.1 (DESCRIPT.F90) PROGRAM descript ! ! Demonstrate the use of format descriptors. ! IMPLICIT NONE INTEGER i ! PRINT 1000,17,3,-3,567 PRINT 1001,99 PRINT 1002,1.41,-1.41 PRINT 1003,2.9979e8 PRINT 1004,2.9979e8 PRINT 1005,-2.9979e8 PRINT 1006,'Fortran' PRINT 1007,'Fortran' PRINT 1008, 'Fortran' PRINT 1009,.TRUE. PRINT 1010,1994,1995,1996 ! carriage control (system may ignore) i=3 PRINT*,'No blank as first character.',i PRINT*,' Blank as first character.',i PRINT*,'++ sign as first character.',i PRINT*,'11 as first character.',i PRINT*,'00 as first character.',i ! 1000 FORMAT(4i4) 1001 FORMAT(i5.4) 1002 FORMAT(2f6.2) 1003 FORMAT(e13.4) 1004 FORMAT(3en13.4) 1005 FORMAT(3es13.4) 1006 FORMAT(a) 1007 FORMAT(a10) 1008 FORMAT(a4) 1009 FORMAT(L3) 1010 FORMAT(3g5.0) ! END ! Program 5.2 (POLAR3.F90) PROGRAM polar3 ! ! Convert polar coordinates to Cartesian and check the ! results by converting them back to polar coordinates. ! Demonstrates formatted output and creation of an output ! data file. ! IMPLICIT NONE REAL X,Y,r,theta,pi,DegToRad ! OPEN(1,file='polar3.out') ! pi=4.0*ATAN(1.0) DegToRad=pi/180.0 PRINT *,' Give polar coordinates r, theta (deg): ' READ *,r,theta WRITE(1,*)' For r and theta = ',r,theta ! X=r*COS(theta*DegToRad) Y=r*SIN(theta*DegToRad) ! PRINT 1000,X,Y WRITE(1,1000)X,Y ! ! Recalculate values of r and theta... ! r=SQRT(X*X+Y*Y) theta=ATAN2(Y,X)/DegToRad ! PRINT 1001,r,theta WRITE(1,1001)r,theta CLOSE(1) ! ! FORMAT statements... ! 1000 FORMAT(1x,'x and y: ',2f6.2) 1001 FORMAT(1x,'recalculated r and theta: ',2f7.3) END ! see Chapter Five (POLAR4.F90) PROGRAM polar4 ! ! Convert polar coordinates to Cartesian and check the ! results by converting them back to polar coordinates. ! Demonstrates formatted output and creation of an output ! data file. ! IMPLICIT NONE REAL X,Y,r,theta,pi,DegToRad ! OPEN(1,file='LPT1') ! pi=4.0*ATAN(1.0) DegToRad=pi/180.0 PRINT *,' Give polar coordinates r, theta (deg): ' READ *,r,theta WRITE(1,*)' For r and theta = ',r,theta ! X=r*COS(theta*DegToRad) Y=r*SIN(theta*DegToRad) ! PRINT 1000,X,Y WRITE(1,1000)X,Y ! ! Recalculate values of r and theta... ! r=SQRT(X*X+Y*Y) theta=ATAN2(Y,X)/DegToRad ! PRINT 1001,r,theta WRITE(1,1001)r,theta CLOSE(1) ! ! FORMAT statements... ! 1000 FORMAT(1x,'x and y: ',2f6.2) 1001 FORMAT(1x,'recalculated r and theta: ',2f7.3) END ! Program 5.3 (STARMAG.F90) PROGRAM StarMag ! ! Calculates absolute stellar magnitude based on relative magnitude ! and distance. ! IMPLICIT NONE REAL abs_mag,rel_mag,parsecs ! WRITE(6,"(' To calculate absolute stellar magnitude:')") WRITE(6,"(' Give relative magnitude and distance in parsecs: ')"& ,advance='no') READ *,rel_mag,parsecs ! abs_mag=rel_mag+5.0-5.0*LOG10(parsecs) ! WRITE(6,1000)rel_mag,parsecs,abs_mag ! 1000 FORMAT(' A star with relative magnitude ',f6.2,/& ' at a distance of ',f5.1,' parsecs'/& ' has an absolute magnitude of ',f6.2) END ! Program 5.4 (RELMASS3.F90) PROGRAM RelMass2 ! ! Calculates relativistic mass and speed of an electron. ! See Schaum's Outline of Theory and Problems of College Physics. ! IMPLICIT NONE REAL rest_mass,relativistic_mass ! kg REAL voltage ! volts REAL speed ! m/s REAL e ! electronic charge, Coulomb REAL c ! speed of light, m/s INTEGER u ! output unit PARAMETER (e=1.602e-19, c=2.9979e8, rest_mass=9.109e-31) u=6 ! PRINT *,' Give electron gun voltage in volts: ' READ *,voltage ! relativistic_mass=(voltage*e+rest_mass*c**2)/c**2 speed=c*SQRT(1.-(rest_mass/relativistic_mass)**2) ! WRITE(u,1000)voltage WRITE(u,1001)rest_mass WRITE(u,1002)relativistic_mass,speed,relativistic_mass/rest_mass WRITE(u,1003)speed/c ! 1000 FORMAT(' For an electron gun voltage of: ',es10.4,' V') 1001 FORMAT(' rest mass of electron: ',es10.4,' kg') 1002 FORMAT(' relativistic mass and speed: ',es10.4,' kg',& es12.4,' m/s'/& ' ratio of relativistic to rest mass: ',es10.4) 1003 FORMAT(' ratio of speed to speed of light: ',es10.4) END ! Program 6.3 (GRADES.F90) PROGRAM Grades ! ! Converts a numerical grade to a letter grade. Demonstrates ! compound IF... statements. ! IMPLICIT NONE REAL grade ! PRINT *,' Give a numerical grade 0-100: ' READ *,grade ! IF (grade .GE. 90.0) THEN PRINT *,' Letter grade: A' ELSE IF (grade .GE. 80.0) THEN PRINT *,' Letter grade: B' ELSE IF (grade .GE. 70.0) THEN PRINT *,' Letter grade: C' ELSE IF (grade .GE. 60.0) THEN PRINT *,' Letter grade: D' ELSE PRINT *,' Letter grade: F' END IF ! END ! Program 6.4 (GRADES2.F90) PROGRAM Grades2 ! ! Converts a numerical grade to a letter grade. Demonstrates ! SELECT CASE statements. ! IMPLICIT NONE REAL grade ! PRINT *,' Give a numerical grade 0-100: ' READ *,grade ! SELECT CASE (NINT(grade)) CASE (90:100) PRINT *,' Letter grade: A' CASE (80:89) PRINT *,' Letter grade: B' CASE (70:79) PRINT *,' Letter grade: C' CASE (60:69) PRINT *,' Letter grade: D' CASE DEFAULT PRINT *,' Letter grade: F' END SELECT ! END ! Program 6.5 (TRIGTABL.F90) PROGRAM TrigTabl ! ! Generate a table of trig values. Demonstrates count-controlled loops. ! IMPLICIT NONE REAL angle,deg_to_rad INTEGER i !loop counter ! deg_to_rad=4.0*ATAN(1.0)/180.0 PRINT *,' x sin(x) cos(x) tan(x)' PRINT *,' ------------------------------' ! DO i=0,180,5 angle=REAL(i)*deg_to_rad IF (i .NE. 90) THEN PRINT 1000,i,SIN(angle),COS(angle),TAN(angle) ELSE PRINT 1001,i,SIN(angle),COS(angle) END IF END DO ! 1000 FORMAT(1x,i3,3f9.4) 1001 FORMAT(1x,i3,2f9.4,' undef.') END ! Program 6.6 (ARCTAN.F90) PROGRAM tan_1 ! ! File name ARCTAN.F90. ! Uses conditional loop to estimate tan^1(x) from its series ! expansion. ! IMPLICIT NONE REAL x,term,arctan,sign INTEGER denominator REAL,PARAMETER :: error_limit=1e-7 ! PRINT *,' Give x, where x^2<1' READ *,x ! arctan=x term=x denominator=1 sign=1 PRINT *,' Intermediate values...' DO WHILE (ABS(term)>error_limit) sign=-sign denominator=denominator+2 term=sign*x**denominator/REAL(denominator) arctan=arctan+term PRINT *,denominator,term,arctan END DO PRINT *,' Estimated = ',arctan,' Intrinsic = ',ATAN(x) ! END ! Program 6.7 (CURE.F90) PROGRAM Cure ! ! Control temperature increases in a curing oven. ! INTEGER final_temperature,proposed_increase,current_temperature INTEGER, PARAMETER :: room_temperature=20 INTEGER counter ! PRINT *,' Give final oven temperature...' READ *,final_temperature counter=0 current_temperature=room_temperature PRINT *,' Current temperature is: ',current_temperature DO counter=counter+1 PRINT *,' Give proposed temperature increase...' READ *,proposed_increase IF ((current_temperature + proposed_increase) > & final_temperature) THEN PRINT *,' This increase is too large!' ELSE current_temperature = current_temperature + proposed_increase PRINT *,' Current temperature is: ',current_temperature PRINT *, ' This is how far you have to go: ', & final_temperature - current_temperature END IF IF (current_temperature >= final_temperature) EXIT END DO PRINT *, ' It took ',counter, & ' increases to reach the final oven temperature.' ! END ! Program 6.8 (LOOP.F90) PROGRAM loop ! IMPLICIT NONE INTEGER i,j ! DO i=1,5 WRITE(6,'(1x,i1)')i WRITE(6,'(1x)',advance='no') DO j=1,4 WRITE(6,'(i1)',advance='no')j END DO WRITE(6,*) END DO END ! Program 6.10 (REFRACT3.F90) PROGRAM refract3 ! ! Creates table of refracted angles for light ray in air incident on ! water, glass, and diamond. ! IMPLICIT NONE REAL air_index,water_index,glass_index,diamond_index REAL water_angle,glass_angle,diamond_angle REAL angle,DegToRad REAL ni,nr,incident,Refract ! for statement function INTEGER i PARAMETER (air_index=1.00,water_index=1.33,glass_index=1.50) PARAMETER (diamond_index=2.42) ! Function to calculate refracted angle... Refract(ni,nr,incident)=ASIN(ni/nr*SIN(incident)) ! DegToRad=4.0*ATAN(1.0)/180.0 WRITE(*,*)' Refracted angle........' WRITE(*,*)'incident water glass diamond' WRITE(*,1001)water_index,glass_index,diamond_index WRITE(*,*)'--------------------------------' DO i=0,90,10 angle=REAL(i)*DegToRad water_angle=Refract(air_index,water_index,angle)/DegToRad glass_angle=Refract(air_index,glass_index,angle)/DegToRad diamond_angle=Refract(air_index,diamond_index,angle)/DegToRad WRITE(*,1000)i,water_angle,glass_angle,diamond_angle END DO ! 1000 FORMAT(1x,i5,3f9.2) 1001 FORMAT(1x,'angle',3(3x,'(',f4.2,')')) END ! Program 6.11 (OSCILLAT.F90) PROGRAM Oscillat ! ! Generates resonant frequency table for LC circuit. ! IMPLICIT NONE REAL inductance ! Henrys INTEGER capacitance ! pico Farads REAL f ! frequency, kHz INTEGER row ! loop index INTEGER u REAL, PARAMETER :: pi=3.1415927 ! u=6 ! OPEN(u,file='oscillat.out') WRITE (u,"('Frequency, kHz')") WRITE (u,"(' C (pF)')") WRITE (u,"(' L (H)',10i7)") (capacitance,capacitance=2,20,2) WRITE (u,"('------------------------------------------& &-----------------------------------')") DO row=1,7 inductance=(row+1)*.0005 WRITE (u,'(f7.4)',ADVANCE='NO') inductance DO capacitance=2,20,2 ! NOTE: Express frequency in kHz by dividing by 1000. ! Convert pico Farads to Farads by multiplying by 10^-12. f=1./(2.*pi*SQRT(inductance*REAL(capacitance)*1e-12))/1000. WRITE (u,'(i7)',ADVANCE='NO') NINT(f) END DO WRITE (u,'()') ! line feed END DO CLOSE(u) ! END ! Program 6.12 (EXPOSE.F90) PROGRAM Expose ! ! Generate a random radiation exposure history for a sample ! so that the total exposure doesn't exceed a specified maximum. ! IMPLICIT NONE REAL current_exposure !the proposed current exposure REAL max_single !maximum single exposure REAL cum_exposure !cumulative exposure REAL max_total !maximum allowed total exposure INTEGER n_exposures !number of exposures REAL x !0<=x<1 INTEGER u INTEGER Count(1) !current value of system clock ! u=6 PRINT *,' What is the total allowed exposure (0-1000)?' READ *,max_total PRINT *,' What is the largest allowed single exposure (0-500)?' READ *,max_single WRITE(u,*) & ' max single and max total are ',max_single,max_total ! CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) CALL Random_Number(x) current_exposure=max_single*x cum_exposure=0.0 n_exposures=0 ! ! Start pre-test loop... DO WHILE ((cum_exposure+current_exposure)<=max_total) n_exposures=n_exposures+1 cum_exposure=cum_exposure+current_exposure WRITE(u,1000)x,n_exposures,current_exposure,cum_exposure ! Get a new exposure value to try... CALL Random_Number(x) current_exposure=Max_single*x END DO ! End pre-test loop... WRITE(u,1001)current_exposure ! 1000 FORMAT(f10.5,i3,2f8.1) 1001 FORMAT( & ' The next proposed exposure of ',f5.1,' is too large.'/& ' Terminate the experiment.') ! END ! Program 6.13 (BEAM2.F90) PROGRAM Beam2 ! ! Calculates beam deflection for four different support/loading systems. ! IMPLICIT NONE REAL elasticity !lb/in^2 REAL moment_of_inertia !in^4 REAL length !ft REAL load !lb REAL deflection !in INTEGER systemID !1 - supported at each end, concentrated load !2 - supported at each end, distributed load !3 - supported one end, concentrated at free end !4 - supported one end, distributed CHARACTER YesNo ! 10 PRINT *,' Give elasticity (lb/in^2) and moment of inertia (in^4).' READ *,elasticity, moment_of_inertia PRINT *,' Give the beam length in ft.' READ *,length PRINT *,' Choose one of these support/loading systems: ' PRINT *,' 1 - supported at each end, concentrated load' PRINT *,' 2 - supported at each end, uniformly distributed load' PRINT *,' 3 - supported at one end, concentrated load at free end' PRINT *,' 4 - supported at one end, distributed load' READ *,systemID SELECT CASE (systemID) CASE (1,3) PRINT *,' Give the concentrated force.' CASE (2,4) PRINT *,' Give the distributed weight.' CASE DEFAULT STOP 'Program terminated because of input error.' END SELECT READ *,load ! length=length*12.0 SELECT CASE (systemID) CASE (1) deflection=-load*length**3/(48.0*elasticity*moment_of_inertia) CASE (2) deflection=& -5.0*load*length**3/(384.0*elasticity*moment_of_inertia) CASE (3) deflection=-load*length**3/(3.0*elasticity*moment_of_inertia) CASE(4) deflection=-load*length**3/(8.0*elasticity*moment_of_inertia) END SELECT ! PRINT 1000,deflection PRINT *,'More? (y/n)' READ *,YesNo IF (YesNo=='y') GO TO 10 ! 1000 FORMAT(1x,es10.3) END ! Program 7.2 (CIRCLSUB.F90) PROGRAM CirclSub ! ! Calculate area and circumference of a circle, using ! a subroutine. ! IMPLICIT NONE REAL radius,area,circumference ! PRINT *,' What is the radius of the circle?' READ *,radius CALL CircleStuff(radius,area,circumference) ! PRINT 1000,area,circumference ! 1000 FORMAT(1x,2f10.3) END ! SUBROUTINE CircleStuff(radius,area,circumference) ! ! Do area and circumference calculations. ! IMPLICIT NONE REAL radius,area,circumference,pi INTENT(IN) radius INTENT(OUT) area,circumference PARAMETER (pi=3.14159) ! area=pi*radius*radius circumference=2.0*pi*radius ! END ! Program 7.4 (CIRCLSB2.F90) MODULE CircleSubs ! CONTAINS ! SUBROUTINE CircleStuff(radius,area,circumference) ! ! Do area and circumference calculations. ! IMPLICIT NONE REAL radius,area,circumference,pi INTENT(IN) radius INTENT(OUT) area,circumference PARAMETER (pi=3.14159) ! area=pi*radius*radius circumference=2.0*pi*radius ! END SUBROUTINE CircleStuff ! END MODULE !---------------------------------------------------- ! PROGRAM CirclSub2 ! ! File name CIRCLSB2.F90. ! Calculate area and circumference of a circle, using ! a subroutine. Uses MODULE to enforce intent. ! USE CircleSubs IMPLICIT NONE REAL radius,area,circumference ! PRINT *,' What is the radius of the circle?' READ *,radius CALL CircleStuff(radius,area,circumference) ! PRINT 1000,area,circumference ! 1000 FORMAT(1x,2f10.3) END ! Program 7.5 (CIRCLFUN.F90) MODULE CircleFunctions ! CONTAINS !--------------------------------- REAL FUNCTION Area(radius) ! ! Do area calculation. ! IMPLICIT NONE REAL, PARAMETER :: pi=3.14159 REAL, INTENT(IN) :: radius ! Area=pi*radius*radius ! END FUNCTION Area !------------------------------------------ REAL FUNCTION Circumference(radius) ! ! Do circumference calculation. ! IMPLICIT NONE REAL, PARAMETER :: pi=3.14159 REAL, INTENT(IN) :: radius ! Circumference=2.0*pi*radius ! END FUNCTION Circumference !--------------------------------- END MODULE !===================== PROGRAM CirclFun ! ! Calculate area and circumference of a circle, using ! two functions. ! USE CircleFunctions, ONLY : Area,Circumference IMPLICIT NONE REAL radius ! PRINT *,' What is the radius of the circle?' READ *,radius ! PRINT 1000,Area(radius),Circumference(radius) ! 1000 FORMAT(1x,2f10.3) END ! Program 7.8 (EXT_FUNC.F90) MODULE ExternalFunctions CONTAINS !-------------------------- REAL FUNCTION f_of_x(x) IMPLICIT NONE REAL, INTENT(IN) :: x f_of_x=x**2 END FUNCTION f_of_x !------------------------------------- SUBROUTINE Print_f(lo,hi,step,f) IMPLICIT NONE INTEGER, INTENT(IN) :: lo,hi REAL, INTENT(IN) :: step REAL x INTEGER i ! Explicit function interface definition... INTERFACE REAL FUNCTION f(x) REAL, INTENT(IN) :: x END FUNCTION f END INTERFACE ! x=0.-step PRINT *,' x f(x)' DO i=lo,hi x=x+step PRINT 1000,x,f(x) END DO 1000 FORMAT(1x,2f10.5) END SUBROUTINE Print_f !--------------------------------- SUBROUTINE f_of_x_sub(x,s) IMPLICIT NONE REAL, INTENT(IN) :: x REAL, INTENT(OUT) :: s s=SQRT(x) END SUBROUTINE f_of_x_sub !------------------------------------- SUBROUTINE Print_s(lo,hi,step,s) IMPLICIT NONE INTEGER, INTENT(IN) :: lo,hi REAL, INTENT(IN) :: step REAL x,y INTEGER i ! Explicit subroutine interface definition... INTERFACE SUBROUTINE s(x,y) REAL, INTENT(IN) :: x REAL, INTENT(OUT) :: y END SUBROUTINE END INTERFACE ! x=0.-step PRINT *,' x f(x)' DO i=lo,hi x=x+step CALL s(x,y) PRINT 1000,x,y END DO 1000 FORMAT(1x,2f10.5) END SUBROUTINE Print_s !--------------------------------- END MODULE ExternalFunctions !================================= PROGRAM Ext_Func ! ! Demonstrate how to pass a function or subroutine ! to a subprogram. ! USE ExternalFunctions, ONLY : f_of_x,f_of_x_sub, & Print_f,Print_s IMPLICIT NONE INTRINSIC SIN ! CALL Print_f(1,11,.1,SIN) CALL Print_f(1,11,.1,f_of_x) !uses a user-defined function CALL Print_s(1,11,.1,f_of_x_sub) !uses a subroutine ! END ! Program 7.9 (RELMASS3.F90) MODULE ElectronConstants IMPLICIT NONE REAL, PARAMETER :: rest_mass=9.109e-31 !kg REAL, PARAMETER :: e=1.602e-19 ! electronic charge, Coulomb REAL, PARAMETER :: c=2.9979e8 ! speed of light, m/s END MODULE ElectronConstants !--------------------------------- MODULE ElectronSubs ! CONTAINS ! SUBROUTINE Rel_E(voltage,rel_mass,speed) ! ! Calculates relativistic mass and speed of an electron. ! See Schaum's Outline of Theory and Problems of College Physics. ! USE ElectronConstants, Only : rest_mass,e,c IMPLICIT NONE ! REAL rel_mass ! kg REAL voltage ! volts REAL speed ! m/s INTENT(IN) voltage INTENT(OUT) rel_mass,speed ! rel_mass=(voltage*e+rest_mass*c**2)/c**2 speed=c*SQRT(1.-(rest_mass/rel_mass)**2) ! END SUBROUTINE Rel_E ! END MODULE ElectronSubs ! PROGRAM RelMass3 ! ! Driver for subroutine to calculate relativistic mass and speed ! of electron. ! USE ElectronSubs, Only : Rel_E USE ElectronConstants, Only : rest_mass,speed_of_light => c IMPLICIT NONE REAL voltage,relativistic_mass,speed ! PRINT *,' Give electron gun voltage in volts: ' READ *,voltage CALL Rel_E(voltage,relativistic_mass,speed) ! WRITE(6,1000)voltage WRITE(6,1001)rest_mass WRITE(6,1002)relativistic_mass,speed,relativistic_mass/rest_mass WRITE(6,1003)speed/speed_of_light ! 1000 FORMAT(' For an electron gun voltage of: ',es10.4,' V') 1001 FORMAT(' rest mass of electron: ',es10.4,' kg') 1002 FORMAT(' relativistic mass and speed: ',es10.4,' kg',& es12.4,' m/s'/& ' ratio of relativistic to rest mass: ',es10.4) 1003 FORMAT(' ratio of speed to speed of light: ',es10.4) END ! Program 7.10 (UNITS.F90) PROGRAM Units ! IMPLICIT NONE REAL Inches_to,PSI_to,Watts_to ! PRINT*,' Test Inches_to for 6"...' PRINT*,Inches_to('feet',6.0),Inches_to('ft',6.0),' feet' PRINT*,Inches_to('cm',6.0),Inches_to('centimeters',6.0),' cm' PRINT*,Inches_to('m',6.0),Inches_to('meters',6.0),' m' PRINT*,Inches_to('yd',6.0),Inches_to('yards',6.0),' yd' PRINT*,' Test PSI_to for 14.7 PSI...' PRINT*,PSI_to('Pa',14.7),PSI_to('Pascal',14.7),' Pa' PRINT*,PSI_to('mb',14.7),PSI_to('millibars',14.7),' mb' PRINT*,PSI_to('atm',14.7),PSI_to('atmospheres',14.7),' atm' PRINT*,PSI_to('cm-Hg',14.7),' cm-Hg' PRINT*,PSI_to('in-Hg',14.7),' in-Hg' PRINT*,' Test Watts_to for 1 kW...' PRINT*,Watts_to('Btu_per_s',1000.0),' Btu/s' PRINT*,Watts_to('cal_per_s',1000.0),' cal/s' PRINT*,Watts_to('ft-lb_per_s',1000.0),'ft-lb/s' PRINT*,Watts_to('hp',1000.0),Watts_to('horsepower',1000.0),' hp' END ! REAL FUNCTION Inches_to(what,value) ! IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: what REAL, INTENT(IN) :: value ! SELECT CASE (what) CASE ('ft','feet') Inches_to=value/12.0 CASE ('cm','centimeters') Inches_to=value*2.54 CASE ('m','meters') Inches_to=value*.0254 CASE ('yd','yards') Inches_to=value/36.0 CASE DEFAULT STOP 'ERROR in InchesTo: Requested conversion does not exist.' END SELECT END ! REAL FUNCTION PSI_to(what,value) ! IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: what REAL, INTENT(IN) :: value ! SELECT CASE (what) CASE ('Pa','Pascal') PSI_to=value/6895.0 CASE ('mb','millibars') PSI_to=value*6895.0/1.0e5*1000.0 CASE ('atm','atmospheres') PSI_to=value*6895.0/1.013e5*1000.0 CASE ('cm-Hg') PSI_to=value*6895.0/1333.0 CASE ('in-Hg') PSI_to=value*6895/1333.0/2.54 CASE DEFAULT STOP 'ERROR in PSI_to: The requested conversion is not available.' END SELECT END ! REAL FUNCTION Watts_to(what,value) IMPLICIT NONE CHARACTER*(*), INTENT(IN) :: what REAL, INTENT(IN) :: value ! SELECT CASE (what) CASE ('Btu_per_s') Watts_to=value/1054.0 CASE ('cal_per_s') Watts_to=value/4.184 CASE ('ft-lb_per_s') Watts_to=value/1.356 CASE ('hp','horsepower') Watts_to=value/746.0 CASE DEFAULT PRINT *,' The requested conversion is not available.' Watts_to=0.0 END SELECT END ! Program 7.11 (PLOTTER.F90) MODULE PlotSubs CONTAINS !------------------------------------ SUBROUTINE Plot(lo,hi,step,f,u) ! Create a crude plot of a function. IMPLICIT NONE REAL, INTENT(IN) :: lo,hi,step INTEGER, INTENT(IN) :: u !output unit REAL f,largest,smallest,value,scale REAL max_x INTEGER i,j,n,where PARAMETER (max_x=50.) !# of spaces for plot along x-axis n=NINT((hi-lo)/step) !# of points to plot smallest=1e10 largest=-1e10 ! Find smallest and largest values. DO i=0,n value=f(lo+REAL(i)*step) !evaluate the function IF (valuelargest) largest=value !and largest values END DO scale=max_x/(largest-smallest) !calculate scaling factor IF (u/=6) OPEN(u,file='c:\ftn90.dir\plot.out') DO i=0,n ! Where does the plotting symbol go? where=NINT((f(lo+REAL(i)*step)-smallest)*scale) ! Print the step. WRITE(u,1000,advance='no')i ! Print where-1 blank characters. WRITE(u,1010,advance='no')(' ',j=1,where-1) ! Finally, print the plotting character. WRITE(u,*)'*' END DO IF (u/=6) CLOSE(u) 1000 FORMAT(i3) 1010 FORMAT(50a1) RETURN END SUBROUTINE Plot !------------------------ END MODULE PlotSubs !======================== PROGRAM Plot_it ! ! Create crude plot of a user-specified function. ! USE PlotSubs, ONLY : Plot IMPLICIT NONE INTEGER u CHARACTER*1 where INTERFACE REAL FUNCTION F_of_x(x) REAL, INTENT(IN) :: x END END INTERFACE ! PRINT *,' Output to (s)creen or (f)ile? ' READ *,where IF (where=='f') THEN u=1 ELSE u=6 END IF CALL Plot(0.,10.,.5,F_of_x,u) END ! REAL FUNCTION F_of_x(x) REAL, INTENT(IN) :: x F_of_x=x**2-50. RETURN END ! see Chapter Seven (INTENT.F90) PROGRAM Intent_test ! ! MS-DOS name INTENT.F90 ! Check response of program to INTENT violations when subroutines ! are or are not included in a MODULE. ! USE One IMPLICIT NONE REAL x,y x=3.0 CALL Change(x,y) CALL Change2(x,y) PRINT *,x ! END ! SUBROUTINE Change(x,y) IMPLICIT NONE REAL, INTENT(IN) :: x REAL, INTENT(OUT) :: y ! x=2. !NAG F90 catches this violation, ! y=5. !but not this one. RETURN END ! Program 8.1 (RANDTEST.F90) PROGRAM RandTest ! ! Generate count histogram from random number generator. ! (Demonstrate 1-D arrays.) ! IMPLICIT NONE INTEGER A(0:9),i,index INTEGER Count(1) ! for random number generator REAL x ! ! Initialize array. ! A = 0 ! ! Generate count histogram. ! CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) DO 10 i=1,1000 CALL Random_Number(x) index = INT(x*10.) A(index) = A(index) + 1 10 CONTINUE ! PRINT*,A END ! Program 8.2 (SINCOS.F90) PROGRAM SinCos ! ! Illustrate the use of elemental functions for array assignments. ! IMPLICIT NONE REAL X(0:18),Y(0:18),Z(0:18),pi INTEGER i ! pi=4.*ATAN(1.) ! ! Create array of angles... DO i=0,18 X(i)=5*i*pi/180.0 END DO ! Calculate trig functions for entire array... Y=SIN(X) Z=COS(X) ! PRINT 1000,(i*5,Y(i),Z(i),i=0,18) 1000 FORMAT(1x,i3,2f6.3) END ! Program 8.3 (ROWCOL.F90) PROGRAM rowcol ! IMPLICIT NONE INTEGER B(3,4),row,col ! DO row=1,3 DO col=1,4 B(row,col)=row*10+col END DO END DO PRINT 1000,((B(row,col),col=1,4),row=1,3) PRINT * PRINT 1000,B 1000 FORMAT(1x,4i4) END ! Program 8.4 (OZONE.F90) PROGRAM OzoneData ! ! File name OZONE.F90. ! Process monthly ozone data ! (short version for 3 days with 5 hours in each day.) ! IMPLICIT NONE REAL ozone(31,24) INTEGER n_hours,n_days,hours,days DATA ozone(1,1:5)/3.2,3.3,3.1,2.9,3.4/, & ozone(2,1:5)/2.9,2.8,2.7,3.0,4.0/, & ozone(3,1:5)/2.8,2.6,2.5,2.9,3.1/, & n_hours,n_days/5,3/ ! ! Get daily averages... PRINT 1010,((ozone(days,hours),hours=1,n_hours), & SUM(ozone(days,1:n_hours))/n_hours,days=1,n_days) ! Get hourly averages... PRINT 1020,(SUM(ozone(1:n_days,hours))/n_days,hours=1,n_hours) ! 1010 FORMAT(1x,5f5.1,f5.2) 1020 FORMAT(1x,5f5.2) END ! Program 8.5 (UNION1.F90) MODULE SetStuff ! CONTAINS !-------------------------------- SUBROUTINE Union(A,B,C,n_A,n_B,n_C,n_set) IMPLICIT NONE INTEGER, INTENT(IN) :: n_A,n_B,n_C INTEGER, INTENT(OUT) :: n_set INTEGER, INTENT(IN) :: A(n_A) INTEGER, INTENT(IN) :: B(n_B) INTEGER, INTENT(OUT) :: C(n_C) INTEGER i,j LOGICAL add ! C=-1 !Initialize union "set" to negative value. C(1:n_A)=A !Now set union to A... n_set=n_A !and initialize union counter to # values in A. DO i=1,n_B !Search through B... add=.TRUE. j=1 DO WHILE ((j <= n_A) .AND. (add)) IF (B(i) == A(j)) add=.FALSE. !looking for duplicates... j=j+1 END DO IF (add) THEN !add new values if not duplicates... n_set=n_set+1 !increment union counter. C(n_set)=B(i) END IF END DO END SUBROUTINE Union !------------------------ END MODULE SetStuff !======================== PROGRAM UnionDemo ! ! Calculate union of two proper sets stored in arrays. ! USE SetStuff, ONLY : Union IMPLICIT NONE INTEGER :: A(10),B(10),C(20) INTEGER n_A,n_B,n_C,n_set,i ! n_A=4 n_B=3 n_C=n_A+n_B ! A=(/1,5,9,3,0,0,0,0,0,0/) B=(/5,3,8,0,0,0,0,0,0,0/) CALL Union(A,B,C,n_A,n_B,n_C,n_set) PRINT *,' Elements in union: ',n_set PRINT *,(C(i),i=1,n_C) ! END ! Program 8.6 (UNION2.F90) MODULE SetStuff ! CONTAINS !-------------------------------- SUBROUTINE Union(A,B,C,n_A,n_B,n_C,n_set) IMPLICIT NONE INTEGER, INTENT(IN) :: n_A,n_B,n_C INTEGER, INTENT(OUT) :: n_set INTEGER, INTENT(IN) :: A(n_A),B(n_B) INTEGER, INTENT(OUT) :: C(n_C) INTEGER i,j LOGICAL add ! C=-1 !Initialize union "set" to negative value. C(1:n_A)=A !Now set union to A... n_set=n_A !and initialize union counter to # values in A. DO i=1,n_B !Search through B... add=.TRUE. j=1 DO WHILE ((j <= n_A) .AND. (add)) IF (B(i) == A(j)) add=.FALSE. !looking for duplicates... j=j+1 END DO IF (add) THEN !add new values if not duplicate... n_set=n_set+1 !Increment union counter. C(n_set)=B(i) END IF END DO END SUBROUTINE Union !------------------------ END MODULE SetStuff !======================== PROGRAM UnionDemo ! ! Calculate union of two proper sets stored in arrays. ! USE SetStuff, ONLY : Union IMPLICIT NONE INTEGER, ALLOCATABLE :: A(:),B(:),C(:) INTEGER n_A,n_B,n_C,n_set ! n_A=4 n_B=3 n_C=n_A+n_B ALLOCATE(A(n_A),B(n_B),C(n_C)) ! A=(/1,5,9,3/) B=(/5,3,8/) CALL Union(A,B,C,n_A,n_B,n_C,n_set) PRINT *,' Elements in union: ',n_set PRINT *,C ! END ! Program 8.7 (MAT_SQ1.F90) MODULE MatrixStuff !- CONTAINS !- SUBROUTINE GetSquare(A,row,col) IMPLICIT NONE INTEGER r,c INTEGER, INTENT(IN) :: row,col REAL,INTENT(INOUT) :: A(row,col) ! PRINT *,'from GetSquare...' DO r=1,row PRINT *,(A(r,c),c=1,col) END DO A=A*A ! END SUBROUTINE GetSquare !- END MODULE MatrixStuff !=========================== PROGRAM Mat_sq1 ! USE MatrixStuff, ONLY : GetSquare IMPLICIT NONE REAL, ALLOCATABLE :: A(:,:) INTEGER row,col ! ALLOCATE(A(2,3)) A=RESHAPE( (/3.3,2.9,4.7,1.7,5.4,0.6/), (/2,3/) ) DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO CALL GetSquare(A,2,3) PRINT * DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO ! END ! Program 8.7(a) (MAT_SQ2.F90) MODULE MatrixStuff !- CONTAINS !- SUBROUTINE GetSquare(A,row,col) IMPLICIT NONE INTEGER r,c INTEGER, INTENT(IN) :: row,col REAL,INTENT(INOUT) :: A(row,col) ! PRINT*,'from GetSquare...' DO r=1,row PRINT *,(A(r,c),c=1,col) END DO A=A*A ! END SUBROUTINE GetSquare !- END MODULE MatrixStuff !=========================== PROGRAM Mat_Sq2 ! USE MatrixStuff, ONLY : GetSquare IMPLICIT NONE REAL A(5,5) INTEGER row,col DATA A/25*0/ ! A(1,1)=3.3 A(1,2)=4.7 A(1,3)=5.4 A(2,1)=2.9 A(2,2)=1.7 A(2,3)=0.6 DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO CALL GetSquare(A,2,3) PRINT * DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO ! END ! Program 8.7(b) (MAT_SQ3.F90) MODULE MatrixStuff !- CONTAINS !- SUBROUTINE GetSquare(A,row,col,maxrow,maxcol) IMPLICIT NONE INTEGER r,c INTEGER, INTENT(IN) :: row,col,maxrow,maxcol REAL, INTENT(INOUT) :: A(maxrow,maxcol) ! PRINT*,'from GetSquare...' DO r=1,row PRINT *,(A(r,c),c=1,col) END DO A=A*A ! END SUBROUTINE GetSquare !- END MODULE MatrixStuff !=========================== PROGRAM Mat_Sq3 ! USE MatrixStuff, ONLY : GetSquare IMPLICIT NONE REAL A(5,5) INTEGER row,col DATA A/25*0/ ! A(1,1)=3.3 A(1,2)=4.7 A(1,3)=5.4 A(2,1)=2.9 A(2,2)=1.7 A(2,3)=0.6 DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO CALL GetSquare(A,2,3,5,5) PRINT * DO row=1,2 PRINT *,(A(row,col),col=1,3) END DO ! END ! Program 8.8 (ALLOCAT.F90) PROGRAM Allocat ! ! Demonstrate syntax for allocatable arrays. ! IMPLICIT NONE REAL, ALLOCATABLE :: A(:,:) INTEGER i,j,n ! n=2 ALLOCATE(A(n,2*n)) DO i=1,n DO j=1,2*n A(i,j)=i*j END DO END DO PRINT *,((A(i,j),J=1,2*n),i=1,n) DEALLOCATE(A) ALLOCATE(A(n,3*n)) DO i=1,n DO j=1,3*n A(i,j)=i*j END DO END DO PRINT *,((A(i,j),J=1,3*n),i=1,n) ! END ! Program 8.10 (DEGDAYS.F90) PROGRAM DegDays ! ! Calculate heating degree days for a week. ! Demonstrate arrays of records. ! IMPLICIT NONE TYPE temp_data_type REAL high,low,heating_deg END TYPE temp_data_type INTEGER days,n_days PARAMETER (n_days=7) TYPE(temp_data_type) temp_data(n_days) REAL heating_deg_days,average ! ! Gather and process data. ! heating_deg_days=0. DO days=1,n_days PRINT 1000,days READ *,temp_data(days)%high,temp_data(days)%low temp_data(days)%heating_deg=65.-(temp_data(days)%high & +temp_data(days)%low)/2. IF (temp_data(days)%heating_deg < 0.) & !set neg. values to 0 temp_data(days)%heating_deg=0. heating_deg_days=heating_deg_days+temp_data(days)%heating_deg END DO average=heating_deg_days/n_days ! ! Display results ! PRINT *,' Day high low heating degrees' PRINT *,'-------------------------------' DO days=1,n_days WRITE(*,1010,advance='no') & days,temp_data(days)%high,temp_data(days)%low,& temp_data(days)%heating_deg IF (temp_data(days)%heating_deg > average) THEN WRITE(*,*)' *' ELSE WRITE(*,*) END IF END DO PRINT 1020,average,heating_deg_days ! 1000 FORMAT(' Enter high and low temperature (F) for day ',i1) 1010 FORMAT(1x,i4,2f5.1,f6.1) 1020 FORMAT(' Average heating degree days = ',f5.1/& ' Total heating degree days = ',f5.1) END ! Program 8.11 (VECTOROP.F90) !----------------------------------------------- MODULE VectorSubs CONTAINS ! SUBROUTINE Dot(a,b,DotProduct,angle) ! ! Calculate the scalar (dot) product of two 3-D vectors and the ! angle between them. ! IMPLICIT NONE INTEGER i REAL, INTENT(IN) :: a(3),b(3) REAL, INTENT(OUT) :: DotProduct,angle ! DotProduct=0. DO i=1,3 DotProduct=DotProduct+a(i)*b(i) END DO angle=ACOS(DotProduct/SQRT(a(1)**2+a(2)**2+a(3)**2)/ & SQRT(b(1)**2+b(2)**2+b(3)**2)) END SUBROUTINE Dot ! SUBROUTINE Cross(a,b,c) ! ! Calculate the cross product of two vectors. ! IMPLICIT NONE REAL, INTENT(IN) :: a(3),b(3) REAL, INTENT(OUT) :: c(3) ! c(1)=a(2)*b(3)-b(2)*a(3) c(2)=a(3)*b(1)-b(3)*a(1) c(3)=a(1)*b(2)-b(1)*a(2) END SUBROUTINE Cross ! END MODULE !------------------------------------------------- PROGRAM VectorOperations ! ! File name: VECTOROP.F90 ! Calculate and display vector dot and cross products. ! USE VectorSubs IMPLICIT NONE REAL DotProd,angle,pi INTEGER i REAL a(3),b(3),c(3) ! pi=4.*ATAN(1.) PRINT *,' Give 3 components for a, 3 for b: ' READ *,(a(i),i=1,3),(b(i),i=1,3) CALL Dot(a,b,DotProd,angle) PRINT 1000,DotProd,angle*180.0/pi CALL Cross(a,b,c) PRINT 1010,c ! 1000 FORMAT(1x,'dot product and angle (deg): ',2f8.3) 1010 FORMAT(1x,'cross product: ',3f8.3) END ! Program 8.12 (SIERPINS.F90) PROGRAM Sierpinski ! ! Program name SIERPINS.F90 ! One-dimensional cellular automata with rule that generates a ! Sierpinski triangle. ! IMPLICIT NONE INTEGER size PARAMETER (size=60) LOGICAL a(60),old_a(60),a_1,a0,a1 CHARACTER*1,b(60) INTEGER i,j,cycle,n_cycles DATA n_cycles,a,b/15,size*.FALSE.,size*' '/ ! cycle=0 a(size/2)=.TRUE. !start with a single live cell... b(size/2)='*' WRITE(*,1000)cycle,b ! ! Generate more cycles... ! DO i=1,n_cycles old_a=a DO j=2,size-1 a_1=old_a(j-1) a0=old_a(j) a1=old_a(j+1) a(j)=(a_1.AND.(.NOT. a0).AND.(.NOT. a1)) & .OR. ((.NOT. a_1).AND.a1) .OR. (a0.AND.a1) END DO DO j=1,size IF (a(j)) THEN b(j)='*' ELSE b(j)=' ' END IF END DO WRITE(*,1000)i,b END DO ! 1000 FORMAT(1x,'(',i2,')',60a1) END ! Program 8.13 (PROB.F90) MODULE ProbFunctions ! CONTAINS !----------------------------- INTEGER FUNCTION Fact(x) ! ! Calculate x! as long as x! not too large for default integer type. ! IMPLICIT NONE INTEGER, INTENT(IN) :: x INTEGER prod,i ! prod=1 DO i=2,x prod=prod*i END DO Fact=prod RETURN END FUNCTION Fact !---------------------------- INTEGER FUNCTION C(n,k) ! ! Calculate combinations of n things taken k at a time. ! IMPLICIT NONE INTEGER, INTENT(IN) :: n,k ! C=Fact(n)/Fact(k)/Fact(n-k) RETURN END FUNCTION C !----------------------------- END MODULE ProbFunctions !============================= PROGRAM prob ! ! Evaluate probabilities by parsing a string expression and ! performing the implied calculations. ! USE ProbFunctions, ONLY : C IMPLICIT NONE CHARACTER*80 a INTEGER i,n,k,length,sign,left REAL probability,prob_a ! PRINT*,'Type expression to be evaluated (no syntax checking).' PRINT*,'Example: 1-c(10,0,.2)-c(10,1,.2)' READ(*,1020)a ! length=LEN_TRIM(a) probability=0. IF (a(1:1)=='1') probability=1. sign=1 !a leading + sign is optional DO i=1,length IF (a(i:i)=='+') sign=1 IF (a(i:i)=='-') sign=-1 IF (a(i:i)=='(') left=i IF (a(i:i)==')') THEN READ(a(left+1:i-1),*)n,k,prob_a probability=probability+sign*C(n,k)*prob_a**k*(1.-prob_a)**(n-k) PRINT 1030,n,k,C(n,k) END IF END DO PRINT 1010,probability ! 1000 FORMAT(a1) 1010 FORMAT(1x,'probability = ',f10.4) 1020 FORMAT(a80) 1030 FORMAT(1x,'C(',i2,',',i2,') = ',i5) END MODULE One CONTAINS SUBROUTINE Change2(x,y) IMPLICIT NONE REAL, INTENT(IN) :: x REAL, INTENT(OUT) :: y x=2. !NAG F90 catches this violatoin ! y=5. !AND this one for subroutine in MODULE. END SUBROUTINE END MODULE One !-------------------------------- ! Program 9.2 (GRADES_F.F90) PROGRAM Grades_F ! ! Calculate average grade from a file of names and grades. ! Demonstrates use of sequential access text files. ! IMPLICIT NONE CHARACTER*6 name,heading*15 REAL grade,avg INTEGER n_grades DATA avg,n_grades/0.0,0/ ! OPEN(1,file='c:\ftn90\source\grades.dat',action='read') ! READ(1,1000)heading PRINT *,heading READ(1,1000)heading PRINT *,heading ! ! Start READ loop ! 10 READ(1,1010,END=999)name,grade WRITE(*,1020)name,grade avg=avg+grade n_grades=n_grades+1 GO TO 10 ! ! End READ loop ! 999 CLOSE(1) IF (n_grades > 0) THEN avg=avg/n_grades ELSE avg=0. END IF WRITE(*,1030)avg ! 1000 FORMAT(a15) 1010 FORMAT(a6,f9.1) 1020 FORMAT(1x,a6,f9.1) 1030 FORMAT(1x,'Average grade = ',f6.1) END ! Program 9.3 (AVG_TEMP.F90) PROGRAM avg_temp ! ! Process high and low temperatures for one month. ! IMPLICIT NONE TYPE hi_lo_type CHARACTER*8 date REAL hi,lo,daily_average END TYPE hi_lo_type ! TYPE(hi_lo_type) hi_lo(31) INTEGER n_days,d REAL monthly_average ! OPEN(1,file='c:\ftn90\source\january.dat') ! ! Initialize ... monthly_average=0.0 n_days=0 ! Read data file... 10 n_days=n_days+1 READ(1,1000,end=20)hi_lo(n_days)%date,hi_lo(n_days)%hi,hi_lo(n_days)%lo hi_lo(n_days)%daily_average= & (hi_lo(n_days)%hi+hi_lo(n_days)%lo)/2. monthly_average=monthly_average+hi_lo(n_days)%daily_average PRINT *,hi_lo(n_days)%date,hi_lo(n_days)%hi,hi_lo(n_days)%lo GO TO 10 ! End of READ loop... 20 CONTINUE CLOSE(1) ! Decrement loop counter by 1... n_days=n_days-1 ! Calculate average... monthly_average=monthly_average/n_days ! Print contents days for avg. T < monthly avg. T WRITE(*,1002)monthly_average WRITE(*,1003) WRITE(*,1005) DO d=1,n_days IF (hi_lo(d)%daily_average < monthly_average) & WRITE(*,1010)hi_lo(d)%date, & NINT(hi_lo(d)%hi),NINT(hi_lo(d)%lo), & hi_lo(d)%daily_average END DO ! 1000 FORMAT(a8,2f4.0) 1002 FORMAT(' Monthly average temperature = ',f6.2) 1003 FORMAT(' These days have below-average temperature:') 1005 FORMAT(1x,' Date hi low avg'/1x,'----------------------') 1010 FORMAT(1x,a8,2i4,f6.1) END ! Program 9.4 (BAROM.F90) PROGRAM barometer ! ! File name BAROM.F90 ! Interpret file of daily barometric pressure readings. ! IMPLICIT NONE CHARACTER*8 date REAL at_6am,at_noon,at_6pm,at_midnight CHARACTER*1 flag_6am,flag_noon,flag_6pm,flag_midnight ! OPEN(1,file='c:\ftn90\source\barom.dat') ! 10 READ(1,1000,end=20)date,at_6am,flag_6am,at_noon,flag_noon, & at_6pm,flag_6pm,at_midnight,flag_midnight WRITE(*,1010)date,at_6am,flag_6am,at_noon,flag_noon, & at_6pm,flag_6pm,at_midnight,flag_midnight GO TO 10 ! 20 CLOSE(1) ! 1000 FORMAT(a8,4(f6.0,a1)) 1010 FORMAT(1x,a8,4(f8.2,1x,a1)) END ! Program 9.5 (READTEST.F90) PROGRAM ReadTest ! ! Test program for formatted READ. ! IMPLICIT NONE INTEGER a,b REAL x,y ! OPEN(1,file='c:\ftn90\source\readtest.dat') READ(1,1000)a,b READ(1,1010)x READ(1,1020)y PRINT*, a,b PRINT*,x,y CLOSE(1) ! 1000 FORMAT(2i5) 1010 FORMAT(f6.0) 1020 FORMAT(f5.3) END ! Program 9.6 (BAROM2.F90) PROGRAM barom2 ! ! Interpret file of barometric pressure readings when data are missing. ! IMPLICIT NONE CHARACTER*8 date CHARACTER*6 temp REAL barom_value CHARACTER*1 flag CHARACTER*80 buffer !temporary string storage INTEGER i ! OPEN(1,file='c:\ftn90\source\barom2.dat') ! 10 READ(1,1000,end=20)buffer ! Read and write date... date=buffer(1:8) WRITE(*,1010,advance='no')date ! Check and interpret the four barometer fields... DO i=1,4 temp=buffer(3+7*i:9+7*i) IF (temp == '------') THEN WRITE(*,1020,advance='no') ELSE !do internal read on substring READ(temp,1025)barom_value,flag WRITE(*,1030,advance='no')barom_value,flag END IF END DO WRITE(*,1040) !Write a "line feed" GO TO 10 ! 20 CLOSE(1) ! 1000 FORMAT(a80) !to hold the temporary string 1010 FORMAT(1x,a8) !to read the date 1020 FORMAT(1x,'-------') !output for missing barometer field 1025 FORMAT(f5.0,a1) !for internal read on substring 1030 FORMAT(f6.2,1x,a1) !output for data 1040 FORMAT() !line feed after advance='no' END ! Program 9.7 (SMOOTH.F90) MODULE Forecasting ! IMPLICIT NONE REAL a(100),b(100) REAL b1,b2,b3,rand DATA b1,b2,b3/5.,15.,50./ INTEGER size CONTAINS !------------------------- SUBROUTINE ExpSmooth(weight) ! ! Apply exponential smoothing with specified weight. IMPLICIT NONE REAL, INTENT(IN) :: weight INTEGER i ! b(1)=a(1) DO i=2,size b(i)=b(i-1)+weight*(a(i-1)-b(i-1)) END DO END SUBROUTINE ExpSmooth !----------------------------- SUBROUTINE MakeData ! ! Create array = b1+b2*x+b3*random ! INTEGER Count(1),i ! CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) ! size=20 DO i=1,size CALL Random_Number(rand) a(i)=b1+b2*REAL(i)+b3*(rand-.5) ! PRINT 1000,a(i) END DO 1000 FORMAT(f8.3) END SUBROUTINE MakeData !---------------------------- END MODULE Forecasting !============================ PROGRAM Smooth ! ! Demonstrate exponential smoothing. ! USE Forecasting, ONLY : ExpSmooth,MakeData,a,b,size,b1,b2 IMPLICIT NONE INTEGER i,u REAL weight,sum,diff_sq CHARACTER*20 filename,date*9 ! PRINT *,' Give input file name or "none":' READ *,filename PRINT *,' Give output destination 6 (monitor) or 1 (file):' READ *,u PRINT *,' Give smoothing factor 0-1:' READ *,weight IF (filename /= 'none') THEN OPEN(1,file=filename,action='read') READ(1,*) !read past two header lines READ(1,*) size=1 10 READ(1,*,end=20)date,a(size) PRINT *,size,a(size) size=size+1 GO TO 10 20 CLOSE(1) size=size-1 ELSE CALL MakeData !get array a END IF CALL ExpSmooth(weight) IF (u==1) OPEN(1,file='smooth.out',action='write') sum=0. DO i=1,size diff_sq=(a(i)-b(i))**2 sum=sum+diff_sq WRITE(u,1000)i,a(i),b(i),diff_sq END DO WRITE(u,1010)sum IF (u==1) CLOSE(1) 1000 FORMAT(i3,3f10.3) 1010 FORMAT(' Sum of squares = ',f10.3) END ! Program 9.8 (WATER.F90) MODULE WaterBillDefinitions ! TYPE MeterCode REAL MinCharge INTEGER allowance END TYPE TYPE CustomerClass CHARACTER*2 ClassCode INTEGER min,max REAL Per1000Gallons END TYPE ! CONTAINS !----------------------------------------- SUBROUTINE MeterCodeData(MeterData) ! ! Read H2O_RATE.DAT and extract billing information that ! depends on meter code. ! IMPLICIT NONE TYPE (MeterCode), INTENT(OUT) :: MeterData(10) INTEGER i,code ! ! Read 4 header lines... DO i=1,4 READ(1,*) END DO ! ! Read meter code data... DO i=1,10 READ(1,*)code,MeterData(code)%MinCharge, & MeterData(code)%allowance END DO END SUBROUTINE MeterCodeData !-------------------------------------------- SUBROUTINE CustomerClassData(ClassData) ! ! Read H2O_RATE.DAT and extract rate information that depends ! on customer class. ! IMPLICIT NONE TYPE (CustomerClass) ClassData(14) INTEGER i CHARACTER a*78,b*10 ! ! Read 4 header lines... DO i=1,4 READ(1,*) END DO ! ! Read class data... DO i=1,14 READ(1,1000)a ! PRINT *,a ClassData(i)%ClassCode=a(26:27) READ(a(33:40),*)ClassData(i)%min b=a(45:52) IF (b(5:8)=='over') THEN ClassData(i)%max=HUGE(1) ELSE READ(a(45:52),*)ClassData(i)%max END IF READ(a(57:61),*)ClassData(I)%Per1000Gallons END DO ! 1000 FORMAT(a78) END SUBROUTINE CustomerClassData !------------------------------------------------------------- SUBROUTINE CalculateBill(code,meter,gallons,MeterData, & RateData,QuarterlyBill,index) ! IMPLICIT NONE CHARACTER, INTENT(IN) :: code INTEGER, INTENT(IN) :: gallons,meter TYPE (MeterCode), INTENT(IN) :: MeterData(10) TYPE (CustomerClass), INTENT(IN) :: RateData(14) REAL, INTENT(OUT) :: QuarterlyBill INTEGER, INTENT(OUT) :: index INTEGER MinimumGallons ! ! Based on meter code... QuarterlyBill=MeterData(meter)%MinCharge MinimumGallons=MeterData(meter)%allowance ! ! Based on customer class and consumption... SELECT CASE (code) CASE ('r') index=1 IF (gallons>RateData(1)%max) index=2 IF (gallons>MinimumGallons) & QuarterlyBill=QuarterlyBill+ & MIN(RateData(1)%max-MinimumGallons,gallons-MinimumGallons)/1e3* & RateData(1)%Per1000Gallons+ & MAX(gallons-RateData(1)%max,0)/1e3* & RateData(2)%Per1000Gallons CASE ('u') index=3 IF (gallons>RateData(3)%max) index=4 IF (gallons>RateData(4)%max) index=5 IF (gallons>MinimumGallons) & QuarterlyBill=QuarterlyBill+ & MIN(RateData(3)%max-MinimumGallons,gallons-MinimumGallons)/1e3* & RateData(3)%Per1000Gallons+ & MIN(RateData(4)%max-RateData(3)%max, & MAX(gallons-RateData(3)%max,0))/1e3* & RateData(4)%Per1000Gallons+ & MAX(gallons-RateData(4)%max,0)/1e3* & RateData(5)%Per1000Gallons CASE ('c') index=6 IF (gallons>RateData(6)%max) index=7 IF (gallons>RateData(7)%max) index=8 IF (gallons>RateData(8)%max) index=9 IF (gallons>MinimumGallons) & QuarterlyBill=QuarterlyBill+ & MIN(RateData(6)%max-MinimumGallons,gallons-MinimumGallons)/1e3* & RateData(6)%Per1000Gallons+ & MIN(RateData(7)%max-RateData(6)%max, & MAX(gallons-RateData(6)%max,0))/1e3* & RateData(7)%Per1000Gallons+ & MIN(RateData(8)%max-RateData(7)%max, & MAX(gallons-RateData(7)%max,0))/1e3* & RateData(8)%Per1000Gallons+ & MAX(gallons-RateData(8)%max,0)/1e3* & RateData(9)%Per1000Gallons CASE ('i') index=10 IF (gallons>RateData(10)%max) index=11 IF (gallons>RateData(11)%max) index=12 IF (gallons>RateData(12)%max) index=13 IF (gallons>RateData(13)%max) index=14 IF (gallons>MinimumGallons) & QuarterlyBill=QuarterlyBill+ & MIN(RateData(10)%max-MinimumGallons,gallons-MinimumGallons)/1e3* & RateData(10)%Per1000Gallons+ & MIN(RateData(11)%max-RateData(10)%max, & MAX(gallons-RateData(10)%max,0))/1e3* & RateData(11)%Per1000Gallons+ & MIN(RateData(12)%max-RateData(11)%max, & MAX(gallons-RateData(11)%max,0))/1e3* & RateData(12)%Per1000Gallons+ & MIN(RateData(13)%max-RateData(12)%max, & MAX(gallons-RateData(12)%max,0))/1e3* & RateData(13)%Per1000Gallons+ & MAX(gallons-RateData(13)%max,0)/1e3* & RateData(14)%Per1000Gallons END SELECT END SUBROUTINE CalculateBill ! END MODULE WaterBillDefinitions !==================================== PROGRAM WaterBill ! ! File name WATER.F90. ! Billing program for water service. ! USE WaterBillDefinitions, ONLY : MeterCode,CustomerClass, & MeterCodeData,CustomerClassData,CalculateBill IMPLICIT NONE TYPE (MeterCode) MeterData(10) TYPE (CustomerClass) RateData(14) INTEGER i,meter,gallons,RateIndex REAL QuarterlyBill CHARACTER code ! OPEN(1,file='c:\ftn90\source\h2o_rate.dat',action='read') ! Get meter code data. CALL MeterCodeData(MeterData) DO i=1,10 PRINT 1000,i,MeterData(i)%MinCharge,MeterData(i)%allowance END DO ! ! Get customer class data. CALL CustomerClassData(RateData) DO i=1,14 PRINT 1010,RateData(i)%ClassCode,RateData(i)%min, & RateData(i)%max,RateData(i)%Per1000Gallons END DO CLOSE(1) ! ! Read test data file. OPEN(1,file='c:\ftn90\source\water.dat',action='read') 10 READ(1,1050,end=900)code,meter,gallons CALL CalculateBill(code,meter,gallons,MeterData,RateData, & QuarterlyBill,RateIndex) PRINT 1060,code,meter,gallons,QuarterlyBill, & RateData(RateIndex)%ClassCode GO TO 10 900 CLOSE(1) PRINT *,'Test data processed successfully.' ! 1000 FORMAT(i3,f10.2,i6) 1010 FORMAT(a3,2i12,f8.3) 1050 FORMAT(a1,i2,i10) 1060 FORMAT(a2,i3,i10,f10.2,a3) END ! Program 9.9 (MERGE.F90) MODULE MergeData ! IMPLICIT NONE TYPE fields INTEGER x END TYPE fields ! CONTAINS !-------------------------------------------- SUBROUTINE ReadOne(unit,value,end_file) ! ! Read one value and return to program; set an e-o-f flag. ! IMPLICIT NONE INTEGER, INTENT(IN) :: unit TYPE (fields), INTENT(OUT) :: value LOGICAL, INTENT(OUT) :: end_file ! end_file=.false. READ(unit,*,END=20)value%x RETURN 20 end_file=.true. END SUBROUTINE ReadOne !----------------------------- SUBROUTINE ReadAll(unit) ! ! Read to end of file. ! IMPLICIT NONE INTEGER, INTENT(IN) :: unit TYPE (fields) value ! 10 READ(unit,*,END=20)value%x PRINT*,value%x GO TO 10 20 CONTINUE END SUBROUTINE ReadAll !----------------------------------------- SUBROUTINE MergeLists(unit_a,unit_b) ! IMPLICIT NONE INTEGER, INTENT(IN) :: unit_a,unit_b TYPE (fields) list_a,list_b LOGICAL end_a,end_b DATA end_a,end_b/2*.false./ ! CALL ReadOne(unit_a,list_a,end_a) CALL ReadOne(unit_b,list_b,end_b) 10 CONTINUE IF (list_a%x .LT. list_b%x) THEN PRINT*,list_a%x IF (.NOT. end_a) CALL ReadOne(unit_a,list_a,end_a) IF (end_b) PRINT*,list_b%x ELSE IF (list_a%x .EQ. list_b%x) THEN PRINT*,list_a%x PRINT*,list_b%x IF (.NOT. end_a) CALL ReadOne(unit_a,list_a,end_a) IF (.NOT. end_b) CALL ReadOne(unit_b,list_b,end_b) ELSE PRINT*,list_b%x IF (.NOT. end_b) CALL ReadOne(unit_b,list_b,end_b) IF (end_a) PRINT*,list_a%x END IF IF ((.NOT. end_a) .AND. (.NOT. end_b)) GO TO 10 IF (.NOT. end_a) THEN PRINT*,list_a%x CALL ReadAll(unit_a) END IF IF (.NOT. end_b) THEN PRINT*,list_b%x CALL ReadAll(unit_b) END IF END SUBROUTINE MergeLists !------------------------------ END MODULE MergeData !========================= PROGRAM merge ! ! Merge two sorted lists of integers.. ! USE MergeData, ONLY : MergeLists IMPLICIT NONE INTEGER unit_a,unit_b DATA unit_a,unit_b/1,2/ ! PRINT *,' opening list a...' OPEN(unit_a,file='c:\ftn90\source\lista.dat') PRINT *,' opening list b...' OPEN(unit_b,file='c:\ftn90\source\listb.dat') ! CALL MergeLists(unit_a,unit_b) ! CLOSE(unit_a) CLOSE(unit_b) ! END ! Program 9.10 (JAN_QCD.F90) PROGRAM January_QCD ! ! Create a quote-and-comma-delimited ASCII text file for ! use by a spreadsheet. ! IMPLICIT NONE CHARACTER*8 date REAL hi,lo ! OPEN(1,file='c:\ftn90\source\january.dat',action='read') OPEN(2,file='c:\ftn90\source\january.qcd',action='write') ! 10 READ(1,1000,end=900)date,hi,lo PRINT 1010,'"',date,'"',',',hi,',',lo WRITE(2,1010)'"',date,'"',',',hi,',',lo GO TO 10 900 CLOSE(1) CLOSE(2) 1000 FORMAT(a8,2f4.0) 1010 FORMAT(1x,a1,a8,a1,a1,f4.0,a1,f4.0) END ! Program 10.1 - 10.4 (SEARCH.F90) MODULE Setup IMPLICIT NONE CHARACTER(20) a(100),target ! CONTAINS !--------------------------------- SUBROUTINE GetList(size) ! IMPLICIT NONE INTEGER, INTENT(INOUT) :: size ! OPEN(1,file="c:\ftn90\source\search.dat") size=1 10 READ(1,*,END=900)a(size) size=size+1 GO TO 10 900 size=size-1 CLOSE(1) ! PRINT *,'# of items in list = ',size END SUBROUTINE GetList !------------------------------- END MODULE Setup !====================== MODULE SearchSubs CONTAINS ! SUBROUTINE FindFirst(size,where) ! ! Linear search of a list for first occurrence of ! specified target value. ! USE Setup, ONLY : a,target IMPLICIT NONE INTEGER, INTENT(IN) :: size INTEGER, INTENT(OUT) :: where ! where=1 DO WHILE ((a(where) /= target) .AND. (where < size)) where=where+1 END DO IF (a(where) /= target) where=0 ! END SUBROUTINE FindFirst !---------------------------------------- SUBROUTINE FindAll(size,how_many) ! ! Linear search of a list for all occurrences of ! specified target value. ! USE Setup, ONLY : a,target IMPLICIT NONE INTEGER i INTEGER, INTENT(IN) :: size INTEGER, INTENT(OUT) :: how_many ! how_many=0 DO i=1,size IF (a(i) .EQ. target) THEN PRINT *,target,' at position ',i how_many=how_many+1 END IF END DO ! END SUBROUTINE FindAll !---------------------------------------- SUBROUTINE Binary(low,high,where) ! ! Binary search of an ordered list for one occurrence of ! specified target value. ! USE Setup, ONLY : a,target IMPLICIT NONE INTEGER mid,lo,hi INTEGER, INTENT(OUT) :: where INTEGER, INTENT(IN) :: low,high ! lo=low !assign these values locally so... hi=high !INTENT(IN) won't be violated. where=0 DO WHILE ((lo .LE. hi) .AND. (where .EQ. 0)) mid=(lo+hi)/2 IF (a(mid) .EQ. target) THEN where=mid ELSE IF (a(mid) .GT. target) THEN hi=mid-1 ELSE lo=mid+1 END IF END DO ! END SUBROUTINE Binary !---------------------------- END MODULE SearchSubs !=========================== PROGRAM Search ! ! Driver program to test linear and binary search algorithms. ! USE Setup, ONLY : a,target,GetList USE SearchSubs, ONLY : FindAll,FindFirst,Binary IMPLICIT NONE CHARACTER YesNo,choice INTEGER size,where,how_many ! CALL GetList(size) !Build list to search. ! 10 CONTINUE !START LOOP PRINT *,' What name would you like to look for?' READ *,target PRINT *,'Linear search for first (f), all (a), or binary (b)?' READ *,choice SELECT CASE (choice) CASE ('a') !linear search for all occurrences CALL FindAll(size,how_many) IF (how_many .GT. 0) THEN PRINT *,how_many,' occurrences found' ELSE PRINT *,target,' not found' END IF CASE ('b') !binary search for one occurrence CALL Binary(1,size,where) IF (where .GT. 0) THEN PRINT *,target,' found at position',where ELSE PRINT *,target,' not found' END IF CASE ('f') !linear search for first occurrence CALL FindFirst(size,where) IF (where .GT. 0) THEN PRINT *,target,' found at position',where ELSE PRINT *,target,' not found' END IF END SELECT PRINT *,' Try again (y/n)?' READ *,YesNo IF (YesNo .EQ. 'y') GO TO 10 ! END ! Program 10.5 - 10.7 (SORT.F90) MODULE SortSetup IMPLICIT NONE INTEGER array(100),save_value ! (see also declarations for a, b, and temp in Swap) ! CONTAINS ! SUBROUTINE Swap(a,b) ! IMPLICIT NONE INTEGER, INTENT(INOUT) :: a,b INTEGER temp ! temp=a a=b b=temp END SUBROUTINE Swap ! SUBROUTINE GetList(size) IMPLICIT NONE INTEGER, INTENT(OUT) :: size INTEGER i INTEGER Count(1) ! for random number generator REAL x ! CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) DO i=1,size CALL Random_Number(x) array(i)=x*100+1 END DO END SUBROUTINE GetList ! SUBROUTINE PrintList(size) IMPLICIT NONE INTEGER, INTENT(IN) :: size INTEGER i ! DO i=1,size PRINT *,array(i) END DO END SUBROUTINE PrintList ! END MODULE SortSetup ! MODULE SortSubs CONTAINS ! SUBROUTINE Selection(size) ! USE SortSetup, ONLY : a => array,Swap IMPLICIT NONE INTEGER, INTENT(IN) :: size INTEGER current_cell,i,smallest ! DO current_cell=1,size-1 smallest=current_cell DO i=current_cell+1,size IF (a(i) .LT. a(smallest)) smallest=i END DO IF (a(current_cell) .NE. a(smallest)) & CALL Swap(a(current_cell),a(smallest)) END DO END SUBROUTINE Selection ! SUBROUTINE Insertion(size) ! Assumes size>=3 USE SortSetup, ONLY : a => array,save_value,Swap IMPLICIT NONE INTEGER, INTENT(IN) :: size INTEGER current,smallest,test ! ! Find smallest element... smallest=1 DO current=2,size IF (a(current) array,Swap IMPLICIT NONE INTEGER, INTENT(IN) :: lower,upper INTEGER, INTENT(OUT) :: mid INTEGER left,right ! left=lower right=upper ! ! Optionally, remove the "!" from CALL Swap to use a value from the ! middle of the list as the partitioning value. This will increase ! the efficiency of the algorithm when the list isn't in random order. ! (Because the algorithm works best when the partitioning value is ! the median value in the list.) ! CALL Swap(a(lower),a((lower+upper)/2)) DO WHILE (left .LT. right) DO WHILE (a(right) .GT. a(lower)) right=right-1 END DO DO WHILE ((left .LT. right) .AND. (a(left) .LE. a(lower))) left=left+1 END DO IF (left .LT. right) CALL Swap(a(left),a(right)) END DO mid=right CALL Swap(a(mid),a(lower)) END SUBROUTINE Partition ! RECURSIVE SUBROUTINE QuickSort(lower,upper) USE SortSetup, ONLY : a => array,Swap IMPLICIT NONE INTEGER, INTENT(IN) :: lower,upper INTEGER mid ! IF (lower .LT. upper) THEN CALL Partition(lower,upper,mid) CALL QuickSort(lower,mid-1) CALL QuickSort(mid+1,upper) END IF ! END SUBROUTINE QuickSort ! END MODULE SortSubs ! PROGRAM Sort ! ! Test and demonstrate several sorting algorithms. ! USE SortSetup, ONLY : a => array,GetList,PrintList USE SortSubs IMPLICIT NONE INTEGER size CHARACTER choose,YesNo ! size=20 ! 10 PRINT *,' Choose sorting algorithm: ' PRINT *,' Insertion (i), Quicksort (q), Selection (s)' READ *,choose CALL GetList(size) SELECT CASE (choose) CASE ('i') CALL Insertion(size) CASE ('q') CALL QuickSort(1,size) CASE ('s') CALL Selection(size) END SELECT CALL PrintList(size) PRINT *,' Again (y/n)?' READ *,YesNo IF (YesNo == 'y') GO TO 10 ! END ! Program 10.8 (FACTORAL.F90) MODULE FactorialSubs !------------------------- CONTAINS !------------------------------ INTEGER FUNCTION Fact1(n) IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER i,factorial ! factorial=1 IF (n .GT. 1) THEN DO i=2,n factorial=factorial*i END DO END IF Fact1=factorial END FUNCTION Fact1 !------------------------------ INTEGER FUNCTION Fact2(n) IMPLICIT NONE INTEGER, INTENT(IN) :: n INTEGER factorial,temp ! temp=n ! use temp to avoid violating INTENT(IN) on n factorial=1 IF (temp > 1) THEN DO WHILE (temp > 1) factorial=factorial*temp temp=temp-1 END DO END IF Fact2=factorial END FUNCTION Fact2 !----------------------------------------------------------- RECURSIVE INTEGER FUNCTION Factorial(n) RESULT(n_Fact) IMPLICIT NONE INTEGER, INTENT(IN) :: n ! IF (n <= 1) THEN n_Fact=1 ELSE n_Fact=n*Factorial(n-1) END IF END FUNCTION Factorial !----------------------------- END MODULE FactorialSubs !============================= PROGRAM n_Factorial ! ! File name FACTORAL.F90 ! Use factorial function to demonstrate recursive algorithms. ! USE FactorialSubs IMPLICIT NONE INTEGER n ! PRINT *,' Give an integer: ' READ *,n PRINT *,' Using a DO... loop, ',n,'! = ',Fact1(n) PRINT *,' Using a WHILE... loop, ',n,'! = ',Fact2(n) PRINT *,' Using a recursive function, ',n,'! = ',Factorial(n) END ! Program 10.9 (QUIKSORT.F90) MODULE Quick_Sort ! INTEGER a(100),pivot,temp CONTAINS !------------------------------------------------- SUBROUTINE Partition(lower,upper,left,right) ! IMPLICIT NONE INTEGER,INTENT(IN) :: lower,upper INTEGER, INTENT(INOUT) :: left,right ! left=lower ! Start at bottom and top of list. right=upper pivot=a((lower+upper)/2) ! Begin post-test loop. 10 CONTINUE DO WHILE (a(left)left) GO TO 10 ! End post-test loop. END SUBROUTINE Partition !------------------------------------------------ RECURSIVE SUBROUTINE QuickSort(lower,upper) ! IMPLICIT NONE INTEGER left,right INTEGER, INTENT(IN) :: lower,upper ! CALL Partition(lower,upper,left,right) IF (lower a, QuickSort IMPLICIT NONE INTEGER i,n ! PRINT *,' How many values to sort (<=100)? ' READ *,n PRINT*,'Type ',n,' integers...' READ*,(x(i),i=1,n) PRINT*,(x(i),i=1,n) ! CALL QuickSort(1,n) PRINT*,(x(i),i=1,n) END ! Program 10.10 (INSERT2.F90) MODULE DataDef IMPLICIT NONE INTEGER size,where,max_size REAL a(22),value PARAMETER (max_size=22) ! END MODULE DataDef ! MODULE InsertSubs CONTAINS ! SUBROUTINE GetPosition ! USE DataDef, ONLY : a,size,where IMPLICIT NONE ! IF (a(size) .LE. a(1)) THEN where=1 ELSE where=size-1 DO WHILE (a(where) .GT. a(size)) where=where-1 END DO where=where+1 END IF ! PRINT*,' Put ', a(size),' in position ',where END SUBROUTINE GetPosition ! SUBROUTINE Add ! ! Inserts a value in an array at a specified location. ! USE DataDef, ONLY : a,size,value,where IMPLICIT NONE ! size=size+1 a(size)=value IF (a(size) .LT. a(size-1)) THEN CALL GetPosition a(where+1:size)=a(where:size-1) a(where)=value ENDIF ! END SUBROUTINE Add ! SUBROUTINE Remove ! ! Removes an element at specified location. ! USE DataDef, ONLY : a,size,where IMPLICIT NONE ! a(where:size-1)=a(where+1:size) size=size-1 ! END SUBROUTINE Remove ! END MODULE InsertSubs ! !********************************************************************* ! PROGRAM Insert ! ! Demonstrate an insertion algorithm that keeps a list sorted. ! USE DataDef USE InsertSubs IMPLICIT NONE INTEGER i CHARACTER choice ! ! Create an array of integers in ascending order. ! size=20 DO i=1,size a(i)=2*i END DO CALL PrintList choice=' ' ! DO WHILE (choice .NE. 'q') PRINT*,'(a)dd or (r)emove value, or (q)uit?' READ*,choice SELECT CASE (choice) CASE ('a') IF (size .LT. max_size) THEN PRINT*,' Give a number to add to this list:' READ*,value CALL Add CALL PrintList ELSE PRINT*,"Can't add a number to this list." END IF CASE ('r') PRINT*,' Give position of a number to remove from this list'& ,' between 1 and ',size READ*,where IF ((where .GE. 1) .and. (where .LE. size)) THEN CALL Remove CALL PrintList ELSE PRINT*,'Position is out of range.' END IF CASE ('q') CASE DEFAULT PRINT*,'INPUT ERROR - THIS CHOICE NOT AVAILABLE' END SELECT END DO ! CONTAINS ! internal subroutine SUBROUTINE PrintList IMPLICIT NONE INTEGER j DO j=1,size WRITE(*,1000,advance='no')a(j) END DO PRINT * 1000 FORMAT(f5.1) END SUBROUTINE PrintList ! END PROGRAM ! Program 10.11 (LEGENDRE.F90) MODULE Poly ! CONTAINS ! RECURSIVE REAL FUNCTION Legendre(n,x) RESULT(LegendrePoly) IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x ! IF (n==0) THEN LegendrePoly=1. ELSE IF (n==1) THEN LegendrePoly=x ELSE LegendrePoly=(2.*REAL(n)-1.)/REAL(n)*x*Legendre(n-1,x)- & (REAL(n)-1.)/REAL(n)*Legendre(n-2,x) END IF END FUNCTION Legendre !------------------------------ END MODULE Poly !================================ PROGRAM LegendreTest ! ! MS-DOS file LEGENDRE.F90 ! Calculate Legendre polynomials using a recursive relationship. ! USE Poly, ONLY : Legendre IMPLICIT NONE INTEGER i REAL x x=.5 PRINT *,'n=0-7 from recursive function, x=0.5' DO i=0,7 PRINT *,Legendre(i,x) END DO ! Check against the first 8 functions (from tables) PRINT *,'n=0-7 from tabulated polynomials, x=0.5' PRINT *,0. PRINT *,1. PRINT *,(3.*x**2-1.)/2. PRINT *,(5.*x**3-3.*x)/2. PRINT *,(35.*x**4-30*x**2+3.)/8. PRINT *,(63.*x**5-70.*x**3+15.*x)/8. PRINT *,(231.*x**6-315.*x**4+105.*x**2-5.)/16. PRINT *,(429.*x**7-693.*x**5+315.*x**3-35.*x)/16. ! END ! Program 11.1 (STATS.F90) MODULE DescriptiveStats ! CONTAINS !-------------------------------------------------- SUBROUTINE NormalStats(a,n,flag,avg,std_dev) ! ! Basic descriptive statistics for normally distributed data. ! Calculates either sample or population statistics. ! Sets std_dev = -1 if error condition is detected. ! IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: a(n) CHARACTER, INTENT(IN) :: flag ! 'P' or 'S' statistics REAL, INTENT(OUT) :: avg,std_dev REAL sum_,sum_sq,variance ! sum_=SUM(a) sum_sq=DOT_PRODUCT(a,a) SELECT CASE (flag) CASE ('p','P') variance=(sum_sq-sum_**2/n)/n CASE ('s','S') variance=(sum_sq-sum_**2/n)/(n-1) CASE DEFAULT PRINT *," FROM NormalStats: FLAG ERROR, 'P' assumed" variance=(sum_sq-sum_**2/n)/n END SELECT IF (variance < 0.) THEN !an error condition exists PRINT *,' FROM NormalStats: NEGATIVE VARIANCE', variance std_dev=-1 ELSE std_dev=SQRT(variance) END IF avg=sum_/n ! END SUBROUTINE NormalStats !--------------------------------------------------------- SUBROUTINE LinearRegression(x,y,n,flag,a,b,s_yx,r) ! ! For data to be represented by y=a+bx, calculates linear ! regression coefficients, sample standard error of y on x, and ! sample correlation coefficient. Sets r=0 if error condition ! exists. If the intercept coefficient a is set to 0 on input, ! the regression is forced through (0,0). ! IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n),y(n) CHARACTER, INTENT(IN) :: flag !'P' or 'S' statistics REAL, INTENT(OUT) :: b,s_yx,r REAL, INTENT(INOUT) :: a REAL avg,std_dev REAL sum_x,sum_y,sum_xy,sum_xx,sum_yy,temp ! sum_x=SUM(x) sum_y=SUM(y) sum_xy=DOT_PRODUCT(x,y) sum_xx=DOT_PRODUCT(x,x) sum_yy=DOT_PRODUCT(y,y) IF (a /= 0.) THEN !calculate full regression temp=n*sum_xx-sum_x**2 a=(sum_y*sum_xx-sum_x*sum_xy)/temp b=(n*sum_xy-sum_x*sum_y)/temp s_yx=SQRT((sum_yy-a*sum_y-b*sum_xy)/n) ELSE !just calculate slope b=sum_y/sum_x s_yx=SQRT((sum_yy-2.*b*sum_xy+b*b*sum_xx)/n) END IF SELECT CASE (flag) CASE ('p','P') CASE ('s','S') s_yx=s_yx*SQRT(REAL(n)/REAL(n-2)) CASE DEFAULT PRINT *," FROM LinearRegression: FLAG ERROR, 'P' assumed" END SELECT ! Use NormalStats to get standard deviation of y. CALL NormalStats(y,n,flag,avg,std_dev) IF (std_dev > 0.) THEN temp=1.-(s_yx/std_dev)**2 IF (temp >= 0.) THEN r=SQRT(temp) ELSE ! an error condition exists r=0. PRINT *,'FROM LinearRegression: ERROR CONDITION',temp END IF ELSE ! an error condition exists r=0. END IF END SUBROUTINE LinearRegression !------------------------------------ SUBROUTINE NormalArray(a,n) ! ! Generates an array of normal random numbers from ! pairs of uniform random numbers in range 0<=x<1. ! IMPLICIT NONE INTEGER, INTENT(IN) :: n REAL, INTENT(OUT) :: a(n) INTEGER i INTEGER Count(1) ! for random number generator REAL pi,u1,u2 PARAMETER (pi=3.1415927) ! CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) CALL Random_Number(a) !fills array with uniform random DO i=1,n,2 u1=a(i) u2=a(i+1) IF (u1 == 0.) u1=1e-15 ! u1 must not be 0 a(i) =SQRT(-2.0*LOG(u1))*COS(2.0*pi*u2) a(i+1)=SQRT(-2.0*LOG(u1))*SIN(2.0*pi*u2) END DO IF (MOD(n,2) /= 0) THEN !there's one extra element IF (a(n) == 0.) a(n)=1e-15 a(n)=SQRT(-2.0*LOG(a(n)))*SIN(2.0*pi*a(n)) END IF ! END SUBROUTINE NormalArray !------------------------------- END MODULE !====================== PROGRAM Stats ! USE DescriptiveStats IMPLICIT NONE REAL x(500),y(500),avg,std_dev REAL a,b ! intercept and slope for linear regression REAL s_yx ! standard error of estimate of y on x REAL corr ! correlation coefficient INTEGER n ! # of points in array INTEGER i DATA n/100/ ! ! Test basic statistics... CALL NormalArray(x,n) CALL NormalStats(x,n,'p',avg,std_dev) PRINT *,' population mean and std dev: ',avg,std_dev CALL NormalStats(x,n,'s',avg,std_dev) PRINT *,' sample mean and std dev: ',avg,std_dev ! Test linear regression... CALL NormalArray(y,n) ! Create a linear relationship with "noise"... DO i=1,n x(i)=i y(i)=2.*i+10.*y(i) END DO ! Set a /= 0 for full regression analysis. a=1. CALL LinearRegression(x,y,n,'s',a,b,s_yx,corr) PRINT *,' FOR FULL REGRESSION...' PRINT *,' regression coefficients: ',a,b PRINT *,' standard error of estimate of y on x : ',s_yx PRINT *,' correlation coefficient: ',corr ! Set a=0 for regression forced through (0,0) a=0. CALL LinearRegression(x,y,n,'s',a,b,s_yx,corr) PRINT *,' FOR REGRESSION FORCED THROUGH (0,0)...' PRINT *,' regression coefficients: ',a,b PRINT *,' standard error of estimate of y on x : ',s_yx PRINT *,' correlation coefficient: ',corr ! END ! Program 11.2 (FALLING.F90) MODULE Numerical_Differentiation ! CONTAINS !-------------------------------------------------- REAL FUNCTION Stirling(x1,x2,x3,y1,y2,y3) ! IMPLICIT NONE REAL, INTENT(IN) :: x1,x2,x3,y1,y2,y3 ! Stirling=((y2-y1)/(x2-x1)+(y3-y2)/(x3-x2))/2. END FUNCTION Stirling !----------------------------------------- END MODULE Numerical_DIfferentiation !========================================= PROGRAM Falling ! ! Driver for numerical differentiation routines. ! USE Numerical_Differentiation, ONLY : Stirling IMPLICIT NONE TYPE fall_data REAL true_time REAL true_distance REAL measured_time END TYPE fall_data TYPE(fall_data) fall(0:20) REAL g,x,true_speed,speed INTEGER i,n,Count(1) PARAMETER (g=9.8) !m/s**2 (gravitational acceleration) ! OPEN(1,file='c:\ftn90\source\falling.dat') ! ! Get data... ! (Read past 1 header line.) READ(1,*) CALL System_Clock(Count(1)) CALL Random_Seed(Put=Count) n=-1 10 n=n+1 READ(1,*,end=900)fall(n)%true_time,fall(n)%true_distance CALL Random_Number(x) fall(n)%measured_time=fall(n)%true_time+x*.4-.2 PRINT 1000,fall(n)%true_time,fall(n)%true_distance, & fall(n)%measured_time GO TO 10 900 CLOSE(1) n=n-1 ! Calculate numerical derivative... PRINT *,' true meas. true true' PRINT *,' time time distance speed speed speed/(true speed)' DO i=1,n-1 speed=Stirling(fall(i-1)%measured_time, & fall(i )%measured_time, & fall(i+1)%measured_time, & fall(i-1)%true_distance, & fall(i )%true_distance,& fall(i+1)%true_distance) true_speed=Stirling(fall(i-1)%true_time, & fall(i )%true_time, & fall(i+1)%true_time, & fall(i-1)%true_distance, & fall(i )%true_distance,& fall(i+1)%true_distance) PRINT 1010,fall(i)%true_time,fall(i)%measured_time, & fall(i)%true_distance,true_speed,speed, & speed/true_speed END DO ! 1000 FORMAT(1x,f6.2,2f10.3) 1010 FORMAT(1x,2f6.2,3f10.3,f10.2) END ! Program 11.3 (NUM_INT.F90) MODULE NumericalIntegration ! CONTAINS !------------------------------------------------ REAL FUNCTION SimpsonsRule(F,x1,x2,n_steps) ! ! Does numerical integration using Simpson's Rule. ! IMPLICIT NONE REAL, INTENT(IN) :: x1,x2 REAL sum_odd,sum_even INTEGER, INTENT(IN) :: n_steps REAL dx INTEGER i !------------------------------ INTERFACE REAL FUNCTION F(x) REAL, INTENT(IN) :: x END FUNCTION F END INTERFACE !------------------------------ dx=(x2-x1)/n_steps sum_odd=0. DO i=1,n_steps-1,2 sum_odd=sum_odd+F(x1+REAL(i)*dx) END DO sum_even=0. DO i=2,n_steps-2,2 sum_even=sum_even+F(x1+REAL(i)*dx) END DO SimpsonsRule=(F(x1)+F(x2)+4.*sum_odd+2.*sum_even)*dx/3. ! END FUNCTION SimpsonsRule !------------------------------ END MODULE NumericalIntegration !==================================== MODULE FunctionDefinition ! CONTAINS !-------------------------- REAL FUNCTION FofX(x) ! ! This is the cumulative probability function for the ! normal function. Its integral doesn't have an analytic form. ! IMPLICIT NONE REAL, INTENT(IN) :: x REAL pi pi=4.*ATAN(1.) FofX=EXP(-x*x/2.0)/SQRT(2.0*pi) RETURN END FUNCTION FofX !---------------------------------- END MODULE FunctionDefinition !================================== PROGRAM Num_Int ! ! Performs numerical integration on specified function. ! USE NumericalIntegration, ONLY : SimpsonsRule USE FunctionDefinition, ONLY : FofX IMPLICIT NONE REAL z PRINT *,' Integrate normal function to z = ' PRINT *,' (Appropriate values in range -10 to +10.)' READ *,z IF (z /= 0.0) THEN PRINT *," With Simpson's Rule..." PRINT 1000,SimpsonsRule(FofX,0.,z,100)+.5 ELSE PRINT 1000,0.5 END IF ! 1000 FORMAT(1x,f8.5) END ! Program 11.4 (GAMMA.F90) MODULE Gamma_Calc ! CONTAINS ! REAL FUNCTION Gamma(n) IMPLICIT NONE REAL x,n,z,f,dx,x_max,sum INTEGER i,n_terms f(x,z)=EXP(-x)*x**(z-1.0) DATA dx,x_max/0.001,10.0/ ! n_terms=NINT(x_max/dx) PRINT*,n_terms sum=0.0 DO i=2,n_terms-2,2 sum=sum+2.0*f(i*dx,n) END DO DO i=1,n_terms-1,2 sum=sum+4.0*f(i*dx,n) END DO sum=(sum+f(x_max,n))*dx/3.0 Gamma=sum ! END FUNCTION Gamma !----------------------- END MODULE !=========================== PROGRAM Gamma_Function ! ! Evaluate the Gamma function. ! USE Gamma_Calc, ONLY : Gamma IMPLICIT NONE REAL n,x ! PRINT*,'Give 0 ABS(a(PivotRow,row)) ) PivotRow=i END DO ! ! Swap pivot row if required. IF (PivotRow /= row) THEN PRINT*,'swapping pivot row...' DO col=row,n_cols temp=a(PivotRow,col) a(PivotRow,col)=a(row,col) a(row,col)=temp END DO PRINT*,'swap done...' CALL PrintMatrix(a,n_rows,n_cols) END IF END IF !IF (row < n_rows) ... ! ! Divide by pivot. DivideBy=a(row,row) a(row,row:n_cols)=a(row,row:n_cols)/DivideBy ! ! Reduce the (row)th column to 0 IF (row < n_rows) THEN DO i=row+1,n_rows a(i,row+1:n_cols) = a(i,row+1:n_cols)-a(row,row+1:n_cols)*a(i,row) a(i,row)=0. END DO END IF ! ! Print reduced matrix. CALL PrintMatrix(a,n_rows,n_cols) END DO ! ! Backsolve for solutions. solutions(n_rows)=a(n_rows,n_cols) DO row=n_rows-1,1,-1 solutions(row)=a(row,n_cols) DO i=row+1,n_rows solutions(row)=solutions(row)-a(row,i)*solutions(i) END DO END DO ! RETURN END SUBROUTINE GaussianElimination ! SUBROUTINE PrintMatrix(a,n_rows,n_cols) ! IMPLICIT NONE INTEGER, INTENT(IN) :: n_rows,n_cols REAL, INTENT(IN) :: a(n_rows,n_cols) INTEGER rows,cols ! PRINT*,'from PrintMatrix...' DO rows=1,n_rows DO cols=1,n_cols WRITE(*,1000,advance='no')a(rows,cols) END DO WRITE(*,1010) END DO ! 1000 FORMAT(f7.2) 1010 FORMAT() ! RETURN END SUBROUTINE PrintMatrix ! END MODULE !-------------------------------------------------------- PROGRAM Gauss ! ! Driver program for LinearSystemSubs, ! including Gaussian elimination with partial pivoting. ! USE LinearSystemSubs IMPLICIT NONE REAL, DIMENSION(:,:), ALLOCATABLE :: a REAL, DIMENSION(:), ALLOCATABLE :: solutions INTEGER n_rows,n_cols,rows,cols ! OPEN(1,file='c:\ftn90\source\gauss.dat') ! READ(1,*)n_rows ! number of equations n_cols=n_rows+1 ALLOCATE(a(n_rows,n_cols),solutions(n_rows)) ! DO rows=1,n_rows READ(1,*)(a(rows,cols),cols=1,n_cols) PRINT*,(a(rows,cols),cols=1,n_cols) END DO ! ! Print original matrix. CALL PrintMatrix(a,n_rows,n_cols) ! ! Solve system. CALL GaussianElimination(a,n_rows,n_cols,solutions) WRITE(*,1020)(solutions(rows),rows=1,n_rows) ! 1000 FORMAT(f7.2) 1010 FORMAT() 1020 FORMAT(e15.5) END ! Program 11.6 (ROOTS.F90) MODULE RootSubs ! TYPE RootType REAL root,f_root,interval END TYPE RootType ! CONTAINS !---------------------------------------------------- SUBROUTINE Bisection(F,XL,XR,n_intervals,roots) ! ! Set up intervals for finding roots with bisection method. ! IMPLICIT NONE REAL, INTENT(IN) :: XL,XR INTEGER, INTENT(IN) :: n_intervals TYPE (RootType), INTENT(OUT) :: roots(n_intervals) REAL dx,xa,xb INTEGER i !------------------------------- INTERFACE REAL FUNCTION F(x) REAL, INTENT(IN) :: x END FUNCTION F END INTERFACE !------------------------------ ! dx=(XR-XL)/n_intervals DO i=1,n_intervals xa=XL+REAL(i-1)*dx xb=XL+REAL(i)*dx IF (F(xa)*F(xb) < 0.) THEN CALL Bisect(F,xa,xb,roots(i)%root,roots(i)%f_root, & roots(i)%interval) ELSE roots(i)%root=0.; roots(i)%f_root=0.; roots(i)%interval=0. END IF END DO ! END SUBROUTINE Bisection !---------------------------------------------------------- SUBROUTINE Bisect(F,xa,xb,x_mid,f_mid,final_interval) ! IMPLICIT NONE REAL, INTENT(INOUT) :: xa,xb REAL, INTENT(OUT) :: x_mid,f_mid,final_interval LOGICAL hit REAL, PARAMETER :: epsilon_x=1e-5,epsilon_f=1e-5 !------------------------------- INTERFACE REAL FUNCTION F(x) REAL, INTENT(IN) :: x END FUNCTION F END INTERFACE !------------------------------ x_mid=(xb+xa)/2. hit=.false. f_mid=F(x_mid) DO WHILE (((xb-xa) > epsilon_x).AND.(ABS(f_mid) > epsilon_f).AND. & (.NOT. hit)) IF (f_mid == 0.) THEN hit=.true. ELSE IF (F(xa)*f_mid < 0) THEN xb=x_mid ELSE IF (F(xb)*f_mid < 0) THEN xa=x_mid ELSE PRINT *,' Unexplained error!' END IF x_mid=(xb+xa)/2. f_mid=F(x_mid) final_interval=xb-xa END DO ! END SUBROUTINE Bisect !-------------------------- END MODULE RootSubs !============================== MODULE FunctionDefinition ! CONTAINS !----------------------- REAL FUNCTION F(x) ! IMPLICIT NONE REAL x ! F=5.*x**3-2.*x**2+3. END FUNCTION F !---------------------------------- END MODULE FunctionDefinition !================================== PROGRAM GetRoots ! ! MS-DOS file name ROOTS.F90 ! USE RootSubs, ONLY : RootType,Bisection USE FunctionDefinition, ONLY : F IMPLICIT NONE INTEGER, PARAMETER :: n_intervals=10 TYPE (RootType) roots(n_intervals) INTEGER i ! CALL Bisection(F,-10.,0.,n_intervals,roots) DO i=1,n_intervals PRINT *,i,roots(i)%root,roots(i)%f_root,roots(i)%interval END DO ! END PROGRAM ! Program 11.7 (CIRCUIT.F90) MODULE RungeKutta ! CONTAINS !------------------------------------ REAL FUNCTION AofT(x,v,D,K,M,F) ! Calculate acceleration for "mass and spring" problem. IMPLICIT NONE REAL, INTENT(IN) :: x,v,D,K,M,F AofT=-(-F+D*v+K*x)/M END FUNCTION AofT !-------------------------- SUBROUTINE MassAndSpring(x,v,D,K,M,F,dt) ! Calculate motion for mass and spring problem with constant force term. IMPLICIT NONE REAL,INTENT(INOUT) :: x,v REAL, INTENT(IN) :: D,K,M,F,dt REAL k1_x,k2_x,k3_x,k4_x,k1_v,k2_v,k3_v,k4_v ! ! Runge-Kutta coefficients... k1_x=v k2_x=v+AofT(x,v,D,K,M,F)*dt/2. k3_x=v+AofT(x,v,D,K,M,F)*dt/2. k4_x=v+AofT(x,v,D,K,M,F)*dt k1_v=AofT(x,v,D,K,M,F) k2_v=AofT(x+v*dt/2.,v+k1_v*dt/2.,D,K,M,F) k3_v=AofT(x+v*dt/2.,v+k2_v*dt/2.,D,K,M,F) k4_v=AofT(x+v*dt,v+k3_v*dt,D,K,M,F) ! Propagate solution... x=x+(k1_x+2.*k2_x+2.*k3_x+k4_x)*dt/6. v=v+(k1_v+2.*k2_v+2.*k3_v+k4_v)*dt/6. END SUBROUTINE MassAndSpring !--------------------------------- END MODULE RungeKutta !========================== PROGRAM Circuit ! ! Use Runge-Kutta method to solve LRC circuit problems. ! Variable equivalences with mass-and-spring problem: ! V => force F ! q => displacement x ! i => velocity v ! L => mass m ! R => damping constant D ! 1/C => spring constant K USE RungeKutta, ONLY : MassAndSpring IMPLICIT NONE REAL i,q,L,C,R,V,t,dt,t_final REAL k1_i,k2_i,k3_i,k4_i INTEGER j,n CHARACTER*1 choice ! ! Choose circuit type... ! PRINT *,' Specify [o]scillating or [n]o oscillating term...' READ *,choice SELECT CASE (choice) CASE ('o','O') ! Ld^2q/dt^2+Rdq/dt+q/C=V 10 PRINT *,' Give V, L, R, C (4L/C-R^2) > 0:' READ *,V,L,R,C PRINT *,' one period at t=',4*3.14159*L/SQRT(4.*L/C-R*R),' s' PRINT *,' Give t_final and number of points:' READ *,t_final,n IF ((4.*L/C-R*R) <= 0.) GO TO 10 q=0.; i=C*V*R/2./L; t=0. !Initial values dt=t_final/REAL(n) PRINT *,' time q i analytic q' DO j=1,n CALL MassAndSpring(q,i,R,1./C,L,V,dt) t=t+dt PRINT 1000,t,q,i,C*V*(1.-EXP(-R*t/2./L)*COS(SQRT(4.*L/C-R*R)*t/2./L)) END DO 1000 FORMAT(1x,es12.6,3es12.4) ! CASE ('n','N') ! Ldi/dt+Ri=V (no oscillating term)... PRINT *,' Give V, L, R:' READ *,V,L,R PRINT *,' time constant at t=',L/R,' s' PRINT *,' Give t_final and number of points:' READ *,t_final,n i=0.; t=0. !Initial values dt=t_final/REAL(n) PRINT *,' time i analytic 1' DO j=1,n ! Runge-Kutta coefficients... k1_i=(V-R*i )/L k2_i=(V-R*(i+k1_i*dt/2.))/L k3_i=(V-R*(i+k2_i*dt/2.))/L k4_i=(V-R*(i+k3_i*dt) )/L ! Propagate solution... i=i+(k1_i+2.*k2_i+2.*k3_i+k4_i)*dt/6. t=t+dt PRINT 1000,t,i,V/R*(1.-EXP(-R*t/L)) END DO ! CASE DEFAULT PRINT *,' No such choice. Try again...' END SELECT ! END ! Program 12.1 (DERIVATV.F90) MODULE Numerical_Differentiation ! CONTAINS !-------------------------------------------------- REAL FUNCTION derivative_2(x1,x2,x3,y1,y2,y3) ! IMPLICIT NONE REAL, INTENT(IN) :: x1,x2,x3,y1,y2,y3 ! derivative_2=((y2-y1)/(x2-x1)+(y3-y2)/(x3-x2))/2. RETURN END FUNCTION derivative_2 !----------------------------------------- END MODULE Numerical_DIfferentiation ! PROGRAM Derivative ! ! Driver for numerical differentiation routines. ! File name DERIVATV.F90. ! USE Numerical_Differentiation IMPLICIT NONE TYPE distance_data REAL true_time REAL true_distance REAL true_speed REAL measured_time END TYPE distance_data TYPE(distance_data) distance(0:20) REAL g INTEGER i PARAMETER (g=9.8) !m/s**2 (gravitational acceleration) ! OPEN(1,file='derivatv.dat') ! ! Get data... ! (Read past 2 header lines.) READ(1,*) READ(1,*) DO i=0,20 READ(1,*)distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time PRINT *,distance(i)%true_time,distance(i)%true_distance END DO CLOSE(1) ! Calculate numerical derivative... DO i=0,20 IF ((i==0) .OR. (i==20)) THEN PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time ELSE PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time,& derivative_2(distance(i-1)%true_time, & distance(i) %true_time, & distance(i+1)%true_time, & distance(i-1)%true_distance, & distance(i) %true_distance, & distance(i+1)%true_distance) END IF END DO ! 1000 FORMAT(1x,f4.1,5f10.3) END ! Program 12.1(a) (DERIVAT2.F90) INCLUDE 'c:\ftn90\source\deri_mod.f90' PROGRAM Derivative ! ! Driver for numerical differentiation routines. ! File name DERIVAT2.F90. ! USE Numerical_Differentiation IMPLICIT NONE TYPE distance_data REAL true_time REAL true_distance REAL true_speed REAL measured_time END TYPE distance_data TYPE(distance_data) distance(0:20) REAL g INTEGER i PARAMETER (g=9.8) !m/s**2 (gravitational acceleration) ! OPEN(1,file='derivatv.dat') ! ! Get data... ! (Read past 2 header lines.) READ(1,*) READ(1,*) DO i=0,20 READ(1,*)distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time PRINT *,distance(i)%true_time,distance(i)%true_distance END DO CLOSE(1) ! Calculate numerical derivative... DO i=0,20 IF ((i==0) .OR. (i==20)) THEN PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time ELSE PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time,& derivative_2(distance(i-1)%true_time, & distance(i) %true_time, & distance(i+1)%true_time, & distance(i-1)%true_distance, & distance(i) %true_distance, & distance(i+1)%true_distance) END IF END DO ! 1000 FORMAT(1x,f4.1,5f10.3) END ! see DERIVAT2.F90 (DERI_MOD.F90) MODULE Numerical_Differentiation ! CONTAINS ! REAL FUNCTION derivative_2(x1,x2,x3,y1,y2,y3) ! IMPLICIT NONE REAL, INTENT(IN) :: x1,x2,x3,y1,y2,y3 ! derivative_2=((y2-y1)/(x2-x1)+(y3-y2)/(x3-x2))/2. RETURN END FUNCTION derivative_2 ! END MODULE Numerical_DIfferentiation ! ! see Chapter Twelve (DERIVAT3.F90) PROGRAM Derivative ! ! Driver for numerical differentiation routines. ! File name DERIVAT3.F90. ! USE Numerical_Differentiation !separately compiled module IMPLICIT NONE TYPE distance_data REAL true_time REAL true_distance REAL true_speed REAL measured_time END TYPE distance_data TYPE(distance_data) distance(0:20) REAL g INTEGER i PARAMETER (g=9.8) !m/s**2 (gravitational acceleration) ! OPEN(1,file='derivatv.dat') ! ! Get data... ! (Read past 2 header lines.) READ(1,*) READ(1,*) DO i=0,20 READ(1,*)distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time PRINT *,distance(i)%true_time,distance(i)%true_distance END DO CLOSE(1) ! Calculate numerical derivative... DO i=0,20 IF ((i==0) .OR. (i==20)) THEN PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time ELSE PRINT 1000,distance(i)%true_time,distance(i)%true_distance, & distance(i)%true_speed,distance(i)%measured_time,& derivative_2(distance(i-1)%true_time, & distance(i) %true_time, & distance(i+1)%true_time, & distance(i-1)%true_distance, & distance(i) %true_distance, & distance(i+1)%true_distance) END IF END DO ! 1000 FORMAT(1x,f4.1,5f10.3) END ! Program 12.2 (JULIAN_T.F90) PROGRAM Julian_t ! ! Tests a function to calculate Julian time to the nearest ! second for a specified GMT on a specified calendar day. ! IMPLICIT NONE INTEGER mon,day,year,hh,mm,ss INTEGER Calendar_to_Julian !function to get Julian date REAL DayFraction INTEGER, PARAMETER :: long=Selected_Real_Kind(14) REAL(kind=long) JulianTime ! PRINT *,' Give a calendar date in the format mm dd yyyy:' READ *,mon,day,year JulianTime=REAL(Calendar_to_Julian(mon,day,year)) PRINT 1000,JulianTime PRINT *,' Give clock time as hh mm ss, GMT:' READ *,hh,mm,ss JulianTime= JulianTime + DayFraction(hh-12,mm,ss) PRINT 1010,JulianTime 1000 FORMAT(' The Julian date at 12h00m00s GMT is ',f10.1) 1010 FORMAT(' The Julian time is ',f16.7) ! END ! INTEGER FUNCTION Calendar_to_Julian(mon,day,year) ! ! Converts Gregorian calendar date to a Julian date. ! Test value: 01/01/1995 = 2449719 (Greenwich noon). ! IMPLICIT NONE INTEGER, INTENT(IN) :: mon,day,year INTEGER temp ! temp=(mon-14)/12 Calendar_to_Julian=day-32075+1461*(year+4800+temp)/4+ & 367*(mon-2-temp*12)/12- & 3*((year+4900+temp)/100)/4 RETURN END FUNCTION Calendar_to_Julian ! REAL FUNCTION DayFraction(hh,mm,ss) ! ! Converts hours, minutes, and seconds into a fraction of a day. ! IMPLICIT NONE INTEGER, INTENT(IN) :: hh,mm,ss ! DayFraction=REAL(hh)/24.+REAL(mm)/1440.+REAL(ss)/86400. RETURN END FUNCTION DayFraction ! Program 12.3 (ARAYFUNC.F90) PROGRAM ArrayFunctions ! ! MS-DOS file name ARAYFUNC.F90. ! Demonstrate array manipulation functions. ! IMPLICIT NONE REAL a(3),b(3) REAL x(3,4),y(4,3),z(3,3) INTEGER row,col ! Here's one way to initialize an array. DATA x/1.1,2.1,3.1,1.2,2.2,3.2,1.3,2.3,3.3,1.4,2.4,3.4/ DATA y/10.1,20.1,30.1,40.1,10.2,20.2,30.2,40.2,10.3,20.3,30.3,40.3/ ! a=(/1.,2.,3./) b=(/2.,3.,4./) PRINT *,'a=',a PRINT *,'b=',b PRINT *,'DOT_PRODUCT(a,b)=',DOT_PRODUCT(a,b) ! Here's an alternate way to assign values to an array. Don't use both! x=RESHAPE((/1.1,2.1,3.1,1.2,2.2,3.2,1.3,2.3,3.3,1.4,2.4,3.4/),(/3,4/)) PRINT *,'LBOUND(x)=',LBOUND(x),'UBOUND(x)=',UBOUND(x), & 'SHAPE(x)=',SHAPE(x),'SIZE(x)=',SIZE(x) y=RESHAPE( (/10.1,20.1,30.1,40.1,10.2,20.2,30.2,40.2, & 10.3,20.3,30.3,40.3/), (/4,3/) ) PRINT *,'LBOUND(y,1)=',LBOUND(y,1),'LBOUND(y,2)=',LBOUND(y,2), & 'SHAPE(y)=',SHAPE(y),'SIZE(y,1)=',SIZE(y,1), & 'SIZE(y,2)=',SIZE(y,2) ! Using the "mask" parameter. PRINT *,'ALL(mask=(x>30.0))',ALL(mask=(x>30.0)) PRINT *,'ANY(mask=(x> 3.0))',ANY(mask=(x> 3.0)) PRINT *,'COUNT(mask=(x>3.0))',COUNT(mask=(x>3.0)) ! Finding max and min values. PRINT *,'MAXVAL(x)=',MAXVAL(x),' MINVAL(x)=',MINVAL(x) ! Finding product and sum of array elements. PRINT *,'Product across columns of x for each row.' PRINT *,'PRODUCT(x,2)',PRODUCT(x,2) PRINT *,'Product across columns of x for each row for which an element' PRINT *,'is greater than 3.' PRINT *,'PRODUCT(x,2,mask=(x>3.0))',PRODUCT(x,2,mask=(x>3.0)) PRINT *,'Sum of all elements in y.' PRINT *,'SUM(y)=',SUM(y) z=MATMUL(x,y) PRINT *,'x' PRINT 1000,((x(row,col),col=1,4),row=1,3) PRINT *,'y' PRINT 1010,((y(row,col),col=1,3),row=1,4) PRINT *,'z' PRINT 1010,((z(row,col),col=1,3),row=1,3) ! You can apply DOT_PRODUCT to one dimension of a vector. PRINT *,DOT_PRODUCT(x(1,1:4),y(1:4,1)) ! 1000 FORMAT(1x,4f8.2) 1010 FORMAT(1x,3f8.2) END ! Program 12.4 (TEM_PHIL.F90) PROGRAM temp_phil ! ! file name TEM_PHIL.F90 ! Determine average annual and monthly temperatures, using ! array reduction functions. ! IMPLICIT NONE INTEGER year,month,yr REAL temp(60:90,12),monthly_avg(12),yearly_avg(60:90) ! OPEN(1,file='c:\ftn90\source\tem_phil.dat',action='read') ! READ(1,*) DO yr=60,90 READ(1,*)year,(temp(year,month),month=1,12) END DO 900 CLOSE(1) ! ! Get monthly and yearly averages. monthly_avg=SUM(temp,1)/31. yearly_avg=SUM(temp,2)/12. PRINT *,' Monthly averages...' PRINT 1000,monthly_avg PRINT *,' Yearly averages...' PRINT 1000,yearly_avg ! 1000 FORMAT(1x,6f10.1) END ! Program 12.5 (FILETEST.F90) PROGRAM FileTest ! ! Test syntax of various file options. ! IMPLICIT NONE INTEGER i,n,a,b,c,n_bytes ! INQUIRE(file='c:\ftn90\source\filetest.dat',recl=n_bytes) PRINT *,' Record length: ',n_bytes OPEN(1,file='c:\ftn90\source\filetest.dat',action='read') INQUIRE(iolength=n_bytes)a,b,c OPEN(2,file='c:\ftn90\source\filetest.bin',form='unformatted', & action='readwrite',access='direct',recl=n_bytes) ! Read text file and use it to create direct access binary file. n=0 PRINT *,' original file...' 10 READ(1,*,end=900)a,b,c PRINT *,a,b,c n=n+1 WRITE(2,rec=n)a,b,c GO TO 10 900 CLOSE(1) ! Read binary file backwards. PRINT *,' binary file...' DO i=n,1,-1 READ(2,rec=i)a,b,c PRINT *,a,b,c END DO CLOSE(2) ! ! Now create a formatted file for direct access. OPEN(1,file='c:\ftn90\source\filetest.dat',action='read') OPEN(2,file='c:\ftn90\source\filetest.fix',action='write', & status='replace') n=0 PRINT *,' Create fixed-record length file.' 20 READ(1,*,end=910)a,b,c PRINT *,a,b,c n=n+1 WRITE(2,1000)a,b,c !Do NOT use list-directed output. GO TO 20 910 CLOSE(1) CLOSE(2) ! Open this formatted file for random access. PRINT *,' Open for random access.' ! Find the actual record length. INQUIRE(file='c:\ftn90\source\filetest.fix',recl=n_bytes) PRINT *,' Actual record length = ',n_bytes ! Open file with n_bytes+2 to allow for CR/LF characters. OPEN(1,file='c:\ftn90\source\filetest.fix',action='read', & recl=n_bytes+2,access='direct',form='formatted') DO i=1,n READ(1,1000,rec=i)a,b,c !List-directed input not allowed. PRINT *,a,b,c END DO CLOSE(1) 1000 FORMAT(3i4) END ! Program 12.5 (SRCHFILE.F90) MODULE BinSearchSub CONTAINS !------------------------------------------------ SUBROUTINE Bin_Search(u,n_rec,target,where) ! ! Binary search of an ordered list for one occurrence of ! specified target value. ! IMPLICIT NONE INTEGER mid,lo,hi INTEGER x,where INTEGER, INTENT(IN) :: u,n_rec,target ! lo=1 hi=n_rec where=0 DO WHILE ((lo .LE. hi) .AND. (where .EQ. 0)) mid=(lo+hi)/2 READ(u,rec=mid)x IF (x .EQ. target) THEN where=mid ELSE IF (x .GT. target) THEN hi=mid-1 ELSE lo=mid+1 END IF END DO ! END SUBROUTINE Bin_Search !------------------------------ END MODULE BinSearchSub !============================ PROGRAM SearchFile ! ! MS-DOS file name SRCHFILE.F90 ! Binary search on unformatted direct-access file. ! USE BinSearchSub IMPLICIT NONE INTEGER what,i,u,n_rec,n_bytes,where DATA u,n_rec/1,20/ ! INQUIRE(iolength=n_bytes)what OPEN(u,file='c:\ftn90.dir\srchfile.bin',access='direct', & form='unformatted',recl=n_bytes,action='readwrite') ! ! Fill file with even integers. DO i=1,n_rec WRITE(1,rec=i)2*i END DO ! 10 PRINT *,' Search for what, 999 to quit?' READ *,what CALL Bin_Search(1,n_rec,what,where) IF (what/=999) PRINT *,where IF (what/=999) GO TO 10 ! CLOSE(1) END ! Program 12.7 (COMPLEX.F90) PROGRAM Complex ! ! Demonstrate operations with COMPLEX data. ! IMPLICIT NONE REAL x,y,a,b COMPLEX c1,c2 ! x=1. y=2. a=0.5 b=3. c1=CMPLX(x,y) c2=CMPLX(a,b) PRINT *,c1,c2,c1+c2,c1*c2,c1/c2 PRINT *,REAL(c1),AIMAG(c1) END ! Program 12.8 (QUADRAT2.F90) PROGRAM Quadrat2 ! ! Calculates solutions to the quadratic equation. ! IMPLICIT NONE REAL a,b,c ! coefficients of ax^2+bx+c=0 COMPLEX root1,root2 ! two roots REAL discriminant ! PRINT*,' Give coefficients of ax^2+bx+c: ' READ*,a,b,c ! ! Test for existence of one or more roots. ! discriminant=b*b-4.0*a*c SELECT CASE (discriminant>0.) CASE (.TRUE.) root1=CMPLX((-b+SQRT(discriminant))/2./a,0.) root2=CMPLX((-b-SQRT(discriminant))/2./a,0.) PRINT*,' The two real roots are ',root1,root2 CASE (.FALSE.) IF (ABS(discriminant)<1e-7) THEN !assume discriminant=0 root1=-b/2./a PRINT*,' The single real root is ',root1 ELSE root1=CMPLX(-b/2./a,SQRT(-discriminant)/2/a) root2=CMPLX(-b/2./a,-SQRT(-discriminant)/2./a) PRINT*,' The two complex roots are ',root1,root2 ENDIF END SELECT END ! Program 12.9 (CIRCLCOM.F90) MODULE CircleFunctions ! CONTAINS ! REAL FUNCTION Area(radius) ! ! Do area calculation. ! IMPLICIT NONE REAL radius,pi INTENT(IN) radius COMMON pi ! Area=pi*radius*radius ! END FUNCTION Area ! REAL FUNCTION Circumference(radius) ! ! Do circumference calculation. ! IMPLICIT NONE REAL radius,pi INTENT(IN) :: radius COMMON pi ! Circumference=2.0*pi*radius ! END FUNCTION Circumference END MODULE ! PROGRAM CirclCom ! ! Calculate area and circumference of a circle, using ! two functions. ! USE CircleFunctions, ONLY : Area,Circumference IMPLICIT NONE REAL radius,pi COMMON pi ! pi=4.*ATAN(1.) PRINT *,' What is the radius of the circle?' READ *,radius ! PRINT 1000,Area(radius),Circumference(radius) ! 1000 FORMAT(1x,2f10.3) END ! Program 12.10 (COM_TEST.F90) MODULE com_blk_subs CONTAINS !------------------------ SUBROUTINE com_test IMPLICIT NONE REAL A(10),B(10) INTEGER i COMMON A,B ! DO i=1,10 PRINT 1000,A(i),B(i) END DO 1000 FORMAT(2f10.2) END SUBROUTINE com_test !-------------------------- SUBROUTINE com_test_2 IMPLICIT NONE INTEGER A(20) INTEGER i COMMON A ! DO i=1,20 PRINT 1000,A(i) END DO 1000 FORMAT(i20) END SUBROUTINE com_test_2 !---------------------------- END MODULE com_blk_subs !============================ PROGRAM com_blk ! ! Demonstrates use (and possible misuse) of COMMON blocks. ! USE com_blk_subs, ONLY : com_test,com_test_2 IMPLICIT NONE REAL A(20) INTEGER i COMMON A ! DO i=1,20 A(i)=i*10. END DO CALL com_test CALL com_test_2 END ! ! Appendix 4 (DATETIME.F90) PROGRAM DateTime ! ! Demonstrate use of intrinsic date and time procedures. ! IMPLICIT NONE CHARACTER date*8,time*10,zone*5 INTEGER values(8),count,count_rate,count_max ! CALL DATE_AND_TIME(date,time,zone,values) PRINT 1000,date,time,zone,values ! CALL SYSTEM_CLOCK(count,count_rate,count_max) PRINT 1010,count,count_rate,count_max 1000 FORMAT(a9,a11,a6,8i10) 1010 FORMAT (3i10) END