wiki:fortran

Version 2 (modified by Leon Kos, 4 years ago) (diff)

Initial set

FORTRAN 95 examples

Console1/Console1/Console1.f90

real x,dx,s
!open(1,file='y1.dat', status='unknown')
a=-0.5
b=0.5
n=100
dx=(b-a)/n 
S1=0
S2=0
S3=0
do i=0,n
       S1=S1+(f(i*dx+a)+f((i+1)*dx+a))/2*dx
       S2=S2+ f(i*dx+a)*dx 
if(mod(i,2)==0) then 
       S3=S3+dx/3*2*f(i*dx+a)
   else 
       S3=S3+dx/3*4*f(i*dx+a)
   end if       
             
end do 
print *, 'Trap = ', S1 , 'Rec=', S2 , 'Sim=', S3
pause

contains

function f(x) 
f=1/sqrt(1-x**2)
end function f
end

Console10/Console10/Console10.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

external fct, out
real aux(8,4)
common /a/  dt, k, n, vm, A, om, rl, x, y, r, u(1), u(2), u(3), u(4)

real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)


integer :: nd = 2
real :: dt = 2*3.14159265 
real :: t
real :: c = 3.0e10
!B0=me*c*om/ez 
B0 = 200.0
E0 = 5000
eq = 4.8e-10
me = 9.1e-28
om = eq*B0/me*c
A = 2.0
vm = A*om
n = 1 
k = 1
b = B0/B0
rl = c/om
u(1) = vm/vm
u(2) = vm/vm
u(3) = x/rl
u(4) = y/rl



 
open (unit=2, file = '1_rk.dat', status='unknown') 

write(*,*) 'x=', u(3), 'vx=', u(1), 'y=', u(4), 'vy=', u(2)
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih


end

subroutine fct (t, u, du)

real u(3), du(3), u(4), du(4)
common /a/ dt, k, n, vm, A, om, rl, x, y, r,  u(1), u(2), u(3), u(4)
  
du(1) = -u(3)
du(2) = u(4)
du(3) = u(1)
du(4) = u(2)              
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)

common /a/ dt, k, n, om, vm, A, rl, x, y, r,  u(1), u(2), u(3), u(4)

if (t.ge.k*(dt/20.)) then

 write (4,1) t, u(1), u(2), u(3), u(4)   
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               


Console11/Console11/Console11.f90

!  Console11.f90 
!
!  FUNCTIONS:
!  Console11 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console11
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console11

    implicit none

    ! Variables

    ! Body of Console11
    print *, 'Hello World'

    end program Console11

Console12/Console12/Console12.f90

!  Console12.f90 
!
!  FUNCTIONS:
!  Console12 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console12
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console12

    implicit none

    ! Variables

    ! Body of Console12
    print *, 'Hello World'

    end program Console12

Console13/Console13/Console13.f90

integer :: n=1000
real(8) :: dt, t, vm, b, rl, om, x, y, vx, vy, B0
real(8) :: eq=-4.8e-10, me=9.1e-28, c = 3.0e10
!B0=me*c*om/ez 
B0 = 200.0
!E0 = 5000
!A = 1e-30
x=1.0
y=1.0
vx=1000.0
vy=1000.0
a=100
!b = B0/B0
!om = eq*B0/me*c
!vm = A*om
!print *, vm
!rl = c/om
!vx = vm/vm
!vy = vm/vm
!x = x/rl
!y = y/rl
!print *, rl
t = 1.0e-8
dt = t/float(n)

open (1,file='1',status='unknown')

do i=1,n 
     vx = vx+(B0*eq)/(c*me)*vy*dt
     vy = vy-(B0*eq)/(c*me)*vx*dt
     x = x+vx*dt
     y = y+vy*dt
   write(1,*) x, vx
end do

pause
end program

Console15/Console15/Console15.f90

x = (/1,5/)
print *, x
end

Console16/Console16/Console16.f90

integer :: n=1000
real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 
x = 1.0
y = 1.0
a = 0.5

vx = 1.0e09
vy = 1.0e09
E0=0.0
B0 = 200.0

b=(2*3.14*6)/(1.0e09)

