Program d02qgfe ! D02QGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: d02qgf, d02qwf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Real (Kind=nag_wp), Parameter :: alpha = -0.032_nag_wp Real (Kind=nag_wp), Parameter :: beta = -0.02_nag_wp Integer, Parameter :: neqf = 3, neqg = 1, nin = 5, nout = 6 Integer, Parameter :: liwork = 21 + 4*neqg Integer, Parameter :: lrwork = 23 + 23*neqf + 14*neqg ! .. Local Scalars .. Real (Kind=nag_wp) :: grvcm, hmax, t, tcrit, tinc, tout, & trvcm, tstart Integer :: i, ifail, irevcm, j, kgrvcm, latol, & lrtol, maxstp, yprvcm, yrvcm Logical :: alterg, crit, onestp, root, sophst, & vectol Character (1) :: statef ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: atol(:), rtol(:), rwork(:), y(:) Integer, Allocatable :: iwork(:) ! .. Intrinsic Procedures .. Intrinsic :: cos, tan ! .. Executable Statements .. Write (nout,*) 'D02QGF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) latol, lrtol Allocate (atol(latol),rtol(lrtol),rwork(lrwork),y(neqf),iwork(liwork)) Read (nin,*) hmax, tstart, tcrit, tinc Read (nin,*) statef Read (nin,*) vectol, onestp, crit, sophst Read (nin,*) maxstp Read (nin,*) rtol(1:lrtol), atol(1:latol) Read (nin,*) y(1:neqf) t = tstart ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call d02qwf(statef,neqf,vectol,atol,latol,rtol,lrtol,onestp,crit,tcrit, & hmax,maxstp,neqg,alterg,sophst,rwork,lrwork,iwork,liwork,ifail) Write (nout,*) Write (nout,*) ' T Y(1) Y(2) Y(3)' Write (nout,99999) t, (y(i),i=1,neqf) j = 1 tout = tinc irevcm = 0 revcm: Do ifail = -1 Call d02qgf(neqf,t,y,tout,neqg,root,irevcm,trvcm,yrvcm,yprvcm,grvcm, & kgrvcm,rwork,lrwork,iwork,liwork,ifail) Select Case (irevcm) Case (0) If (ifail==0) Then ! Print solution at current t Write (nout,99999) t, (y(i),i=1,neqf) If (t==tout .And. j<5) Then ! Increment tout and cycle to find solution at this new time. j = j + 1 tout = tout + tinc Cycle revcm End If End If Exit revcm Case (1:7) If (yrvcm==0) Then rwork(yprvcm) = tan(y(3)) rwork(yprvcm+1) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3)) rwork(yprvcm+2) = alpha/y(2)**2 Else rwork(yprvcm) = tan(rwork(yrvcm+2)) rwork(yprvcm+1) = alpha*tan(rwork(yrvcm+2))/rwork(yrvcm+1) + & beta*rwork(yrvcm+1)/cos(rwork(yrvcm+2)) rwork(yprvcm+2) = alpha/rwork(yrvcm+1)**2 End If Case (9:) grvcm = y(1) End Select End Do revcm 99999 Format (1X,F7.4,2X,3(F7.4,2X)) End Program d02qgfe