Run time error after any correction

Hi everyone. I have got strange problem with program showed below. When I start original version of program it works fine- I mean it shows something at least. After puting any correction it shows RUN-TIME ERROR 11- also showed below. Here is the program:

 PROGRAM MODEL_KONWEKCJI_SEANS
  INTEGER i,N
  REAL ROH,MK,PI,MI,NI,RE,PR
  REAL OS,OG,OO,Z,RO,AKP,W,DH
  REAL KD,BF1,BF2,S
  REAL NU0,NUS,NUGW,LAMBDAP,ALKS,ALKGW
  REAL X,X0,TP,TS,TG,TH
  DIMENSION DELTAX(16)
  DIMENSION DG(16)
  DIMENSION DS(16)
  DIMENSION TP(16)
  DIMENSION X(16)
  WRITE (*,*)' PROGRAM DO SYMULACJI WYMIANY CIEPLA NA DRODZE KONWEKCJI'
  PAUSE
  WRITE (*,*) 'Podaj wartosc N,z zakresu 1-16'
  READ(*,*)N
  WRITE (*,*) ' '
    
  !WYRAZENIA GEOMETRYCZNE!
  X0=0
  DELTAX(1)=0.106
  DELTAX(2)=0.108
  DELTAX(3)=0.054
  DELTAX(4)=0.054
  DELTAX(5)=0.049
  DELTAX(6)=0.072
  DELTAX(7)=0.048
  DELTAX(8)=0.031
  DELTAX(9)=0.091
  DELTAX(10)=0.061
  DELTAX(11)=0.121
  DELTAX(12)=0.070
  DELTAX(13)=0.073
  DELTAX(14)=0.104
  DELTAX(15)=0.036
  DELTAX(16)=0.096
  


  DG(1)=0.658
  DG(2)=0.666
  DG(3)=0.668
  DG(4)=0.670
  DG(5)=0.686
  DG(6)=0.687
  DG(7)=0.686
  DG(8)=0.684
  DG(9)=0.679
  DG(10)=0.676
  DG(11)=0.664
  DG(12)=0.652
  DG(13)=0.638
  DG(14)=0.614
  DG(15)=0.603
  DG(16)=0.569
  

  DS(1)=0.276
  DS(2)=0.232
  DS(3)=0.226
  DS(4)=0.222
  DS(5)=0.266
  DS(6)=0.357
  DS(7)=0.357
  DS(8)=0.357
  DS(9)=0.357
  DS(10)=0.357
  DS(11)=0.357
  DS(12)=0.364
  DS(13)=0.376
  DS(14)=0.395
  DS(15)=0.402
  DS(16)=0.420

  !!!!! DANE NIEZALENE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ROH=1.18
  TH=293
  TP(1)=293
  PI=3.141592654
  PR=0.7
  
  DO i=1,N

** X(i)=X(i-1)+DELTAX(i)**

  MK=0.2
  OS=PI*(DS(i))
  OG=PI*(DG(i))
  OO=OS+OG
  Z=TH*ROH																					
  RO=Z/(TP(i))
  AKP=(PI/4)*((DG(i)**2)-(DS(i)**2))
  W=MK/(RO*AKP)
  DH=4*AKP/OO
  MI=6.0634+0.0417*(TP(i))
  NI=(MI/RO)
  RE=(10**6)*(DH*W)/NI
  KD=DS(i)/DG(i)
  BF1=1.6075*(KD**3)-3.273*(KD**2)+2.4078*KD+0.4388
  BF2=-0.8089*(KD**2)+2.2256*KD+0.7013
  S=0.0757*(KD**3)-0.1685*(KD**2)+0.1556*KD+0.0374
  NU0=(0.023*(RE**0.8)*PR)/(1+(2.09/(RE**0.1))*((PR**0.667)-1))
  NUS= NU0*(1+(2.09/(RE**0.1))*(((PR**0.667)-1)))/(1+0.1741*KD+((2.09+0.05*KD)/(RE**0.1))*((PR**0.667)-1))
  NUGW=NU0*(BF1+(BF2/(RE**S))*((PR**0.667)-1)/(1+2.09/(RE**0.1))*((PR**0.667)-1))
  LAMBDAP=0.001*(5.73+0.06638*TP(i))
  ALKS=NUS*LAMBDAP/DH
  ALKGW=NUGW*LAMBDAP/DH
  CP=978.2+0.1178*TP(i)

  END DO

  

