aboutsummaryrefslogtreecommitdiff
path: root/orrs/src/fsrun.f
blob: 5b55c5666f7f03ad6bfe13c8d6be166dfb199b9f (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
      PROGRAM FSRUN
      PARAMETER (NMAX=256)
      DIMENSION ETA(NMAX), F(NMAX), U(NMAX), S(NMAX)
      CHARACTER*1 ANS
C
      LST = 1
      LRE = 1
C
      N = 256
      ETAE = 30.0
      GEO = 1.01
C
      IDEV = 6
      SIZE = 6.0
      IHARD = -999
C
      EWT = 1.0/ETAE
      UWT = 0.5
      CH = 0.02
C
      CALL PLOTS(0,IHARD,IDEV)
      CALL FACTOR(SIZE)
C
      CALL PLOT(0.7,0.1,-3)
C
      CALL NEWPEN(1)
C
      CALL PLOT(0.0,0.0,3)
      CALL PLOT(UWT*1.0,0.0,2)
      CALL PLOT(0.0,0.0,3)
      CALL PLOT(0.0,EWT*ETAE,2)
C
      WRITE(*,*) 'Enter H1, H2, dH'
      READ (*,*) H1,H2,DH
C
      NH = INT((H2-H1)/DH) + 1
C
      open(7,file='hfuns.fs',status='unknown')
c
      CALL NEWPEN(3)
      DO 10 IH=1, NH
        H = H1 + DH*FLOAT(IH-1)
        CALL FS(3,2,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
C---------------------
c        BU = H
c        CALL FS(1,1,BU,H,N,ETAE,GEO,ETA,F,U,S,DELTA)
C---------------------
C
        DSI = 0.0
        THI = 0.0
        TSI = 0.0
        CDN = 0.0
        DO 103 I=1, N-1
          UA = 0.5*(U(I+1) + U(I))
          DETA = ETA(I+1) - ETA(I)
C
          DSI = DSI + (1.0 -    UA)   *DETA
          THI = THI + (1.0 -    UA)*UA*DETA
          TSI = TSI + (1.0 - UA*UA)*UA*DETA
C
          CDN = CDN + (U(I+1) - U(I))**2 / DETA
 103    CONTINUE
C
        HK = DSI/THI
        HS = TSI/THI
C
        CDN = CDN *THI * 2.0/HS
        CFN = S(1)*THI
C
        DSI = DSI*DELTA        
        THI = THI*DELTA        
        TSI = TSI*DELTA        
C
        BUF = (CFN - CDN)/(HK-1.0) / THI**2
        write(*,*) H, BU, THI**2, 0.5*(BU + 1.0) * THI**2
        write(7,*) H, BU, THI**2, 0.5*(BU + 1.0) * THI**2
c
        CALL PLOTON
        CALL PLOT(UWT*U(1),EWT*ETA(1),3)
        DO 105 I=2, N
          CALL PLOT(UWT*U(I),EWT*ETA(I),2)
 105    CONTINUE
        CALL PLOTOF
 10   CONTINUE
c
      close(7)
C
c      WRITE(6,*) 'Hit <cr>'
c      READ (5,8000) ANS
c 8000 FORMAT(A1)
C
      CALL PLOT(0.0,0.0,+999)
      STOP
      END