ray.F [SRC] [CPP] [JOB] [SCAN]
TOOLS / RAY / SRC



   1 | include(dom.inc)
   2 | 
   3 |       PROGRAM ray
   4 | 
   5 | !     ==================================================================!
   6 | !                                                                       !
   7 | !     ray.F   : Program that calculates the luminance that exits a      !
   8 | !               1D gaz column of length "e" using different spectral    !
   9 | !               methods                                                 !
  10 | !                                                                       !
  11 | !     author  : J.AMAYA (April 2008)                                    !
  12 | !                                                                       !
  13 | !     ==================================================================!
  14 | 
  15 |         IMPLICIT NONE
  16 | 
  17 |         include 'dom_constants.h'
  18 | 
  19 |         DOM_REAL      :: e,t_start,t_stop,time
  20 |         DOM_REAL      :: Tw, Ka0, I0, L0, Ltot, Lw, attenuation
  21 |         DOM_REAL      :: ls, kabsmoy, Twallout, Lwallout
  22 |         DOM_REAL      :: Kn, L, Lq, somme, Lnus, Qw
  23 |         DOM_REAL      :: blae, x, planck
  24 |         DOM_INT       :: n_segments, ngbands, j, s, i, k, n
  25 |         DOM_INT       :: nbands, nquad, nintbands
  26 |         CHARACTER*3   :: homosyst
  27 |         CHARACTER*15  :: mediumtype
  28 | 
  29 |         DOM_REAL, allocatable, dimension(:,:)   :: celldata,tau
  30 |         DOM_REAL, allocatable, dimension(:)     :: tempe, Lb, press
  31 |         DOM_REAL, allocatable, dimension(:)     :: XH2O, XCO2, XCO
  32 |         DOM_REAL, allocatable, dimension(:)     :: XO2, XN2, XSOOT
  33 |         DOM_REAL, allocatable, dimension(:)     :: all_WVNB, all_DWVNB
  34 |         DOM_REAL, allocatable, dimension(:)     :: deltanu
  35 |         DOM_REAL, allocatable, dimension(:)     :: wquad
  36 |         DOM_REAL, allocatable, dimension(:,:,:) :: kabs
  37 |         DOM_REAL, allocatable, dimension(:,:)   :: KCO,DCO,KC,DC,KH,DH
  38 | 
  39 | !       Tabulation variables
  40 |         DOM_INT                                    :: nYH, nYC, nYCO, nT
  41 |         DOM_REAL                                   :: DYH, DYC, DYCO, DT
  42 |         DOM_REAL                                   :: MH2O, MCO2, MCO
  43 |         DOM_REAL, ALLOCATABLE,DIMENSION(:,:,:,:,:) :: tabkabs
  44 | 
  45 |         DOM_INT        ::cnttemp
  46 |         character*64   ::outfile,tabfile
  47 |         character*80   ::pathspec
  48 | 
  49 | !       ------------------!
  50 | !       Read choices file !
  51 | !       ------------------!
  52 | 
  53 |         call CPU_TIME(t_start)
  54 | 
  55 |         WRITE(*,*) " >> Reading choices file"
  56 |         OPEN (FILE_CHCS , FILE='ray.choices', FORM='FORMATTED')
  57 |         READ(FILE_CHCS,*) pathspec
  58 |         READ(FILE_CHCS,*) mediumtype
  59 |         READ(FILE_CHCS,*) homosyst
  60 |         READ(FILE_CHCS,*) e
  61 |         READ(FILE_CHCS,*) n_segments
  62 |         READ(FILE_CHCS,*) nquad
  63 |         READ(FILE_CHCS,*) Tw
  64 |         READ(FILE_CHCS,*) Twallout
  65 |         CLOSE(FILE_CHCS)
  66 | 
  67 |       IF (mediumtype.eq.'TAB') THEN
  68 |         deltanu = 1.
  69 |         tabfile = pathspec(1:len_trim(pathspec))//'/table.dat'
  70 |         PRINT*, "Tabulated File : ",tabfile
  71 | 
  72 |         OPEN(12,FILE=tabfile,FORM='FORMATTED')
  73 |         READ(12,*) nquad
  74 |         READ(12,*) DT
  75 |         READ(12,*) DYH, MH2O
  76 |         READ(12,*) DYC, MCO2
  77 |         READ(12,*) DYCO, MCO
  78 |         CLOSE(12)
  79 | 
  80 |         nT=INT((2900-300)/DT)+1
  81 |         nYH=INT((MH2O-0.0)/DYH)+1
  82 |         nYC=INT((MCO2-0.0)/DYC)+1
  83 |         nYCO=INT((MCO-0.0)/DYCO)+1
  84 |         ALLOCATE(tabkabs(nYH,nYC,nYCO,nT,nquad))
  85 |         CALL readtab(tabfile,tabkabs,nYH,nYC,nYCO,nT,nquad)
  86 |         PRINT*, "Table for FS-SNBcK model readed"
  87 |       ENDIF
  88 | 
  89 | !       ---------------------!
  90 | !       Initialize variables !
  91 | !       ---------------------!
  92 | 
  93 |         WRITE(*,*) " >> Initializing data and allocating memory"
  94 |         ls         = e / n_segments
  95 |         nbands     = 371
  96 | 
  97 |         IF (ALLOCATED(celldata)) DEALLOCATE(celldata)
  98 |         IF (ALLOCATED(Lb)) DEALLOCATE(Lb)
  99 |         IF (ALLOCATED(press)) DEALLOCATE(press)
 100 |         IF (ALLOCATED(XH2O)) DEALLOCATE(XH2O)
 101 |         IF (ALLOCATED(XCO)) DEALLOCATE(XCO)
 102 |         IF (ALLOCATED(XCO2)) DEALLOCATE(XCO2)
 103 |         IF (ALLOCATED(XO2)) DEALLOCATE(XO2)
 104 |         IF (ALLOCATED(XN2)) DEALLOCATE(XN2)
 105 |         IF (ALLOCATED(XSOOT)) DEALLOCATE(XSOOT)
 106 |         IF (ALLOCATED(tempe)) DEALLOCATE(tempe)
 107 |         IF (ALLOCATED(all_WVNB)) DEALLOCATE(all_WVNB)
 108 |         IF (ALLOCATED(all_DWVNB)) DEALLOCATE(all_DWVNB)
 109 |         IF (ALLOCATED(KCO)) DEALLOCATE(KCO)
 110 |         IF (ALLOCATED(DCO )) DEALLOCATE(DCO)
 111 |         IF (ALLOCATED(KC)) DEALLOCATE(KC)
 112 |         IF (ALLOCATED(DC )) DEALLOCATE(DC)
 113 |         IF (ALLOCATED(KH)) DEALLOCATE(KH)
 114 |         IF (ALLOCATED(DH)) DEALLOCATE(DH)
 115 |         IF (ALLOCATED(wquad)) DEALLOCATE(wquad)
 116 | 
 117 |         allocate(celldata (n_segments,8))
 118 |         allocate(Lb       (n_segments))
 119 |         allocate(press    (n_segments))
 120 |         allocate(XH2O     (n_segments))
 121 |         allocate(XCO2     (n_segments))
 122 |         allocate(XCO      (n_segments))
 123 |         allocate(XO2      (n_segments))
 124 |         allocate(XN2      (n_segments))
 125 |         allocate(XSOOT    (n_segments))
 126 |         allocate(tempe    (n_segments))
 127 |         allocate(all_WVNB (nbands))
 128 |         allocate(all_DWVNB(nbands))
 129 |         allocate(KCO      (14,nbands))
 130 |         allocate(DCO      (14,nbands))
 131 |         allocate(KC       (14,nbands))
 132 |         allocate(DC       (14,nbands))
 133 |         allocate(KH       (14,nbands))
 134 |         allocate(DH       (14,nbands))
 135 |         allocate(wquad    (nquad))
 136 | 
 137 | !       --------------------!
 138 | !       Temperature profile !
 139 | !       --------------------!
 140 | 
 141 |         DO i = 1, n_segments
 142 |           x = i*ls
 143 |           include 'profile.inc'
 144 |           print*, i,":", tempe(i)
 145 |         ENDDO
 146 | 
 147 |         celldata(:,1) = tempe(:)
 148 |         celldata(:,2) = press(:)
 149 |         celldata(:,3) = XH2O(:)
 150 |         celldata(:,4) = XCO2(:)
 151 |         celldata(:,5) = XCO(:)
 152 |         celldata(:,6) = XO2(:)
 153 |         celldata(:,7) = XN2(:)
 154 |         celldata(:,8) = XSOOT(:)
 155 | 
 156 |         deallocate(press, XH2O, XCO2, XCO, XO2, XN2, XSOOT)
 157 | 
 158 | !       ------------------------------!
 159 | !       Read data from spectral files !
 160 | !       ------------------------------!
 161 | 
 162 |         WRITE(*,*) " >> Reading spectral data"
 163 |         OPEN(UNIT=49,FILE='SPECTRAL/SNBWN')
 164 | 
 165 |         i = 1
 166 |         DO WHILE (i.le. nbands)
 167 | 
 168 |           READ(49,*) all_WVNB(i),all_DWVNB(i)
 169 |           IF (all_WVNB(i).lt.0.) i = nbands
 170 |           i = i + 1
 171 | 
 172 |         ENDDO
 173 | 
 174 |         CLOSE(49)
 175 | 
 176 |         call PARAM(KCO,KC,KH,DCO,DC,DH,pathspec)
 177 | 
 178 | !       ------------------------!
 179 | !       Select sommation limits !
 180 | !       ------------------------!
 181 | 
 182 |         WRITE(*,*) " >> Selecting limits for the integration"
 183 | 
 184 |         IF (mediumtype.eq.'FSCK') THEN
 185 |           ngbands   = nbands
 186 |           nintbands = 1
 187 |         ELSEIF (mediumtype.eq.'SNB-CK') THEN
 188 |           ngbands    = 1
 189 |           nintbands  = nbands
 190 |         ELSEIF (mediumtype.eq.'MALKMUS') THEN
 191 |           ngbands    = 1
 192 |           nintbands  = nbands
 193 |         ELSEIF (mediumtype.eq.'TAB') THEN
 194 |           ngbands   = nbands
 195 |           nintbands = 1
 196 |         ENDIF
 197 | 
 198 |         IF (ALLOCATED(kabs))    DEALLOCATE(kabs)
 199 |         IF (ALLOCATED(deltanu)) DEALLOCATE(deltanu)
 200 | 
 201 |         allocate(kabs   (n_segments, nquad, nintbands))
 202 |         allocate(deltanu(nintbands))
 203 | 
 204 |         IF (mediumtype.eq.'FSCK') THEN
 205 |           deltanu = 1.
 206 |         ELSEIF (mediumtype.eq.'SNB-CK') THEN
 207 |           deltanu(:) = all_DWVNB(:) * 100.
 208 |         ELSEIF (mediumtype.eq.'MALKMUS') THEN
 209 |           deltanu(:) = all_DWVNB(:) * 100.
 210 |         ELSE IF (mediumtype.eq.'TAB') THEN
 211 |           deltanu = 1.
 212 |         ENDIF
 213 | 
 214 | !       ---------------------!
 215 | !       Printing information !
 216 | !       ---------------------!
 217 | 
 218 |         WRITE(*,*) " __// RAY v1.0 \\__ "
 219 |         WRITE(*,*)
 220 |         WRITE(*,*) "   Information:"
 221 |         WRITE(*,*) "   ------------"
 222 |         WRITE(*,*)
 223 |         WRITE(*,*) "   Spectral model       :", mediumtype
 224 |         WRITE(*,*) "   Number of segments   :", n_segments
 225 |         WRITE(*,*) "   Quadrature points    :", nquad
 226 |         WRITE(*,*) "   Total number of bands:", nbands
 227 |         WRITE(*,*) "   Total grouped bands  :", ngbands
 228 |         WRITE(*,*) "   Total integrations   :", nintbands
 229 |         WRITE(*,*) "   Gaz wall thickness   :", e
 230 |         WRITE(*,*) "   Gaz segment thickness:", ls
 231 |         WRITE(*,*) "   Wall temperature     :", Tw
 232 |         WRITE(*,*)
 233 | 
 234 |         IF (nquad.ge.10) THEN
 235 |           outfile=trim(mediumtype)//'_Nq='//CHAR(48+int(nquad/10))      &
 236 |      &           //CHAR(48+modulo(nquad,10))//'_out.dat'
 237 |         ELSE
 238 |           outfile = trim(mediumtype)//'_Nq='//CHAR(48+nquad)//'_out.dat'
 239 |         ENDIF
 240 | 
 241 |         OPEN (UNIT=102 , FILE=outfile, FORM='FORMATTED')
 242 | 
 243 | !        ------------------------!
 244 | !        Calculate spectral data !
 245 | !        ------------------------!
 246 | 
 247 |         WRITE(*,*) " >> Calculate absorption coefficients"
 248 | 
 249 |         IF (mediumtype.eq.'FSCK') THEN
 250 |           Lb(:) = planck(celldata(:,1))
 251 |           call FSCK_CASE(kabs(:,:,1), wquad, Lb, celldata, n_segments,  &
 252 |      &                    all_WVNB, all_DWVNB, KCO,DCO,KC,DC,KH,DH,     &
 253 |      &                      homosyst, ngbands, nbands, nquad)
 254 | 
 255 |         ELSEIF(mediumtype.eq.'TAB') THEN
 256 |           Lb(:) = planck(celldata(:,1))
 257 |           CALL TAB_CASE(celldata,kabs(:,:,1), wquad, tabkabs, DYH,      &
 258 |      &                  DYC, DYCO, nYH, nYC, nYCO, DT, nT, nquad,       &
 259 |      &                  n_segments, nbands, all_WVNB, all_DWVNB, Lb)
 260 | 
 261 |         ELSEIF (mediumtype.eq.'SNB-CK') THEN
 262 |           DO i = 1, nintbands
 263 | !           WRITE(*,*) " Absorption coeff - Doing band: ",i,"/",nbands
 264 |             Lb(:) = blae(all_WVNB(i)*100.,celldata(:,1))/pi
 265 |             call SNB_CASE (kabs(:,:,i), wquad, Lb, celldata,            &
 266 |      &                     n_segments, all_WVNB, all_DWVNB, KCO, DCO,   &
 267 |      &                     KC, DC, KH, DH, homosyst,ngbands,i, nbands,  &
 268 |      &                     nquad)
 269 |           ENDDO
 270 | 
 271 |         ELSEIF (mediumtype.eq.'MALKMUS') THEN
 272 |           IF (ALLOCATED(tau)) DEALLOCATE(tau)
 273 |           ALLOCATE(tau(n_segments,nbands))
 274 |           DO i = 1, nintbands
 275 |             call MALKMUS_CASE(tau(:,i),ls,celldata,n_segments,          &
 276 |      &           all_WVNB, all_DWVNB,KCO,DCO,KC,DC,KH,DH,               &
 277 |      &           i,nbands)
 278 |           ENDDO
 279 | 
 280 |         ENDIF
 281 | 
 282 | !       --------------------------------!
 283 | !       Integrate over all wave lengths !
 284 | !       --------------------------------!
 285 | 
 286 |         L = 0.
 287 |         k = 0
 288 |         kabsmoy = 0.
 289 | 
 290 |         WRITE(*,*) " >> Luminance calculation"
 291 |         DO i = 1, nintbands
 292 | 
 293 | !         -------------------------------!
 294 | !         Integrate over all quadratures !
 295 | !         -------------------------------!
 296 | 
 297 |           WRITE(*,*) " Narrow band ",i,"/",nintbands
 298 |           Lq = 0.
 299 | 
 300 |           IF (mediumtype.eq.'MALKMUS') THEN
 301 | 
 302 |             somme = 0.
 303 |             DO s = 1, n_segments
 304 | !             print*, "     Segment: ", s,"/", n_segments
 305 | 
 306 |               Kn = 1.
 307 |               DO n = s+1, n_segments
 308 |                 Kn = Kn * tau (n,i)
 309 |               ENDDO
 310 | 
 311 |               Lnus  = blae(all_WVNB(i)*100.,celldata(s,1))/pi
 312 |               somme = somme + (1-tau(s,i))*Lnus*Kn
 313 | !             print*, "     - tau(s,i)   :", tau(s,i)
 314 | !             print*, "     - Kn, Lnus   :", Kn, Lnus
 315 | 
 316 |             ENDDO
 317 | 
 318 |             L  = L + deltanu(i) * somme
 319 | !           print*, " - L, deltanu, Lq : ", L, deltanu(i), Lq
 320 | 
 321 |           ELSE
 322 | 
 323 |             DO j = 1, nquad
 324 | 
 325 | !             ----------------------------------------------!
 326 | !             Emission for each segment of the optical path !
 327 | !             ----------------------------------------------!
 328 | !             print*, "   Quadrature point: ", j,"/", nquad
 329 | 
 330 |               somme = 0.
 331 |               DO s = 1, n_segments
 332 | !               print*, "     Segment: ", s,"/", n_segments
 333 | 
 334 |                 Kn = 0.
 335 |                 DO n = s+1, n_segments
 336 |                   Kn = Kn + kabs(n,j,i)*ls
 337 | !                 print*, "       Next segment ", n,"/",n_segments
 338 | !                 print*, "       - kabs(n,j,i):", kabs(n,j,i)
 339 | !                 print*, "       - ls         :", ls
 340 |                 ENDDO
 341 | 
 342 |                 IF (mediumtype.eq.'SNB-CK') THEN
 343 |                   Lnus  = blae(all_WVNB(i)*100.,celldata(s,1))/pi
 344 |                 ELSE IF (mediumtype.eq.'FSCK') THEN
 345 |                   Lnus  = planck(celldata(s,1))
 346 |                 ENDIF
 347 | 
 348 |                 somme = somme + (1-exp(-kabs(s,j,i)*ls))*Lnus*exp(-Kn)
 349 | 
 350 | !               print*, "     - kabs(s,j,i):",kabs(s,j,i)
 351 | !               print*, "     - Kn, Lnus   :", Kn, Lnus
 352 |               ENDDO
 353 | 
 354 | !             ---------------------------------------!
 355 | !             Calculate absorption of wall intensity !
 356 | !             ---------------------------------------!
 357 | 
 358 |               I0   = blae(all_WVNB(i)*100.,Tw)/pi
 359 |               Ka0  = 0.
 360 | 
 361 |               DO s = 1, n_segments
 362 |                 Ka0 = Ka0 + kabs(s,j,i)*ls
 363 |               ENDDO
 364 | 
 365 |               L0 = I0 * exp( -Ka0 ) * wquad(j)
 366 | 
 367 | !             -----------------------------------------!
 368 | !             integrate this spectral quadrature point !
 369 | !             -----------------------------------------!
 370 | 
 371 |               Lq      = Lq + somme * wquad(j)
 372 |               kabsmoy = kabsmoy + wquad(j) * deltanu(i) * kabs(1,j,i)   &
 373 |      &                  * blae(all_WVNB(i)*100.,celldata(1,1))/pi
 374 | 
 375 |               IF (mediumtype.eq.'FSCK') k = k + 1
 376 | 
 377 |             ENDDO
 378 | 
 379 |             L  = L  + deltanu(i) * Lq
 380 |             Lw = Lw + deltanu(i) * L0
 381 | 
 382 | !           print*, " - L, deltanu, Lq : ", L, deltanu(i), Lq
 383 | 
 384 |             IF (mediumtype.eq.'SNB-CK') k = k + 1
 385 | 
 386 |           ENDIF
 387 | 
 388 |         ENDDO
 389 | 
 390 |         kabsmoy     = kabsmoy / planck(celldata(1,1))
 391 |         Ltot        = Lw + L
 392 |         attenuation = Ltot / planck(Tw)
 393 |         Lwallout    = planck(Twallout)
 394 |         Qw          = Lwallout - Ltot
 395 | 
 396 |         WRITE(*,*) " Total wall luminance reaching the exit, Lw =",Lw
 397 |         WRITE(*,*) " Emited/self-absorbed luminance of the gas, L =",L
 398 |         WRITE(*,*) " Luminance at the exit, Ltot = Lw + L = ", Ltot
 399 |         WRITE(*,*) " Out wall total heat flux, Qw =", Qw
 400 | 
 401 |         WRITE(*,*)
 402 |         WRITE(*,*) " kabsmoy = ", kabsmoy
 403 |         WRITE(*,*)
 404 | 
 405 |         WRITE(102,*) celldata(n_segments,1), Lw, L, Ltot,               &
 406 |      &               attenuation, kabsmoy, Qw
 407 | 
 408 | 
 409 |         CLOSE(102)
 410 |         call CPU_TIME(t_stop)
 411 | 
 412 |         PRINT*,'Elapsed CPU time = ',t_stop-t_start
 413 | 
 414 |         outfile="TIME_"//outfile
 415 |         OPEN(103,file=outfile,form='formatted')
 416 |         WRITE(103,*) t_stop-t_start
 417 |         CLOSE(103)
 418 | 
 419 |       END PROGRAM ray


ray.F could be called by:
Makefile [TOOLS/RAY] - 50 - 56 - 80 - 112
ray.F [TOOLS/RAY/SRC] - 3 - 218
ray [TOOLS/SCRIPTS] - 43 - 44 - 45 - 55