!B0 = 1e-25
!E0 = 5000
!A = 1e-30
!b = B0/B0
!om = eq*B0/me*c
!vm = A*om
!print *, vm
!rl = c/om
!vx = vm/vm
!vy = vm/vm
!x = x/rl
!y = y/rl
!print *, rl
t = 1.00e-08
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
   Bz = B0*(1+a*y) 
   !Bz = B0
   x = x + vx*dt
   y = y + vy*dt

   !vx = vx + (vy*eq*Bz)/(c*me)*dt
   !vy = vy - (vx*eq*Bz)/(c*me)*dt

   
   if(i*dt.gt.b) then 
E0 = 5.0
end if
   
  vx = vx +(eq*E0)/(me)*dt+(vy*eq*Bz)/(c*me)*dt
   vy = vy + (eq*E0)/(me)*dt-(vx*eq*Bz)/(c*me)*dt
  ! x = x + vx*dt
   !y = y + vy*dt

   !print *, vx, vy



  

          write (1,*) x, y 

    end do
 pause  
end program

Console17/Console17/Console17.f90

integer :: N = 1500
real(8) :: x = 2.5, y = 0.5, B0 = 500.0, We = 4.0054e-08, b = 2.0, pi = 3.14159, m = 9.1e-28, e = -4.8e-10, c = 2.99e10
real(8) vx, vy, v, dt, t

v = sqrt(2*We/m)
vx = sin(pi/4.0)*v
vy = cos(pi/4.0)*v
t = 5e-9
dt=t/float(N)
!Norma
!tau = t*10e8
!me = m*1e28
!qe = e*1e10
!vl = c*1e-10
!ve = v/c
!vx = sin(pi/4.0)*v
!vy = cos(pi/4.0)*v
!dt = tau/float(N)

open (1,file='mov',status='replace')
do i=1,N
   Bz = B0*(1-b**2*(x**2+y**2))
   x = x + vx*dt
   y = y + vy*dt
   vx = vx + vy/m*Bz*e/c*dt
   vy = vy - vx/m*Bz*e/c*dt
   print *, x, y
   pause
   write(1,*) y, x
end do
end program

Console18/Console18/Console18.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

external fct, out
real aux(8,2)
common /a/ b, lam, f0, dt, k, n

real :: pr(5)=(/0., 2*3.14159*20., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)
real :: b=0.003
integer :: nd = 2
          
real :: dt = 2*3.14159265 
   
! real :: lam = 0.1
real :: f0 = 0.00

n = 1 
k = 1
 
open (unit=2, file = 'osc_rk.dat', status='new') 

write(*,*) 'x=', u(1), 'v=', u(2)
pause
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih

stop
end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ b, lam, f0, dt, k, n

du(1) = -u(2)-2.*b*u(1)                  
du(2) = u(1)
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ b, lam, f0, dt, k, n

if (t.ge.k*(dt/20.)) then

 write (2,1) t,u(1),u(2)  ! 
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

Console19/Console19/Console19.f90

!  Osc_rk_0.f90 
!
!  Lab #3
!  metod Runge-Kuta. Oscillator.

external fct, out
real aux(8,2)
common /a/ b, lam, f0, dt, k, n, c 

real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)
real :: b=0.003
integer :: nd = 2
real::  c = 2.0
real :: dt = 2*3.14159265 

   
! real :: lam = 0.1
real :: f0 = 0.00

n = 1 
k = 1
 
open (unit=2, file = 'osc_rk.dat', status='new') 

write(*,*) 'x=', u(1), 'v=', u(2)
pause
call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=', ih

stop
end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ b, lam, f0, dt, k, n, c

du(1) = c*cos(0.5*t)-u(2)-2.*b*u(1)                  
du(2) = u(1)
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ b, lam, f0, dt, k, n, c

if (t.ge.k*(dt/20.)) then

 write (2,1) t,u(1),u(2)  ! 
  k = k + 1
   else
    end if

