1 | include(dom.inc)
2 |
3 | PROGRAM solinit_interp1D
4 |
5 | ! ================================================================!
6 | ! !
7 | ! solinit_interp1D.F : Creates an initial solution for PRISSMA !
8 | ! !
9 | ! author : J. AMAYA (october 2007) !
10 | ! !
11 | ! ================================================================!
12 |
13 | IMPLICIT NONE
14 |
15 | include 'dom_constants.h'
16 |
17 | DOM_INT :: ntcells, ntnodes, ntfaces, inode
18 | DOM_INT :: i, j
19 | DOM_INT :: npts, P1(1), P2(1), tmp(1)
20 |
21 | CHARACTER*80 :: nodesfile, clnodefile, propfile, file1
22 | CHARACTER*80 :: gdatafile, path, kextfile
23 |
24 | DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: X1
25 | DOM_REAL, ALLOCATABLE, DIMENSION(:) :: Xtmp
26 | DOM_REAL :: interpol(5)
27 |
28 | DOM_REAL :: T_n, P_n, XH2O_n, XCO2_n
29 | DOM_REAL :: XCO_n, XO2_n, XN2_n
30 | DOM_REAL :: XSOOT_n
31 | DOM_REAL :: kscattering_n
32 | DOM_REAL :: kabsorption_n
33 | DOM_REAL, ALLOCATABLE, DIMENSION(:) :: x, y, z
34 |
35 | ! ------------------!
36 | ! Read choices file !
37 | ! ------------------!
38 |
39 | OPEN (FILE_CHCS , FILE='solinit_interp1D.choices', &
40 | & FORM='FORMATTED')
41 | READ (FILE_CHCS,*) file1
42 | READ (FILE_CHCS,*) npts
43 | READ (FILE_CHCS,*) path
44 | READ (FILE_CHCS,*) P_n
45 | CLOSE(FILE_CHCS)
46 |
47 | ALLOCATE(X1(5,npts))
48 |
49 | ! ------------------!
50 | ! Read profile file !
51 | ! ------------------!
52 |
53 | open(UNIT=2,file=file1,form='FORMATTED')
54 | READ(2,*)
55 |
56 | PRINT*, "Line X T/4000 ",&
57 | & " H2O CO2"
58 |
59 | DO i=1, npts
60 | READ(2,*) (X1(j,i),j=1,5)
61 | PRINT*,i,X1(1:4,i)
62 | ENDDO
63 |
64 | close(2)
65 |
66 | ! -----------!
67 | ! Open files !
68 | ! -----------!
69 |
70 | nodesfile = path(1:len_trim(path))//'/Nodelist.in'
71 | clnodefile = path(1:len_trim(path))//'/Cellnodes.in'
72 | propfile = path(1:len_trim(path))//'/Properties.in'
73 | gdatafile = path(1:len_trim(path))//'/Global.in'
74 | kextfile = path(1:len_trim(path))//'/K_Extinction.in'
75 |
76 | OPEN(FILE_NODES, FILE=nodesfile , FORM='UNFORMATTED')
77 | OPEN(FILE_CLNOD, FILE=clnodefile, FORM='UNFORMATTED')
78 | OPEN(FILE_GDATA, FILE=gdatafile , FORM='UNFORMATTED')
79 | OPEN(FILE_PROP , FILE=propfile , FORM='UNFORMATTED')
80 | OPEN(FILE_KEXT , FILE=kextfile , FORM='UNFORMATTED')
81 |
82 | ! ----------------------!
83 | ! Read Global mesh data !
84 | ! ----------------------!
85 |
86 | READ(FILE_GDATA) ntcells, ntnodes, ntfaces
87 | CLOSE(FILE_GDATA)
88 |
89 | ! -----------------!
90 | ! Allocate vectors !
91 | ! -----------------!
92 |
93 | ALLOCATE(x (ntnodes))
94 | ALLOCATE(y (ntnodes))
95 | ALLOCATE(z (ntnodes))
96 |
97 | ! --------------------!
98 | ! Read node positions !
99 | ! --------------------!
100 |
101 | ALLOCATE(Xtmp(npts))
102 |
103 | DO i=1, ntnodes
104 | READ(FILE_NODES) inode, x(i), y(i), z(i)
105 |
106 | Xtmp = abs(X1(1,:)-x(i))
107 | P1=MINLOC(Xtmp)
108 | Xtmp(P1)=1e10
109 | P2=MINLOC(Xtmp)
110 |
111 | IF(P1(1).gt.P2(1)) THEN
112 | tmp = P2
113 | P2 = P1
114 | P1 = tmp
115 | ENDIF
116 |
117 | interpol(:) = X1(:,P1(1)) + (x(i)-X1(1,P1(1)))* &
118 | & (X1(:,P2(1))-X1(:,P1(1)))/ &
119 | & (X1(1,P2(1))-X1(1,P1(1)))
120 |
121 | T_n = interpol(2)*4000.d0
122 | XH2O_n = interpol(3)
123 | XCO2_n = interpol(4)
124 | XCO_n = 0.d0
125 | XO2_n = 0.d0
126 | XN2_n = 0.d0
127 | XSOOT_n = 0.d0
128 |
129 | kscattering_n = 0.
130 | kabsorption_n = 0.
131 |
132 | WRITE(FILE_PROP ) i, T_n, P_n, XH2O_n, XCO2_n, XCO_n, XO2_n, &
133 | & XN2_n, XSOOT_n
134 | WRITE(FILE_KEXT ) i, kabsorption_n, kscattering_n
135 |
136 | ENDDO
137 |
138 | CLOSE(FILE_NODES)
139 | CLOSE(FILE_CLNOD)
140 | CLOSE(FILE_PROP )
141 | CLOSE(FILE_KEXT )
142 |
143 | WRITE(*,*) " Initial solution file successfully written"
144 |
145 | END PROGRAM solinit_interp1D