!!!!!!!!!!!!!!!!!!!!!!!OBLICZENIA TEMPERATUR POW. SILNIKA I POW. GONDOLI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  WRITE (*,*) 'Podaj nr seansu badawczego, S z zakresu 1-4:'
  READ (*,*) S
  
  !SEANS 1!
  IF (N.LT.2.and.S.EQ.1) THEN
TS=475.13*(X(i)**2)-40.929*X(i)+320
TG=(-526.02)*(X(i)**6)+2020.3*(X(i)**5)-2690*(X(i)**4)+1349.7*(X(i)**3)-138.25*(x(i)**2)+9.5664*X(i)+294.01

  end if
  IF (N.GE.2.and.S.EQ.1) THEN
TS=(-45547)*(X(i)**6)+173798*(X(i)**5)-251552*(X(i)**4)+171250*(X(i)**3)-55740*(X(i)**2)+8476.7*X(i)-150.38
TG=(-526.02)*(X(i)**6)+2020.3*(X(i)**5)-2690*(X(i)**4)+1349.7*(X(i)**3)-138.25*(x(i)**2)+9.5664*X(i)+294.01

  END IF

  IF (i.EQ.1) THEN 
  TP=TP(i)
  END IF 
  IF (i.GE.1) THEN 
  TP(i)=TP(i-1)+(TS-TP(i-1))*DELTAX(i)*0.001*(NUS/DH)*(5.73+0.06638*TP(i-1))*OS/(MK*(978.2+0.1178*TP(i-1)))+(TG-TP(i))*DELTAX(i)*0.001*(NUGW/DH)*(5.73+0.06638*TP(i-1))*OG/(MK*(978.2+0.1178*TP(i-1)))

  END IF

!!!!!!!!!!!!WYNIKI!!!!!!!!!!!!!! WRITE (,)'Wyniki obliczen' WRITE (,)'WYNIKI OBLICZEN DLA PRZEKROJU i=',i-1 WRITE (,)'***********************************************' WRITE (,)'X',X(i-1) WRITE (,)'RE',RE WRITE (,)'ALKS,ALKGW',ALKS,ALKGW WRITE (,)'TS',TS WRITE (,)'TG',TG WRITE (,)'TP',TP

  END PROGRAM

and here is the error message:

Run-time Error *** Error 11, Array subscript(s) out-of-bounds