1 format (f8.3, f8.3, f8.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               

Console2/Console2/Console2.f90

real x,dx,s
!open(1,file='y1.dat', status='unknown')
a=0
b=1.57
n=100
dx=(b-a)/n
S=0
do i=0,n
    if(mod(i,2)==0) then 
       S=S+dx/3*2*f(i*dx+a)
   else 
       S=S+dx/3*4*f(i*dx+a)
   end if
end do
 
print *, 'Sim = ', S 
pause

contains

function f(x) 
f=sin(x)*sin(2*x)*sin(3*x)
end function f
end

Console20/Console20/Console20.f90

!  Console20.f90 
!
!  FUNCTIONS:
!  Console20 - Entry point of console application.
!

!****************************************************************************
!
!  PROGRAM: Console20
!
!  PURPOSE:  Entry point for the console application.
!
!****************************************************************************

    program Console20

    implicit none

    ! Variables

    ! Body of Console20
    print *, 'Hello World'

    end program Console20

Console21/Console21/Console21.f90

integer :: n=1000
real(8) :: dt, t, vm, om, x, y, vx, vy, B0,b
real(8) :: eq=4.8e-10, me=9.1e-28, c = 3.0e10, R=1.0 
x = 1.0
y = 1.0
a = 0.5

vx = 1.0e09
vy = 1.0e09
E0 = 0.0
B0 = 200.0

b = (2*3.14)/(1.0e09)

t = 1.00e-08
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
   Bz = B0*(1+a*y) 
   !Bz = B0
   x = x + vx*dt
   y = y + vy*dt

   vx = vx + (vy*eq*Bz)/(c*me)*dt
   vy = vy - (vx*eq*Bz)/(c*me)*dt

   
   if(i*dt.gt.b) then 
E0 =  3.0
   end if
   
  vx = vx + (vy*eq*Bz)/(c*me)*dt
  vy = vy + (eq*E0)/(me)*dt  - (vx*eq*Bz)/(c*me)*dt
  x = x - vx*dt
  y = y + vy*dt

   !print *, vx, vy



  

          write (1,*) x, y 

    end do
 pause  
end program

Console22/Console22/Console22.f90

integer :: n=1000
real(8) :: dt, t, x, y, vx, vy, v, E, W
real(8) ::  mp = 1.67e-24, q = 4.8e-10
x = 1.0
y = 1.0
d = 1.0
a = 30
W0 = 1000

v = sqrt(2*W0/mp)


t = 1.00e-05
dt = t/float(n)


open (1,file='1',status='unknown')

do i=1,n
   
  vx = v*sin(a)
  vy = v*cos(a)
  x = x + vx*dt
   y = y + vy*dt
  
  E = (4*3.14*mp*y)/d
  W = (E*q)/ x
   

   

   !print *, vx, vy



  
          write (1,*) W, t

    end do
 pause  
end program

Console24/Console24/Console24.f90

!  Lab #1
!  metod Runge-Kuta. Rezonans i avtorezonans.

external fct, out
real aux(8,2)
common /a/  dt, k, n, g0, om, B0, B, a 

real :: pr(5)=(/0., 2*3.14159*100., 2*3.14159/50, .0001, .0/)
real :: u(2) = (/0., 0.5/) 
real :: du(2) = (/0.5, 0.5/)

integer :: nd = 2
real :: dt = 2*3.14159265 

real :: t, m0, B, B0, a
real :: f0 = 0.00
f = 2.45e09
pi = 3.14159265
ez = 4.8e-10
m0 = 9.1e-28
c = 3.0e10
om = 2*pi*f
E = 2.
w = 1000.
a = 0.09
g0 = ez*E/(m0*c*om)
B0 = m0*c*om/ez
B = B0*(1+a*t)
b = B/B0
u(1) = pi
u(2) = w/511000.+1
write (*,*) g0
pause


 
n = 1 
k = 1 
 
open (unit=2, file = 'rez_avtorez.dat', status='unknown') 


call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=' ,ih 



end

subroutine fct (t, u, du)

real u(2), du(2), lam 
common /a/ dt, k, n, g0, om, B, B0, a

du(1) = (a*t+1-u(2))/u(2)+g0*sqrt(1./(u(2)**2-1))*sin(u(1))          
du(2) = -g0*sqrt(1.-1./u(2)**2)*cos(u(1))

!write (*,*)  du(1), du(2)
!pause
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(2), du(2), pr(5), lam
common /a/ dt, k, n, g0, om, B, B0, a

if (t.ge.k*(dt/2.)) then

 write (2,1) t,u(1),(u(2)-1)*511. 
  k = k + 1
   else
    end if

1 format (3e10.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(2),DERY(2),AUX(8,2),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39                                                           
   37 IHLF=12                                                           
      GOTO 39                                                           
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END                                                               


Console26/Console26/Console26.f90

!  Lab #1
!  metod Runge-Kuta. Rezonans i avtorezonans.

external fct, out
real aux(8,4)
common /a/  dt, k, n, g0, om, B0, B, a, y, rl 

real :: pr(5)=(/0., 2*3.14159*200., 2*3.14159/50, .0001, .0/)
real :: u(4) = (/0., 0., 0., 0./) 
real :: du(4) = (/0.25, 0.25, 0.25, 0.25/)

integer :: nd = 4
real :: dt = 2*3.14159265 

real :: t, m0, B, B0, a
real :: f0 = 0.00
f = 2.45e09
pi = 3.14159265
ez = 4.8e-10
m0 = 9.1e-28
c = 3.0e10
om = 2*pi*f
E0 = 5.
E = E0
a = 0.00034
B0 = m0*c*om/ez
g0 = E0/B0
B = B0
b = B/B0
rl = c/om
u(1) = 0.
u(2) = 0.
u(3) = 0.
u(4) = 0.
write (*,*) g0, B0, b, E, rl
pause


 
n = 1 
k = 1 
 
open (unit=2, file = 'rez_avtorez.dat', status='unknown') 


call rkgs(pr, u, du, nd, ih, fct, out, aux)

write(2,*) 'ih=' ,ih 



end

subroutine fct (t, u, du)

real u(4), du(4), lam 
common /a/ dt, k, n, g0, om, B, B0, a, y, rl
y=sqrt(1+u(1)**2+u(2)**2)
du(1) = -g0*cos(t)- u(2)*(1+a*t)/y       
du(2) = -g0*sin(t)+u(1)*(1+a*t)/y
du(3) = u(1)/y
du(4) = u(2)/y

!write (*,*)  du(1), du(2), du(3), du(4)
!pause
                                
return 
end  

subroutine  out (t, u, du, ih, nd, pr)
real u(4), du(4), pr(5), lam
common /a/ dt, k, n, g0, om, B, B0, a, y, rl

if (t.ge.k*(dt)) then

 write (2,1) t,(y-1)*511. ,u(3)*rl,u(4)*rl 
  k = k + 1
   else
    end if

1 format (5e12.3)

return
end 


      SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)               
!                                                                       
!                                                                       
      DIMENSION Y(4),DERY(4),AUX(8,4),A(4),B(4),C(4),PRMT(5)            
      DO 1 I=1,NDIM                                                     
    1 AUX(8,I)=.06666667*DERY(I)                                        
      X=PRMT(1)                                                         
      XEND=PRMT(2)                                                      
      H=PRMT(3)                                                         
      PRMT(5)=0.                                                        
      CALL FCT(X,Y,DERY)                                                
!                                                                       
!     ERROR TEST                                                        
      IF(H*(XEND-X))38,37,2                                             
!                                                                       
!     PREPARATIONS FOR RUNGE-KUTTA METHOD                               
    2 A(1)=.5                                                           
      A(2)=.2928932                                                     
      A(3)=1.707107                                                     
      A(4)=.1666667                                                     
      B(1)=2.                                                           
      B(2)=1.                                                           
      B(3)=1.                                                           
      B(4)=2.                                                           
      C(1)=.5                                                           
      C(2)=.2928932                                                     
      C(3)=1.707107                                                     
      C(4)=.5                                                           
!                                                                       
!     PREPARATIONS OF FIRST RUNGE-KUTTA STEP                            
      DO 3 I=1,NDIM                                                     
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=0.                                                       
    3 AUX(6,I)=0.                                                       
      IREC=0                                                            
      H=H+H                                                             
      IHLF=-1                                                           
      ISTEP=0                                                           
      IEND=0                                                            
!                                                                      
!                                                                   
!     START OF A RUNGE-KUTTA STEP                                       
    4 IF((X+H-XEND)*H)7,6,5                                             
    5 H=XEND-X                                                          
    6 IEND=1                                                            
!                                                                      
!     RECORDING OF INITIAL VALUES OF THIS STEP                          
    7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)                                
      IF(PRMT(5))40,8,40                                                
    8 ITEST=0                                                           
    9 ISTEP=ISTEP+1                                                     
