1 | include(dom.inc)
   2 | 
   3 |       PROGRAM dom2ascii
   4 | 
   5 | !       ================================================================!
   6 | !                                                                       !
   7 | !       dom2ascii.F : Converts PRISSMA *.out files into ascii files     !
   8 | !                     and interpolates the data over a line             !
   9 | !                                                                       !
  10 | !       author      : J. AMAYA (mars 2007)                              !
  11 | !                                                                       !
  12 | !       ================================================================!
  13 | 
  14 |         IMPLICIT NONE
  15 | 
  16 |         include 'dom_constants.h'
  17 | 
  18 |         DOM_INT       :: npoints, i, j, l, typeunit, readwall, readkp
  19 |         DOM_INT       :: ntcells, ntnodes, ntfaces, ndirs, nbdyfaces
  20 |         DOM_INT       :: nmaxfaces
  21 |         DOM_INT       :: ios, ic, nf, iw, jw
  22 | 
  23 |         DOM_REAL      :: x1, y1, z1, x2, y2, z2
  24 |         DOM_REAL      :: xd, yd, zd
  25 |         DOM_REAL      :: segmentsize, d,unity
  26 | 
  27 |         DOM_REAL, ALLOCATABLE, DIMENSION(:)   :: x,  y,  z
  28 |         DOM_REAL, ALLOCATABLE, DIMENSION(:)   :: xn, yn, zn
  29 |         DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: xw, yw, zw
  30 |         DOM_REAL, ALLOCATABLE, DIMENSION(:)   :: xfc,yfc,zfc
  31 |         DOM_REAL, ALLOCATABLE, DIMENSION(:)   :: Sr,  H,  G,  Qw, Kp
  32 |         DOM_REAL, ALLOCATABLE, DIMENSION(:)   :: Srn, Hn, Gn, Qwn, Kpn
  33 |         DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: Qr
  34 |         DOM_REAL, ALLOCATABLE, DIMENSION(:,:) :: Qrn
  35 | 
  36 |         CHARACTER*80  :: infilespath, outfilespath
  37 |         CHARACTER*80  :: GFile, SrFile, HFile, QwFile, QrFile, Kpfile
  38 |         CHARACTER*80  :: clfacefile, gdatafile,ccellsfile,cfacesfile
  39 |         CHARACTER*80 :: gfile2,srfile2,hfile2,qwfile2,qrfile2,kpfile2
  40 | 
  41 | !       ------------------!
  42 | !       Read choices file !
  43 | !       ------------------!
  44 | 
  45 |         OPEN (FILE_CHCS , FILE='dom2ascii.choices', FORM='FORMATTED')
  46 |         READ (FILE_CHCS,*) infilespath
  47 |         READ (FILE_CHCS,*) outfilespath
  48 |         READ (FILE_CHCS,*) x1, y1, z1
  49 |         READ (FILE_CHCS,*) x2, y2, z2
  50 |         READ (FILE_CHCS,*) npoints
  51 |         READ (FILE_CHCS,*) typeunit
  52 |         READ (FILE_CHCS,*) readwall
  53 |         CLOSE(FILE_CHCS)
  54 | 
  55 |         IF (typeunit.eq.1) THEN
  56 |           unity=1.
  57 |         ELSE IF (typeunit.eq.2) THEN
  58 |           unity=1000.
  59 |         ENDIF
  60 | 
  61 |         IF (npoints.lt.2) THEN
  62 |           WRITE(*,*) " Fatal error: the minimum number of points is 2"
  63 |           STOP
  64 |         ENDIF
  65 | 
  66 | !       -----------------!
  67 | !       Allocate vectors !
  68 | !       -----------------!
  69 | 
  70 |         ALLOCATE(x(npoints))
  71 |         ALLOCATE(y(npoints))
  72 |         ALLOCATE(z(npoints))
  73 | 
  74 | !       ---------------------------!
  75 | !       Calculate geometrical data !
  76 | !       ---------------------------!
  77 | 
  78 |         d  = sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2))
  79 | 
  80 |         segmentsize = d / real(npoints - 1)
  81 |         xd          = (x2-x1) / d
  82 |         yd          = (y2-y1) / d
  83 |         zd          = (z2-z1) / d
  84 | 
  85 |         DO i=1, npoints
  86 |           x(i) = x1 + (i-1) * segmentsize * xd
  87 |           y(i) = y1 + (i-1) * segmentsize * yd
  88 |           z(i) = z1 + (i-1) * segmentsize * zd
  89 | !         print*, x(i), y(i), z(i)
  90 |         ENDDO
  91 | 
  92 | !       --------------------!
  93 | !       Open PRISSMA files !
  94 | !       --------------------!
  95 | 
  96 |         GFile  = trim(outfilespath)//'/G.out'
  97 |         SrFile = trim(outfilespath)//'/Sr.out'
  98 |         HFile  = trim(outfilespath)//'/H.out'
  99 |         QwFile = trim(outfilespath)//'/Qw.out'
 100 |         QrFile = trim(outfilespath)//'/Qr.out'
 101 |         KpFile = trim(outfilespath)//'/Kp.out'
 102 | 
 103 |         gdatafile  = trim(infilespath)//'/Global.in'
 104 |         clfacefile = trim(infilespath)//'/CLFaces.in'
 105 |         ccellsfile = trim(infilespath)//'/Centercells.in'
 106 |         cfacesfile = trim(infilespath)//'/Centerfaces.in'
 107 | 
 108 |         OPEN(FILE_G ,FILE=Gfile ,FORM='UNFORMATTED')
 109 |         OPEN(FILE_Sr,FILE=SrFile,FORM='UNFORMATTED')
 110 |         OPEN(FILE_H ,FILE=HFile ,FORM='UNFORMATTED')
 111 |         OPEN(FILE_Qw,FILE=QwFile,FORM='UNFORMATTED')
 112 |         OPEN(FILE_Qr,FILE=QrFile,FORM='UNFORMATTED')
 113 | 
 114 |         readkp = 0
 115 |         OPEN(FILE_Kp, FILE = KpFile, IOSTAT = ios, STATUS = 'old',      &
 116 |      &       FORM='UNFORMATTED')
 117 |         IF (ios.eq.0) readkp = 1
 118 | 
 119 |         OPEN(FILE_CLFAC, FILE=clfacefile, FORM='UNFORMATTED')
 120 |         OPEN(FILE_GDATA, FILE=gdatafile , FORM='UNFORMATTED')
 121 |         OPEN(FILE_CCELL, FILE=ccellsfile, FORM='UNFORMATTED')
 122 |         OPEN(FILE_CFACE, FILE=cfacesfile, FORM='UNFORMATTED')
 123 | 
 124 | !       ----------------!
 125 | !       Get global data !
 126 | !       ----------------!
 127 | 
 128 |         READ(FILE_GDATA) ntcells, ntnodes, ntfaces,nmaxfaces,ndirs
 129 |         READ(FILE_CLFAC) nbdyfaces
 130 | 
 131 |         WRITE(*,*) "Dirs : ", ndirs
 132 |         WRITE(*,*) "Cells: ", ntcells
 133 |         WRITE(*,*) "Nodes: ", ntnodes
 134 |         WRITE(*,*) "Faces: ", ntfaces
 135 |         WRITE(*,*) "Boundary faces: ", nbdyfaces
 136 |         WRITE(*,*)
 137 | 
 138 |         CLOSE(FILE_GDATA)
 139 |         CLOSE(FILE_CLFAC)
 140 | 
 141 | !       -----------------!
 142 | !       Allocate vectors !
 143 | !       -----------------!
 144 | 
 145 |         ALLOCATE(xn (ntcells))
 146 |         ALLOCATE(yn (ntcells))
 147 |         ALLOCATE(zn (ntcells))
 148 | 
 149 |         ALLOCATE(xw (nmaxfaces,ntcells))
 150 |         ALLOCATE(yw (nmaxfaces,ntcells))
 151 |         ALLOCATE(zw (nmaxfaces,ntcells))
 152 | 
 153 |         ALLOCATE(xfc(nbdyfaces))
 154 |         ALLOCATE(yfc(nbdyfaces))
 155 |         ALLOCATE(zfc(nbdyfaces))
 156 | 
 157 |         ALLOCATE(Sr (ntcells))
 158 |         ALLOCATE(H  (ntcells))
 159 |         ALLOCATE(G  (ntcells))
 160 |         ALLOCATE(Qw (ntcells))
 161 |         ALLOCATE(Srn(ntcells))
 162 |         ALLOCATE(Hn (ntcells))
 163 |         ALLOCATE(Gn (ntcells))
 164 |         ALLOCATE(Qwn(ntcells))
 165 | 
 166 |         IF (readkp.eq.1) ALLOCATE(Kp (ntcells))
 167 |         IF (readkp.eq.1) ALLOCATE(Kpn(ntcells))
 168 | 
 169 |         ALLOCATE(Qr (3,ntcells))
 170 |         ALLOCATE(Qrn(3,ntcells))
 171 | 
 172 | !       -------------------!
 173 | !       Initialise vectors !
 174 | !       -------------------!
 175 | 
 176 |         H   = 0.
 177 |         G   = 0.
 178 |         Qw  = 0.
 179 |         Sr  = 0.
 180 |         Qr  = 0.
 181 |         Hn  = 0.
 182 |         Gn  = 0.
 183 |         Qwn = 0.
 184 |         Srn = 0.
 185 |         Qrn = 0.
 186 | 
 187 |         IF (readkp.eq.1) Kp  = 0.
 188 |         IF (readkp.eq.1) Kpn = 0.
 189 | 
 190 | !       ---------------------!
 191 | !       Read geometical data !
 192 | !       ---------------------!
 193 | 
 194 |         xw = 0
 195 |         yw = 0
 196 |         zw = 0
 197 | 
 198 |         DO i = 1, ntcells
 199 |           READ(FILE_CCELL) ic,xn(ic),yn(ic),zn(ic)
 200 |           READ(FILE_CFACE) ic,nf,(xw(j,ic),yw(j,ic),zw(j,ic),j=1,nf)
 201 |         ENDDO
 202 | 
 203 |         CLOSE(FILE_CCELL)
 204 |         CLOSE(FILE_CFACE)
 205 | 
 206 | !       ---------------!
 207 | !       Read wall data !
 208 | !       ---------------!
 209 | 
 210 |         IF (readwall.eq.1) THEN
 211 | 
 212 |           DO i=1,nbdyfaces
 213 | 
 214 |             READ(FILE_Qw) iw,jw,Qwn(i)
 215 |             READ(FILE_H ) iw,jw,Hn(i)
 216 | 
 217 |             xfc(i) = xw(iw,jw)
 218 |             yfc(i) = yw(iw,jw)
 219 |             zfc(i) = zw(iw,jw)
 220 | 
 221 |           ENDDO
 222 | 
 223 |           CLOSE(FILE_Qw)
 224 |           CLOSE(FILE_H )
 225 | 
 226 | !         ----------------------!
 227 | !         Interpolate wall data !
 228 | !         ----------------------!
 229 | 
 230 |           CALL INTERP3D(npoints,  x,  y,  z,  Qw,                       &
 231 |      &                  nbdyfaces,xfc,yfc,zfc,Qwn,ntcells)
 232 |           CALL INTERP3D(npoints,  x,  y,  z,  H,                        &
 233 |      &                  nbdyfaces,xfc,yfc,zfc,Hn,ntcells)
 234 | 
 235 |         ENDIF
 236 | 
 237 |         DEALLOCATE(xfc,yfc,zfc,xw,yw,zw)
 238 | 
 239 | !       ---------------!
 240 | !       Read cell data !
 241 | !       ---------------!
 242 | 
 243 |         DO i=1,ntcells
 244 | 
 245 |           READ(FILE_G ) ic,Gn(ic)
 246 |           READ(FILE_Sr) ic,Srn(ic)
 247 |           READ(FILE_Qr) ic,(Qrn(l,ic),l=1,3)
 248 |           IF (readkp.eq.1) THEN
 249 |             READ(FILE_Kp) ic,Kpn(i)
 250 |           ENDIF
 251 | 
 252 |         ENDDO
 253 | 
 254 |         CLOSE(FILE_G )
 255 |         CLOSE(FILE_Sr)
 256 |         CLOSE(FILE_Qr)
 257 |         CLOSE(FILE_Kp)
 258 | 
 259 | !       ----------------------!
 260 | !       Interpolate cell data !
 261 | !       ----------------------!
 262 | 
 263 |         CALL INTERP3D(npoints,x,y,z,G,ntcells,xn,yn,zn,Gn,ntcells)
 264 |         CALL INTERP3D(npoints,x,y,z,Sr,ntcells,xn,yn,zn,Srn,ntcells)
 265 |         CALL INTERP3D(npoints,x,y,z,Qr(1,:),ntcells,xn,yn,zn,Qrn(1,:),  &
 266 |      &                ntcells)
 267 |         CALL INTERP3D(npoints,x,y,z,Qr(2,:),ntcells,xn,yn,zn,Qrn(2,:),  &
 268 |      &                ntcells)
 269 |         CALL INTERP3D(npoints,x,y,z,Qr(3,:),ntcells,xn,yn,zn,Qrn(3,:),  &
 270 |      &                ntcells)
 271 |         IF (readkp.eq.1) THEN
 272 |           CALL INTERP3D(npoints,x,y,z,Kp,ntcells,xn,yn,zn,Kpn,ntcells)
 273 |         ENDIF
 274 | 
 275 | !       ------------------!
 276 | !       Write ascii files !
 277 | !       ------------------!
 278 | 
 279 |         gfile2     = trim(outfilespath)//'/G.dat'
 280 |         srfile2    = trim(outfilespath)//'/Sr.dat'
 281 |         hfile2     = trim(outfilespath)//'/H.dat'
 282 |         qwfile2    = trim(outfilespath)//'/Qw.dat'
 283 |         qrfile2    = trim(outfilespath)//'/Qr.dat'
 284 |         kpfile2    = trim(outfilespath)//'/Kp.dat'
 285 | 
 286 |         OPEN(FILE_G2 ,FILE=gfile2 ,FORM='formatted')
 287 |         OPEN(FILE_Sr2,FILE=srfile2,FORM='formatted')
 288 |         OPEN(FILE_H2 ,FILE=hfile2 ,FORM='formatted')
 289 |         OPEN(FILE_Qw2,FILE=qwfile2,FORM='formatted')
 290 |         OPEN(FILE_Qr2,FILE=qrfile2,FORM='formatted')
 291 |         IF (readkp.eq.1) THEN
 292 |           OPEN(FILE_Kp2,FILE=kpfile2,FORM='formatted')
 293 |         ENDIF
 294 | 
 295 |         DO i=1,npoints
 296 | 
 297 |           WRITE(FILE_H2 ,'(3f16.8,1f18.6)') x(i),y(i),z(i),H(i)/unity
 298 |           WRITE(FILE_Qw2,'(3f16.8,1f18.6)') x(i),y(i),z(i),Qw(i)/unity
 299 |           WRITE(FILE_G2 ,'(3f16.8,1f18.6)') x(i),y(i),z(i),G(i)/unity
 300 |           WRITE(FILE_Sr2,'(3f16.8,1f18.6)') x(i),y(i),z(i),Sr(i)/unity
 301 |           WRITE(FILE_Qr2,'(3f16.8,3f18.6)') x(i),y(i),z(i),             &
 302 |      &                                     (Qr(l,i)/unity,l=1,3)
 303 |           IF (readkp.eq.1) THEN
 304 |             WRITE(FILE_Kp2,'(3f16.8,1f18.6)') x(i),y(i),z(i),Kp(i)/unity
 305 |           ENDIF
 306 | 
 307 |         ENDDO
 308 | 
 309 |         CLOSE(FILE_H2 )
 310 |         CLOSE(FILE_G2 )
 311 |         CLOSE(FILE_Sr2)
 312 |         CLOSE(FILE_Qr2)
 313 |         CLOSE(FILE_Qw2)
 314 |         CLOSE(FILE_Kp2)
 315 | 
 316 |       END PROGRAM dom2ascii