moduLe vars !setting up variables that wiLL be used for the lab, the givens and the counts !alpha,beta, and lamda are given from lab one, parameters have been set !maxit= max iterations, and n is the current iteration count integer,parameter::wp=selected_real_kind(15) real(wp),parameter:: alpha = 0.3204, beta= 2.75*.9, lamda=2.96 end module vars !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!! Program !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !Lab was created by Juan Sanchez for Engr 326 Program LabJuan use vars Implicit none !aLL variables must be declared before the interface integer::maxitP,nP,i,j Real(wp),dimension(8,8):: TP real(wp)::valueeP,ULP,LLP real(wp)::esp1P,percentile,q,percent 1 Format(F8.4,F8.4,F8.4,F8.4,F8.4,F8.4,F8.4,F8.4) 2 Format(F10.6) 3 Format(I4) ! ! ! !this is an interface for the function that will be used Interface function f(x) use vars implicit none real(wp), intent(in) :: x real(wp) :: f end function f end interface !This is the output file OPEN(Unit=1, File="Answers.txt",action="write", Status="replace") !initialize your variables, notice that I initialized mine to values but the user !will replace most of these with their preference. !I am opting to keep the lower limit as the exponent of lamda ULP=0 LLP=exp(lamda) maxitP=1000 esp1P=0.000001 Write(*,*)"Hello welcome to this program, this program is used for finding the design flow for & & a chlorination plant in order to treat a daily use of water to a specified percentile" Write(*,*)"To what percentile would you like this chlorination plant to treat?" Read(*,*) percentile Write(*,*)"After how many iterations would you like to stop this program?" Read(*,*) maxitP Write(*,*)"What do you want your accuracy tolerance to be" Read(*,*) esp1p Write(*,*)"What is your upper limit guess?" Read(*,*)ULP Write(1,*)"Your upper limit guess is" Write(1,2)ULP !So I realized that in order to get the correct value of the upper limit I was guessing and checking !I decided to let a do loop do the guess and check for me, as you can see in the do loop !depending on how far my percentile is I ask it to add or subtract values !As you can see the values I ask it to add are different than those that I ask it to !subtract this is so that I may converge on an answer rather than bounce back and forth do caLL romberg(LLP,ULP,ValueeP,Esp1P,MaxitP,nP,Tp) if (abs(valueeP-percentile) < esp1P) then percent=valueeP*100 write(1,*)"your upper limit is in the" , percent ,"perctile" Q=ulp write(*,*)"your upper limit is in the" , percent ,"perctile" Write(1,*)"Your ideal flow would be",Q,"m^3/day." Write(*,*)"Your ideal flow would be",Q,"m^3/day." exit else if (valueeP>percentile) then if (abs(valueeP-percentile)>(0.05)) then ulp = ulp - 5 else if ((abs(valueeP-percentile)<(0.05)) .or. (abs(valueeP-percentile)>(0.02))) then ulp = ulp - 0.05 else if (abs(valueeP-percentile)<(0.01)) then ulp = ulp - 0.001 end if else if (valueeP<percentile) then if (abs(valueeP-percentile)>(0.05)) then ulp = ulp + 3.5 else if ((abs(valueeP-percentile)<(0.05)) .or. (abs(valueeP-percentile)>(0.02))) then ulp = ulp + 0.03 else if (abs(valueeP-percentile)<(0.01)) then ulp = ulp + 0.001 end if end if end do Write(1,*) Write(1,*)"Your upper limit is" Write(1,2)ULP Write(1,*) Write(1,*)"Your lower limit is" Write(1,2)LLP Write(1,*) Write(1,*)"Your tolerance is" Write(1,2)ESP1P Write(1,*) Write(1,*)"Your max iterations is" Write(1,3)Maxitp Write(1,*) Do i=1,8 Write(1,1) Tp(i,:) Write(*,1) Tp(i,:) end do Write(*,*) Write(1,*)"The number of Iterations",Np Write(*,*)"The number of Iterations",Np ! I could in theory ask ask it to add less or more at a certain tolerance in order ! to get a more accurate value End Program LabJuan !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBROUTINE Romberg(LL,UL,valuee,esp1,maxit,n,T) use vars IMPLICIT NONE Real(wp),intent(in):: LL,esp1 Real(wp),intent(inout)::UL Real(wp),dimension(8,8),intent(out):: T Real(wp),intent(out)::valuee Integer,intent(in):: maxit Integer,intent(out)::n integer::nt Integer::r,i,j Real(wp ):: h,s1 ,s2 !Interface below is for the fucntion Interface Function f(x) IMPLICIT NONE Integer,parameter :: wp=selected_real_kind(15) Real(wp):: f Real(wp),intent(in):: x END FUNCTION END interface !Now inilize variables to zero, so the numbers don’t start out as junk nt=1 h=0 s1=0.0_wp s2=0.0_wp !Now to start the romberg subroutine r=2 !ratio of long to short stepsize h=(UL-LL)/(real((nt))) !stepsize s1=(f(LL)+f(UL))/2.0_wp DO i=1,nt-1 s2=s2+f(LL+real(i)*h) END DO T(1,1)=h*(s1+s2) n =1 !number of iterarions DO nt=2*nt !double the number of trapezoids h=(UL-LL)/nt !calculate the new stepsize DO i=1,nt/2 s2=s2+ f(LL+h*(2*i-1.0_wp )) !update s2 END DO T(n+1,1)=h*(s1+s2) DO j=2,n+1 !Do the extrapolations T(n+1,j)=(r**2*(j -1)* T(n+1,j-1)-T(n ,j-1))/(r**2*(j-1)-1) END DO IF ((abs(T(n+1,n)-T(n+1,n+1))) < esp1) THEN valuee=T(n+1,n+1) EXIT ELSE IF(n >maxit) THEN Write(*,*)"we dont messed up" END IF n=n +1 END DO END SUBROUTINE Romberg !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!! Funcions !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function f(x) !this (PDF) is given from the lab !since the whole equation is a mULtiplication equation it is easiest to break into 3 parts and ! mULtiply them together at the end function f is a composition of the mULtiplication of sections a,b,and c !section a is (1/(A*x*L(B) section b is (ln(x)*y)^B-1, and section c is e^((-ln(x)-y)/a) from the lab !where x> e^y use vars implicit none real(wp), intent(in) :: x real(wp) :: f, a, b, c a = 1/(alpha*x*Gamma(beta)) b = ((log(x) - lamda)**(beta-1))/(alpha**(beta-1)) c = exp(-1*((log(x)-lamda)/alpha)) f = a*b*c end function !function f(x) !implicit none !!this is a test function of homework 1 !integer,parameter::wp=selected_real_kind(15) ! !real(wp),intent(in)::x !real(wp)::f !f = (20-(5*exp(-0.1*x))-(10*exp(-0.2*(x-3)))) !end function