!     START OF INNERMOST RUNGE-KUTTA LOOP                               
      J=1                                                               
   10 AJ=A(J)                                                           
      BJ=B(J)                                                           
      CJ=C(J)                                                           
      DO 11 I=1,NDIM                                                    
      R1=H*DERY(I)                                                      
      R2=AJ*(R1-BJ*AUX(6,I))                                            
      Y(I)=Y(I)+R2                                                      
      R2=R2+R2+R2                                                       
   11 AUX(6,I)=AUX(6,I)+R2-CJ*R1                                        
      IF(J-4)12,15,15                                                   
   12 J=J+1                                                             
      IF(J-3)13,14,13                                                   
   13 X=X+.5*H                                                          
   14 CALL FCT(X,Y,DERY)                                                
      GOTO 10                                                           
!     END OF INNERMOST RUNGE-KUTTA LOOP                                 
!     TEST OF ACCURACY                                                  
   15 IF(ITEST)16,16,20                                                 
!                                                                       
!     IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY   
   16 DO 17 I=1,NDIM                                                    
   17 AUX(4,I)=Y(I)                                                     
      ITEST=1                                                           
      ISTEP=ISTEP+ISTEP-2                                               
   18 IHLF=IHLF+1
      X=X-H                                                             
      H=.5*H                                                            
      DO 19 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
      DERY(I)=AUX(2,I)                                                  
   19 AUX(6,I)=AUX(3,I)                                                 
      GOTO 9                                                            