MODEL_KONWEKCJI_SEANS - in file infred_strs-g_1.0.f95 at line 83 [+0499

Error line is bolded. I am new in FORTRAN. Pleas help me if you know the problem. Probably it is obvious but not for me.

Your immediate error in the change is that your variables are not initialised. DO i=1,N X(i)=X(i-1)+DELTAX(i)

This fails for i=1, as x(0) is not defined. To fix this, I changed the definition of X as 'DIMENSION X(0:16)' and initialised X = 0

The complexity of the equations you have, made me want to tidy them up. If you ever need to come back to this code to change it, the following block of code would have you studying for hours, as to what it means. My advice is to simplify it: NU0=(0.023*(RE0.8)*PR)/(1+(2.09/(RE0.1))((PR**0.667)-1)) NUS= NU0(1+(2.09/(RE0.1))*(((PR0.667)-1)))/(1+0.1741KD+((2.09+0.05KD)/(RE0.1))*((PR0.667)-1)) NUGW=NU0*(BF1+(BF2/(RES))*((PR0.667)-1)/(1+2.09/(RE0.1))*((PR0.667)-1)) LAMBDAP=0.001*(5.73+0.06638*TP(i))

I also tried to simplify your code by providing a generic polynomial function. Hopefully it works ! It does at least make the polynomial definitions easier to check.

There is another major problem as the code to reset TP is outside the DO i loop. This needs to be corrected. IF (i.EQ.1) THEN TP=TP(i) END IF IF (i.GE.1) THEN ( should be i.gt. or at least changed for TP(0) )

Anyway, I have made some changes, but there is more to be done. Most of the changes I have suggested are layout suggestions and not essential. My advice is that you may need to come back and modify the code in 6 months time, so provide a layout that helps with understanding the basis of the equations. Layout is a personal preference and you choose a layout that is easier for you to understand. I hope the changes I have suggested might give you some ideas of a different layout.

John

module polynomial_kinds
   implicit none
   interface polynomial
      module procedure poly4, poly8
   end interface polynomial
   contains
      real*4 function poly4 (coeff, x)
         real*4, dimension(:) :: coeff
         real*4               :: x
!         
         real*8    s
         integer*4 i, n
!
         write (*,*) 'Poly_4',size(coeff)-1
         n = size (coeff)
         s = coeff(n)
         do i = n-1,1,-1
           s = s * x + coeff(i)
         end do
         poly4 = s
      end function poly4
!
      real*8 function poly8 (coeff, x)
         real*8, dimension(:) :: coeff
         real*8               :: x
!         
         real*10   s
         integer*4 i, n
!
         write (*,*) 'Poly_8',size(coeff)-1
         n = size (coeff)
         s = coeff(n)
         do i = n-1,1,-1
           s = s * x + coeff(i)
         end do
         poly8 = s
      end function poly8
!      
end module polynomial_kinds

PROGRAM MODEL_KONWEKCJI_SEANS 
!
use  polynomial_kinds
!
INTEGER i,N 
REAL ROH,MK,PI,MI,NI,RE,PR 
REAL OS,OG,OO,Z,RO,AKP,W,DH 
REAL KD,BF1,BF2,S 
REAL NU0,NUS,NUGW,LAMBDAP,ALKS,ALKGW 
REAL X,X0,TP,TS,TG,TH 
DIMENSION DELTAX(16) 
DIMENSION DG(16) 
DIMENSION DS(16) 
DIMENSION TP(16) 
DIMENSION X(0:16) 
REAL DELTAX,DG,DS,CP
!
WRITE (*,*)' PROGRAM DO SYMULACJI WYMIANY CIEPLA NA DRODZE KONWEKCJI' 
PAUSE 
WRITE (*,*) 'Podaj wartosc N,z zakresu 1-16' 
READ(*,*)N 
WRITE (*,*) ' ' 

!WYRAZENIA GEOMETRYCZNE! 
X0=0 
DELTAX(1)=0.106 
DELTAX(2)=0.108 
DELTAX(3)=0.054 
DELTAX(4)=0.054 
DELTAX(5)=0.049 
DELTAX(6)=0.072 
DELTAX(7)=0.048 
DELTAX(8)=0.031 
DELTAX(9)=0.091 
DELTAX(10)=0.061 
DELTAX(11)=0.121 
DELTAX(12)=0.070 
DELTAX(13)=0.073 
DELTAX(14)=0.104 
DELTAX(15)=0.036 
DELTAX(16)=0.096 
!

rest of changes:

!
DG(1)=0.658 
DG(2)=0.666 
DG(3)=0.668 
DG(4)=0.670 
DG(5)=0.686 
DG(6)=0.687 
DG(7)=0.686 
DG(8)=0.684 
DG(9)=0.679 
DG(10)=0.676 
DG(11)=0.664 
DG(12)=0.652 
DG(13)=0.638 
DG(14)=0.614 
DG(15)=0.603 
DG(16)=0.569 
!
DS(1)=0.276 
DS(2)=0.232 
DS(3)=0.226 
DS(4)=0.222 
DS(5)=0.266 
DS(6)=0.357 
DS(7)=0.357 
DS(8)=0.357 
DS(9)=0.357 
DS(10)=0.357 
DS(11)=0.357 
DS(12)=0.364 
DS(13)=0.376 
DS(14)=0.395 
DS(15)=0.402 
DS(16)=0.420 

!!!!! DANE NIEZALENE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
ROH   = 1.18 
TH    = 293 
TP(1) = 293 
PI    = 4.0 * atan (1.0) ! 3.141592654 
PR    = 0.7 
x     = 0
!
DO i=1,N 
  X(i)  = X(i-1)+DELTAX(i) 
  MK    = 0.2 
  OS    = PI*(DS(i)) 
  OG    = PI*(DG(i)) 
  OO    = OS+OG 
  Z     = TH*ROH 
  RO    = Z/(TP(i)) 
  AKP   = (PI/4)*((DG(i)**2)-(DS(i)**2)) 
  W     = MK/(RO*AKP) 
  DH    = 4*AKP/OO 
  MI    = 6.0634+0.0417*(TP(i)) 
  NI    = (MI/RO) 
  RE    = (10**6)*(DH*W)/NI 
  KD    = DS(i)/DG(i) 
!
!  BF1=1.6075*(KD**3)-3.273*(KD**2)+2.4078*KD+0.4388 
!  BF2=-0.8089*(KD**2)+2.2256*KD+0.7013 
!  S=0.0757*(KD**3)-0.1685*(KD**2)+0.1556*KD+0.0374 
  BF1   = polynomial ( (/ 0.4388, 2.4078, -3.273, 1.6075 /), KD)
  BF2   = polynomial ( (/ 0.7013, 2.2256, -0.8089 /), KD)
  BF2   = polynomial ( (/ 0.7013d0, 2.2256d0, -0.8089d0 /), dble(KD) )
  S     = polynomial ( (/ 0.0374, 0.1556, -0.1685, 0.0757 /), KD)
!
! I would rewrite the following so that the equation layout looks more like the written equation
! this could be done by using spaces or temproary variables
! consider replacing integers with reals, by inserting a '.'
! eg 173798 as 173798  or 10**6 as 10.e6
!
  NU0   = (0.023*(RE**0.8)*PR)/(1+(2.09/(RE**0.1))*((PR**0.667)-1)) 
  NUS   = NU0*(1+(2.09/(RE**0.1))*(((PR**0.667)-1)))/(1+0.1741*KD+((2.09+0.05*KD)/(RE**0.1))*((PR**0.667)-1)) 
  NUGW  = NU0*(BF1+(BF2/(RE**S))*((PR**0.667)-1)/(1+2.09/(RE**0.1))*((PR**0.667)-1)) 
  LAMBDAP = 0.001*(5.73+0.06638*TP(i)) 
  ALKS  = NUS*LAMBDAP/DH 
  ALKGW = NUGW*LAMBDAP/DH 
  CP    = 978.2+0.1178*TP(i) 
!
END DO   !  why does teh do loop end here, as TP is defined below


!!!!!!!!!!!!!!!!!!!!!!!OBLICZENIA TEMPERATUR POW. SILNIKA I POW. GONDOLI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 

WRITE (*,*) 'Podaj nr seansu badawczego, S z zakresu 1-4:' 
READ (*,*) S 

!SEANS 1! 
IF (N.LT.2.and.S.EQ.1) THEN 
!  TS=475.13*(X(i)**2)-40.929*X(i)+320 
!  TG=(-526.02)*(X(i)**6)+2020.3*(X(i)**5)-2690*(X(i)**4)+1349.7*(X(i)**3)-138.25*(x(i)**2)+9.5664*X(i)+294.01 
  TS = polynomial ( (/ 320.0, -40.929, 475.13 /), X(i) )
  TG = polynomial ( (/ 294.01, 9.5664, -138.25, 1349.7, -2690.0, 2020.3, -526.02 /), X(i) )

end if 
IF (N.GE.2.and.S.EQ.1) THEN 
!  TS=(-45547)*(X(i)**6)+173798*(X(i)**5)-251552*(X(i)**4)+171250*(X(i)**3)-55740*(X(i)**2)+8476.7*X(i)-150.38 
!  TG=(-526.02)*(X(i)**6)+2020.3*(X(i)**5)-2690*(X(i)**4)+1349.7*(X(i)**3)-138.25*(x(i)**2)+9.5664*X(i)+294.01 
  TS = polynomial ( (/ -150.38, 8476.7, -55740.0, 171250.0, -251552.0, 173798.0, -45547.0 /), X(i) )
  TG = polynomial ( (/  294.01, 9.5664,  -138.25,   1349.7,   -2690.0,   2020.3,  -526.02 /), X(i) )
END IF 

IF (i.EQ.1) THEN 
  TP = TP(i) 
ELSE
!IF (i.GE.1) THEN 
  TP(i) = TP(i-1)                                                                                       &
        + (TS-TP(i-1))*DELTAX(i)*0.001*(NUS/DH) *(5.73+0.06638*TP(i-1))*OS/(MK*(978.2+0.1178*TP(i-1)))  &
        + (TG-TP(i))  *DELTAX(i)*0.001*(NUGW/DH)*(5.73+0.06638*TP(i-1))*OG/(MK*(978.2+0.1178*TP(i-1))) 
END IF 

!!!!!!!!!!!!WYNIKI!!!!!!!!!!!!!! 
WRITE (*,*)'Wyniki obliczen' 
WRITE (*,*)'WY'
end

When I display these posts, the layout is squeezed to the right. Any ideas why this is the case? It's as if the first post did not finish correctly. The post limit is again looking very unreasonable !

John

Here is a stand alone program that demonstrates POLYNOMIAL working for real4 and real8. Not sure how to provide a real*10 constant argument.

module polynomial_kinds
   implicit none
   interface polynomial
      module procedure poly4, poly8, poly10
   end interface polynomial
   contains
      real*4 function poly4 (coeff, x)
         real*4, dimension(:) :: coeff
         real*4               :: x
!         
         real*10   s
         integer*4 i, n
!
         write (*,*) 'Poly_4',size(coeff)-1
         n = size (coeff)
         s = coeff(n)
         do i = n-1,1,-1
           s = s * x + coeff(i)
         end do
         poly4 = s
      end function poly4
!
      real*8 function poly8 (coeff, x)
         real*8, dimension(:) :: coeff
         real*8               :: x
!         
         real*10   s
         integer*4 i, n
!
         write (*,*) 'Poly_8',size(coeff)-1
         n = size (coeff)
         s = coeff(n)
         do i = n-1,1,-1
           s = s * x + coeff(i)
         end do
         poly8 = s
      end function poly8
!
      real*10 function poly10 (coeff, x)
         real*10, dimension(:) :: coeff
         real*10               :: x
!         
         real*10   s
         integer*4 i, n
!
         write (*,*) 'Poly_8',size(coeff)-1
         n = size (coeff)
         s = coeff(n)
         do i = n-1,1,-1
           s = s * x + coeff(i)
         end do
         poly10 = s
      end function poly10
!      
end module polynomial_kinds

PROGRAM MODEL_KONWEKCJI_SEANS 
!
use  polynomial_kinds
!
INTEGER i,N 
real*8 pi
real*8 bf4
REAL KD,BF1,BF2,S, bf3
DIMENSION DG(6) 
DIMENSION DS(6) 
REAL DG,DS
!
n = 5
!
DG(1)=0.658 
DG(2)=0.666 
DG(3)=0.668 
DG(4)=0.670 
DG(5)=0.686 
DG(6)=0.687 
!
DS(1)=0.276 
DS(2)=0.232 
DS(3)=0.226 
DS(4)=0.222 
DS(5)=0.266 
DS(6)=0.357 
!
PI    = 4.0d0 * atan (1.0d0) ! 3.141592654 
write (*,*) 'pi=',pi
!
DO i=1,N 
  write (*,*) 'i=',i
  KD    = DS(i)/DG(i) 
  write (*,*) 'kd=',kd
!
  BF1   = polynomial ( (/ 0.4388, 2.4078, -3.273, 1.6075 /), KD)
  BF2   = polynomial ( (/ 0.7013, 2.2256, -0.8089 /), KD)
  BF3   = polynomial ( (/ 0.7013d0, 2.2256d0, -0.8089d0 /), dble(KD) )  ! r*4 = r*8
  BF4   = polynomial ( (/ 0.7013d0, 2.2256d0, -0.8089d0 /), dble(KD) )  ! r*8 = r*8
  S     = polynomial ( (/ 0.0374, 0.1556, -0.1685, 0.0757 /), KD)
!
   write (*,*) 'bf1=', bf1
   write (*,*) 'bf2=', bf2
   write (*,*) 'bf3=', bf3
   write (*,*) 'bf4=', bf4
   write (*,*) 's=', s
!
END DO

end

Thank you John for your answer. I appreciate your help. As you pointed there is another big problem in this program with TP(i) Now I wonder how can I deal with that. I can not put it in DO i loop because TP(i) is inside IF THEN- or maybe i am wrong and can do that? This time program shows error:

Run-time Error *** Error 112, Reference to undefined variable, array element or function result (/UNDEF)

MODEL_KONWEKCJI_SEANS - in file infred_strs-g_1.0.f95 at line 132 [+18f0]

I understand why does it appear but how can deal with this if TP(i) depends on TS and TG;/ Any ideas? Maybe the whole idea of this program is wrong? I hope not.

Her is the last few lines of the program which first post didnt display corectly.

I understand why does it appear but how can deal with this if TP(i) depends on TS and TG;/ Any ideas? Maybe the whole idea of this program is wrong? I hope not.

Her is the last few lines of the program which first post didnt display corectly. [quote:17db251ca0]!!!!! DANE NIEZALENE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ROH=1.18 TH=293 TP(1)=293 PI=3.141592654 PR=0.7

DO i=1,N X(i)=X(i-1)+DELTAX(i) MK=0.2 OS=PI*(DS(i)) OG=PI*(DG(i)) OO=OS+OG Z=THROH RO=Z/(TP(i)) AKP=(PI/4)((DG(i)2)-(DS(i)*2)) W=MK/(ROAKP) DH=4AKP/OO MI=6.0634+0.0417(TP(i)) NI=(MI/RO) RE=(106)(DHW)/NI KD=DS(i)/DG(i) BF1=1.6075*(KD3)-3.273*(KD2)+2.4078KD+0.4388 BF2=-0.8089(KD2)+2.2256KD+0.7013 S=0.0757(KD3)-0.1685*(KD2)+0.1556KD+0.0374 NU0=(0.023(RE0.8)PR)/(1+(2.09/(RE**0.1))((PR0.667)-1)) NUS= NU0*(1+(2.09/(RE0.1))(((PR**0.667)-1)))/(1+0.1741KD+((2.09+0.05KD)/(RE**0.1))((PR0.667)-1)) NUGW=NU0*(BF1+(BF2/(RES))((PR0.667)-1)/(1+2.09/(RE0.1))((PR**0.667)-1)) LAMBDAP=0.001*(5.73+0.06638TP(i)) ALKS=NUSLAMBDAP/DH ALKGW=NUGWLAMBDAP/DH CP=978.2+0.1178TP(i)

END DO

!!!!!!!!!!!!!!!!!!!!!!!OBLICZENIA TEMPERATUR POW. SILNIKA I POW. GONDOLI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

WRITE (,) 'Podaj nr seansu badawczego, S z zakresu 1-4:' READ (,) S

!SEANS 1! IF (N.LT.2.and.S.EQ.1) THEN TS=475.13*(X(i)**2)-40.929X(i)+320 TG=(-526.02)(X(i)*6)+2020.3(X(i)*5)-2690(X(i)*4)+1349.7(X(i)3)-138.25(x(i)**2)+9.5664X(i)+294.01

end if IF (N.GE.2.and.S.EQ.1) THEN TS=(-45547)*(X(i)*6)+173798(X(i)*5)-251552(X(i)4)+171250(X(i)3)-55740(X(i)**2)+8476.7X(i)-150.38 TG=(-526.02)(X(i)*6)+2020.3(X(i)*5)-2690(X(i)*4)+1349.7(X(i)3)-138.25(x(i)**2)+9.5664X(i)+294.01

END IF

IF (i.EQ.1) THEN TP=TP(i) END IF IF (i.GT.1) THEN TP(i)=TP(i-1)+(TS-TP(i-1))DELTAX(i)0.001(NUS/DH)(5.73+0.06638TP(i-1))OS/(MK(978.2+0.1178TP(i-1)))+(TG-TP(i))DELTAX(i)0.001(NUGW/DH)(5.73+0.06638TP(i-1))OG/(MK(978.2+0.1178TP(i-1)))

END IF

!!!!!!!!!!!!WYNIKI!!!!!!!!!!!!!! WRITE (,)'Wyniki obliczen' ,ALKGW WRITE (,)'TS',TS WRITE (,)'TG',TG WRITE (,)'TP',TP

  END PROGRAM

The code looks if it is out of order ! I've provided some more questions that might help you find what is going wrong.

  ALKS  = NUS*LAMBDAP/DH         ! why are these calculated ?
  ALKGW = NUGW*LAMBDAP/DH        ! why are these calculated ? 
  CP    = 978.2+0.1178*TP(i)     ! why are these calculated ?
! 
END DO   !  why does the do loop end here, as the following code referrs to 'i' 
!  as the loop ends here, i will have the value N+1. Is that what you want ?
!
!!!!!!!!!!!!!!!!!!!!!!!OBLICZENIA TEMPERATUR POW. SILNIKA I POW. GONDOLI!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
!
WRITE (*,*) 'Podaj nr seansu badawczego, S z zakresu 1-4:' 
READ (*,*) S    !  is this the same S calculated above ?  
!
!SEANS 1!   as the following TS and TG calculation refers to X(i), this should be in a DO loop
IF ( s.eq.1 ) then     !  this is a bad test as S is a real
 IF (N.LT.2) THEN 
  TS = polynomial ( (/ 320.0, -40.929, 475.13 /), X(i) ) 
  TG = polynomial ( (/ 294.01, 9.5664, -138.25, 1349.7, -2690.0, 2020.3, -526.02 /), X(i) ) 
 else
  TS = polynomial ( (/ -150.38, 8476.7, -55740.0, 171250.0, -251552.0, 173798.0, -45547.0 /), X(i) ) 
  TG = polynomial ( (/  294.01, 9.5664,  -138.25,   1349.7,   -2690.0,   2020.3,  -526.02 /), X(i) ) 
 END IF 
END IF 

!  what happens to TS and TG if S /= 1  ?

IF (i.EQ.1) TP = TP(i)      ! I presume this is to initialise all of TP to 293 
!
IF (i.GE.1) THEN            ! you need a special case for i=1. what is the value of TP(i-1) ? Is it 293 ?
  if (i> 1) tpm1 = TP(i-1)
  if (i<=1) tpm1 = 293      ! best estimate : need a value if i=1
!
  y  = DELTAX(i)*0.001/DH * (5.73+0.06638*tpm1)*OS/ (MK*(978.2+0.1178*tpm1)) 
  TP(i) = tpm1                  & 
        + (TS-tpm1 ) *NUS  * y  & 
        + (TG-TP(i)) *NUGW * y
END IF 

! This is the end of using 'i' so where the DO loop could finish 

!!!!!!!!!!!!WYNIKI!!!!!!!!!!!!!! 
WRITE (*,*)'Wyniki obliczen' ,ALKGW 
WRITE (*,*)'TS',TS 
WRITE (*,*)'TG',TG 
WRITE (*,*)'TP',TP 


END PROGRAM