!                                                                       
!     IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE                   
   20 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)21,23,21                                       
   21 CALL FCT(X,Y,DERY)                                                
      DO 22 I=1,NDIM                                                    
      AUX(5,I)=Y(I)                                                     
   22 AUX(7,I)=DERY(I)                                                  
      GOTO 9                                                            
!                                                                       
!     COMPUTATION OF TEST VALUE DELT                                    
   23 DELT=0.                                                           
      DO 24 I=1,NDIM                                                    
   24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))                             
      IF(DELT-PRMT(4))28,28,25                                          
!                                                                       
!     ERROR IS TOO GREAT                                                
   25 IF(IHLF-10)26,36,36                                               
   26 DO 27 I=1,NDIM                                                    
   27 AUX(4,I)=AUX(5,I)                                                 
      ISTEP=ISTEP+ISTEP-4                                               
      X=X-H                                                             
      IEND=0                                                            
      GOTO 18                                                           
!                                                                       
!     RESULT VALUES ARE GOOD                                            
   28 CALL FCT(X,Y,DERY)                                                
      DO 29 I=1,NDIM                                                    
      AUX(1,I)=Y(I)                                                     
      AUX(2,I)=DERY(I)                                                  
      AUX(3,I)=AUX(6,I)                                                 
      Y(I)=AUX(5,I)                                                     
   29 DERY(I)=AUX(7,I)                                                  
      CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)                              
      IF(PRMT(5))40,30,40                                               
   30 DO 31 I=1,NDIM                                                    
      Y(I)=AUX(1,I)                                                     
   31 DERY(I)=AUX(2,I)                                                  
      IREC=IHLF                                                         
      IF(IEND)32,32,39                                                  
!                                                                       
!     INCREMENT GETS DOUBLED                                            
   32 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      IF(IHLF)4,33,33                                                   
   33 IMOD=ISTEP/2                                                      
      IF(ISTEP-IMOD-IMOD)4,34,4                                         
   34 IF(DELT-.02*PRMT(4))35,35,4                                       
   35 IHLF=IHLF-1                                                       
      ISTEP=ISTEP/2                                                     
      H=H+H                                                             
      GOTO 4                                                            
!                                                                       
!     RETURNS TO CALLING PROGRAM                                        
   36 IHLF=11                                                           
      CALL FCT(X,Y,DERY)                                                
      GOTO 39      
   38 IHLF=13                                                           
   39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)                                
   40 RETURN                                                            
      END