1C 2C 3C Portions Copyright (C) 2001 Joseph H. Morrison 4C 5C This file is part of ISAAC. 6C 7C This program is distributed under the terms of the ISAAC Public Source 8C License. This program is distributed WITHOUT ANY WARRANTY; without 9C even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 10C PURPOSE. 11C 12C You should have received a copy of the ISAAC Public Source License 13C with this program. If you did not, you may get a copy of the license 14C at http://isaac-cfd.sourceforge.net. 15C 16C $Revision: 4.7 $ 17C $Author: jhmorr $ 18C $Date: 2001/10/29 03:25:28 $ 19C $State: Exp $ 20C $Log: main.F,v $ 21C Revision 4.7 2001/10/29 03:25:28 jhmorr 22C Updated license information 23C 24C Revision 4.6 2001/06/22 03:47:28 jhmorr 25C Add output of function name file. Modify contact information 26C 27C Revision 4.5 2001/06/20 03:28:43 jhmorr 28C Added test for STOP file to stop a run. 29C 30C Revision 4.4 2001/06/18 03:36:04 jhmorr 31C Added calculation of computational time. 32C 33C Revision 4.3 2001/06/08 06:25:36 jhmorr 34C Fixed error in BCPRSB where S, PROPS are dimensioned but not passed 35C 36C Revision 4.2 2001/06/08 04:56:15 jhmorr 37C Added notice 38C 39C Revision 4.1 1998/04/14 20:33:16 jhmorr 40C Alpha 4.1. 41C 42C 43 44 45 46 47 48 PROGRAM ISAAC 49C 50C Program: ISAAC 51C (Integrated Solution Algorithm for Arbitrary Configurations) 52C Author: Joseph H. Morrison 53C 54C Original code developed under contract to the NASA Langley Research Center 55C while the author worked for: 56C Analytical Services and Materials, Inc. 57C 107 Research Drive 58C Hampton, Virginia 23666 59C 60C 61C Main routine for ISAAC 62C 63C 64C Parameter statements setting up the maximum storage 65C MXBLKS : maximum number of blocks 66C MXLVLS : maximum number of grid levels 67C MXCFS : maximum number of coarse to fine sequences 68C MXBCS : maximum number of boundary conditions per block 69C MXCUTS : maximum number of cuts 70C MXPRT : maximum number of print statements per block 71C MXPTS : maximum number of grid points 72C = summation of (JDIM+3)*(KDIM+3)*(IDIM+3) over all blocks 73C MXPPTS : maximum number of points to store PROPERTIES at 74C = 1 for inviscid calculation 75C = MXPTS for viscous calculation 76C MXQNPT : maximum number of points to store additional TIME levels at 77C = 1 for steady calculation 78C = MXPTS for unsteady calculation 79C QN(1:NTIME) stores Q^(n), Q^(n-1), ... 80C QN(NTIME+1:NTMLVL) stores Q^(m), ... 81C MXRKPT : maximum number of points required for Runge-Kutta algorithm 82C = 1 if not using Runge-Kutta or doing turb. model conversion 83C = MXPTS if using Runge-Kutta or doing turb. model conversion 84C MXNODE : maximum number of grid nodes over entire field 85C = summation of (JDIM+2)*(KDIM+2)*(IDIM+2) over all blocks 86C MXBCPT : maximum number of boundary condition points to be 87C stored in the QBC array (for the PROFILE bc type) 88C MXPROF : maximum number of profile boundary condition segments in QBC array 89C MAXQ : maximum number of variables to store 90C = 5 for perfect gas 91C = 7 for perfect gas with 2 equation turbulence model 92C = 12 for perfect gas with Reynolds stress turbulence 93C MAXF : maximum number of fluxes and residuals to store 94C = 5 for perfect gas 95C = 7 for perfect gas with 2 equation turbulence model 96C = 12 for perfect gas with Reynolds stress turbulence 97C MAXP : maximum number of properties 98C = 2 for inviscid or laminar 99C = 3 for turbulent flow 100C = 4 for Reynolds Stress turbulence model 101C = 5 for Algebraic Reynolds Stress turbulence model 102C MAXW : maximum size of additional work array (MXSECT*MAXW) 103C = 19 for source terms: WORKX(1) = TKE 104C WORKX(2) = ETILDE 105C WORKX(3) = XSIW 106C WORKX(4-5) = FDAMP 107C WORKX(6) = EPSC 108C WORKX(7) = PDIL 109C WORKX( 8-13) = SIJ 110C WORKX(14-19) = WIJ 111C MAXTQ : maximum number of wall function data stored 112C = 2 for (tau_wall and qdot_wall) 113C MAXMDL : maximum size of array describing turbulence model (IMODEL) 114C MXTRSG : maximum number of transition specification segments 115C MXSECT : maximum size of a cross section i.e., 116C maximum of (ixj, jxk or ixk) for all blocks 117C used to dimension temporary space for q(+,-) and f. 118C MXSIZE : maximum size of any block, used to dimension RES array 119C = maximum of (JDIM+3)*(KDIM+3)*(IDIM+3) of all blocks 120C MXSIZ4 : maximum size of any block, used to dimension FHAT for 4th order 121C = 1 for 1st or 2nd order schemes 122C = maximum of (JDIM+3)*(KDIM+3)*(IDIM+3) of all blocks for 4th order 123C MXCRS : maximum size of all coarse levels of all blocks used to store 124C the multigrid forcing function 125C MXRANK : maximum rank of Jacobian blocks 126C = 5 for perfect gas 127C = 7 for perfect gas with 2 equation turbulence model 128C = 12 for perfect gas with Reynolds stress turbulence 129C MXABC : maximum number of storage locations for tridiagonal system to invert 130C = NPL*max(IDIM*JDIM,JDIM*KDIM) for (AF3F) 131C where NPL is the number of secondary planes vectorized over 132C MXRM : maximum amount of storage for the RM array 133C = MXSECT for LGS and MARCH 134C = not necessary for RK, AF3F 135C MXSTAG : maximum number of stages for multistage time integration 136C MXTMLV : maximum number of additional time levels for unsteady calculation 137C = 1 for steady (and set MXQNPT = 1) 138C = 2 for iterative, implicit unsteady 139C = 3 for tau-ts iterative unsteady with 2nd order pseudo sub-iter 140C 141 PARAMETER (MAXQ = 12) 142c-- PARAMETER (MXPTS = 50000) 143c++ PARAMETER (MXPTS = 34000) 144c32^3 PARAMETER (MXPTS = 46700) 145c^^ PARAMETER (MXPTS = 314432) 146 PARAMETER (MXPTS = 56240) 147c%% PARAMETER (MXPTS = 146400) 148C 149c++ PARAMETER (MXSIZE = 34000) 150c32^3 PARAMETER (MXSIZE = 46700) 151c^^ PARAMETER (MXSIZE = 314432) 152 PARAMETER (MXSIZE = 41040) 153c%% PARAMETER (MXSIZE = 109200) 154c## PARAMETER (MXSIZ4 = MXSIZE) 155 PARAMETER (MXSIZ4 = 1) 156c^^ PARAMETER (MXCRS = MXPTS/4) 157 PARAMETER (MXCRS = 15200) 158c%% PARAMETER (MXCRS = 37200) 159c-- PARAMETER (MXSECT = 5650) 160c%% PARAMETER (MXSECT = 5000) 161c$$ PARAMETER (MXSECT = 8210) 162c^^ PARAMETER (MXSECT = 26000) 163 PARAMETER (MXSECT = 8208) 164c%% PARAMETER (MXSECT = 21840) 165C 166 PARAMETER (MXTMLV = 3) 167C 168 PARAMETER (MXBLKS = 20) 169 PARAMETER (MXLVLS = 4) 170 PARAMETER (MXCFS = 4) 171 PARAMETER (MXBCS = 10) 172 PARAMETER (MXCUTS = 100) 173 PARAMETER (MXPRT = 20) 174C 175 PARAMETER (MXPPTS = MXPTS) 176 PARAMETER (MXQNPT = MXPTS) 177 PARAMETER (MXNODE = MXPTS) 178 PARAMETER (MXRKPT = MXPTS) 179C 180 PARAMETER (MXBCPT = 1000) 181 PARAMETER (MXPROF = 10) 182C 183 PARAMETER (MAXF = MAXQ) 184 PARAMETER (MAXP = 5) 185 PARAMETER (MAXW = 19) 186 PARAMETER (MAXTQ = 2) 187 PARAMETER (MAXMDL = 13) 188 PARAMETER (MXTRSG = MXBLKS) 189C 190 PARAMETER (MXRANK = MAXQ) 191c PARAMETER (MXABC = 2*MXSECT) 192 PARAMETER (MXABC = MXSECT) 193 PARAMETER (MXRM = MXSECT) 194 PARAMETER (MXSTAG = 6) 195C 196C Commons 197C 198 include '../header/common.h' 199 include '../header/histry.h' 200C 201C PERMANENT STORAGE FOR ALL ZONES/BLOCKS 202C 203C Storage for flow variables, turbulent eddy viscosity, 204C metrics, volumes, jacobians and grid. 205C 206 COMMON /LOCAL/ Q (MXPTS *MAXQ), 207 1 QN (MXQNPT*MAXQ*MXTMLV), 208 2 PROPS (MXPPTS*MAXP), 209 3 S (MXPTS*4*3), 210 4 VOL (MXPTS), 211 5 DTJ (MXPTS), 212 6 R (MXNODE*3), 213 7 QBC (MXBCPT*MAXQ) 214C 215C Dimensions of domains 216C 217 DIMENSION IDIM (MXBLKS,MXLVLS), 218 1 JDIM (MXBLKS,MXLVLS), 219 2 KDIM (MXBLKS,MXLVLS) 220C 221C Extra dimensions for input of grids 222C 223 DIMENSION IDIMIN(MXBLKS), 224 1 JDIMIN(MXBLKS), 225 2 KDIMIN(MXBLKS) 226C 227C Offsets of domains used to calculate location in total storage 228C of Q,S,... of a given domain. 229C IOFFQN offset for previous time level storage for QN 230C IOFFQC offset for coarse grid storage for QCRS and RESCRS. 231C IOFQBC offset for multiple profile boundary condition segments 232C 233 DIMENSION IOFF (MXBLKS,MXLVLS), 234 1 IOFFS (MXBLKS,MXLVLS), 235 2 IOFFQ (MXBLKS,MXLVLS), 236 3 IOFFP (MXBLKS,MXLVLS), 237 4 IOFFQN(MXBLKS,MXLVLS), 238 5 IOFFQC(MXBLKS,MXLVLS), 239 6 IOFQBC(MXPROF,MXLVLS) 240C 241C Boundary condition data. 242C IBCPRF is re-ordered BC data to read in profile BC segments. 243C 244 DIMENSION IBCDAT(10,MXBCS,MXBLKS,MXLVLS), NBCS(MXBLKS), 245 1 IBCPRF(10,MXPROF,MXLVLS) 246C 247C Cut (domain to domain communication) data. 248C 249 DIMENSION ICUTS (21,MXCUTS,MXLVLS) 250 INTEGER NCUTS 251 CHARACTER*20 CUTNAM(MXCUTS) 252C 253C Printout control data. 254C IPRINT, NPRINT control output at end of iteration cycle 255C IPRMON, NPRMON control output during iteration cycle 256C 257 DIMENSION IPRINT(10,MXPRT,MXBLKS), NPRINT(MXBLKS), 258 1 IPRMON(10,MXPRT,MXBLKS), NPRMON(MXBLKS) 259C 260C Data controlling PLOT3D output 261C 262 LOGICAL IFPLT3, BINPLT 263C 264C Data controlling grid type 265C 266 LOGICAL GRDBIN 267C 268C Data controlling the time integration 269C 270 DIMENSION NITS (MXCFS), MGLVLS(MXCFS), ITSLVL(MXLVLS,MXCFS) 271 DIMENSION RKALPH(MXSTAG) 272C 273C Data controlling the accuracy and scheme used. 274C 275 DIMENSION IORDER(3,MXBLKS,MXLVLS,MXCFS), 276 1 LIMIT (3,MXBLKS), 277 2 RKAP (3,MXBLKS), 278 3 IFLUX (MXBLKS), 279 4 ITORDR(3) 280 LOGICAL DIAG (3,MXBLKS) 281C 282C Data controlling the viscous/inviscid calculation 283C 284 LOGICAL VISCOS(3,MXBLKS) 285C 286C Data controlling calculation of delq 287C 288 DIMENSION IFDELQ(MAXQ) 289C 290C Data controlling the turbulence model 291C 292 DIMENSION IMODEL(MAXMDL), 293 1 IMODLX(MAXMDL) 294C 295C Data controlling the specification of transition location 296C 297 DIMENSION ITRSEG(7,MXTRSG) 298C 299C Storage for L2 Norms 300C For Perfect Gas: 301C R2NORM(1) = L2 Norm of the mass conservation equation 302C R2NORM(2-4) = L2 Norm of the momentum equations 303C R2NORM(5) = L2 Norm of the energy equation 304C R2NORM(6) = L2 Norm of all equations 305C R2NORM -> L2 Norm of transport equations (includes dQ/dt for iterative, 306C implicit scheme) 307C R2DQDT -> L2 Norm of dQ/dt terms in iterative, implicit scheme 308C 309 DIMENSION R2NORM(MAXF+1), 310 1 R2DQDT(MAXF+1), 311 2 NBRPTS(MXLVLS) 312C 313C Storage for forces 314C FORCE -> Forces integrated in x, y, z directions 315C RCNTR -> Location to calculate point vortex and moments about 316C 317 DIMENSION FORCE(3), 318 1 RCNTR(3) 319C 320C Data for diagnostic output 321C 322 REAL RESMAX 323 INTEGER IFRSMX, 324 1 IRESMX(0:3) 325 INTEGER NFIXQ (MAXQ), 326 1 NRELIZ(MAXQ), 327 2 NPRLIM, 328 3 NPRNEG 329 DIMENSION QMIN (MAXQ) 330 CHARACTER FIXTYP*7 331C 332C STORAGE REQUIRED FOR RUNGE-KUTTA 333C Required on all zones/blocks for the size of the finest grid 334C N.B. This storage is not needed unless using the RK routine 335C 336 DIMENSION W (MXRKPT*MAXQ) 337C 338C STORAGE REQUIRED FOR MG FORCING FUNCTION AND STORED RESTRICTED Q 339C Required only for coarse grid levels 340C 341 DIMENSION QCRS (MXCRS*MAXQ), 342 1 RESCRS(MXCRS*MAXF) 343C 344C STORAGE REQUIRED FOR A ZONE/BLOCK 345C 346C Storage for residuals 347C For (RK, AF) : store residual over entire block 348C For (LGS, MARCH) : store residuals over a plane of block 349C 350 DIMENSION RES (MXSIZE*MAXF) 351C 352C Storage for fluxes for fourth order 353C 354 DIMENSION FHAT (MXSIZ4*MAXF) 355C 356C Storage for MUSCL and flux evaluation 357C 358 DIMENSION DQP (MXSECT*MAXQ), 359 1 DQM (MXSECT*MAXQ), 360 2 QP (MXSECT*MAXQ), 361 3 QM (MXSECT*MAXQ), 362 4 TAU (MXSECT*6), 363 5 F (MXSECT*MAXF), 364 6 FM (MXSECT*MAXF), 365 7 STEMP (MXSECT*4), 366 8 WORKX (MXSECT*MAXW) 367C 368C Temporary storage for implicit terms 369C For (RK) : not needed 370C For (AF3F) : NPL * max(IDIM*JDIM, JDIM*KDIM) 371C where NPL is the number of secondary 372C planes to vectorize over 373C For (MARCH) : ABC: JDIM*KDIM*2 (store LU of the current plane) 374C RTEMP,AT: JDIM*KDIM 375C 376 DIMENSION ABC (MXABC*MXRANK*MXRANK*3), 377 1 RTEMP (MXABC*MXRANK), 378 2 AT (MXABC*MXRANK*MXRANK*3) 379C 380C Temporary storage for the I contribution to the implicit scheme 381C for LGS and MARCH. 382C 383 DIMENSION RM (MXRM*MXRANK*MXRANK) 384C 385C Identity Matrix required in implicit work 386C 387 DIMENSION AI (MXRANK*MXRANK) 388C 389C Filenames 390C FRDRST File to read restart from 391C FWRRST File to write restart to 392C GRDFIL Array of files to read grid from 393C FPLT3G File to write PLOT3D Grid file to 394C FPLT3Q File to write PLOT3D Q (conserved variables) to 395C FPLT3F File to write PLOT3D function file to 396C FPLT3FN File to write PLOT3D function name file to 397C FPROFL Array of files to read Profile boundary data from 398C FERROR File name output if error 399C 400 CHARACTER*80 FRDRST, FWRRST, GRDFIL(MXBLKS), 401 1 FPLT3G, FPLT3Q, FPLT3F, FPLT3FN, FPROFL(MXPROF), 402 2 FERROR 403C 404C Switches 405C RESTRT -> If true then do a restart 406C VGNRN -> If true then use Vigneron technique when marching 407C CONTRB -> If true then converting turbulence models on restart 408C MRCHNG -> If true then use marching fluxes in residual 409C INITRS -> If true then initialize residuals to zero (used in MG) 410C IFL2DQ -> If true then calculate L2 norm of dQ/dt term (unsteady) 411C YNCALC -> If true then recalculate Y_normal on restart 412C IFMGFF -> If true then calculate the MG forcing function from the 413C coarse grid residual in the MG cycle 414C 415 LOGICAL RESTRT, VGNRN, CONTRB, MRCHNG, INITRS, IFL2DQ, YNCALC, 416 1 IFMGFF 417C 418 419#ifdef CPU_TIME 420C 421C Storage to calculate execution time 422C 423 REAL*4 TARRAY(2), RUNTIME 424#endif 425 426C 427C Version Number for Output 428C 429 CHARACTER*22 VERSN 430C 431C Local storage to output number of grid points 432C 433 DIMENSION NNODES(MXLVLS,2) 434C 435C 0. PRELIMINARIES 436C Output Version Number 437C 438 VERSN = '$Revision: 4.7 $' 439 WRITE (IOUT,1200) VERSN 440 1200 FORMAT (' ','ISAAC - ',A22,/) 441C 442C Set switch to test for storage errors 443C 444 IERROR = 0 445C 446C Set storage for IMODEL to maximum 447C 448 NMDL = MAXMDL 449C 450C Set storage for WORKX to maximum 451C 452 NW = MAXW 453C 454C Set up the Characteristic Array to define the characteristics 455C of the available flux evaluation schemes. 456C 457 IFCHAR(IFROE,1) = IFDS 458C 459C Zero appropriate storage 460C 461 ITTOT = 0 462C 463C Initialize time to zero 464C 465 TIME = 0.E0 466C 467C Initialize force coefficients 468C 469 CL = 0.E0 470 CD = 0.E0 471C 472C 1. READ INPUT 473C Read input data controlling the run 474C 475 NXQ = MAXQ 476 NXCFS = MXCFS 477 NXLVLS = MXLVLS 478 NXBLKS = MXBLKS 479 NXBCS = MXBCS 480 NXCUTS = MXCUTS 481 NXPRT = MXPRT 482 NXSTAG = MXSTAG 483 NXPROF = MXPROF 484 NXSIZE = MXSIZE 485 NXSECT = MXSECT 486 NXPTS = MXPTS 487 NXPPTS = MXPPTS 488 NXQNPT = MXQNPT 489 NXTRSG = MXTRSG 490 CALL RDDATA (NXQ, NXCFS, NXLVLS, NXBLKS, NXBCS, NXCUTS, NXPRT, 491 1 NCFS, MGLVLS, NBLKS, IDIM, JDIM, KDIM, 492 2 NBCS, IBCDAT, NCUTS, ICUTS, CUTNAM, 493 3 NPRINT, IPRINT, NPRMON, IPRMON, NITMON, 494 4 ITUPJ, ITSLVL, NITS, NSUBIT, NITFO, NITBBC, NITALG, 495 5 NITRST, NXSTAG, NSTAGE, RKALPH, 496 6 RESTRT, FRDRST, FWRRST, GRDFIL, IGRDTP, GRDBIN, 497 7 NXPROF, NPROFL, FPROFL, 498 8 IFPLT3, BINPLT, FPLT3G, FPLT3Q, FPLT3F, FPLT3FN, 499 9 IORDER, LIMIT, RKAP, IFLUX, VISCOS, DIAG, IMODEL, 500 A IMSTRT, IMEND, TOLER, VGNRN, SIGMA, IFDELQ, 501 B CONTRB, IMODLX, QMIN, ITRSMX, IFFORC, IFRCPR, 502 C YNCALC, NXTRSG, NTRSEG, ITRSEG, RCNTR) 503C 504 ITURB = IMODEL( 1) 505 IASM = IMODEL(10) 506C 507C 2. VERIFY INPUT 508C 509C Check input data for BC, CUTS, Print statements 510C 511 CALL CHKDAT (NBLKS, IDIM, JDIM, KDIM, VISCOS, 512 1 NXBCS, NBCS, IBCDAT, NCUTS, ICUTS, CUTNAM, 513 2 NXPRT, NPRINT, IPRINT, IERRCD) 514 IF (IERRCD .GT. INO) THEN 515 IERROR = IERROR + 1 516 ENDIF 517C 518C Check to verify that NQ >= NP. This is required due to temporary 519C arrays in DELQ routines. If NP is ever greater than NQ, need to 520C change the temporary storage. 521C 522 IF (NP .GT. NQ) THEN 523 WRITE (IOUT,8000) 524 IERROR = IERROR + 1 525 ENDIF 526 8000 FORMAT (' ','MAIN : ERROR-> NP is greater than NQ.', 527 1 /' ',15X,'This is not allowed due to temporary storage', 528 2 /' ',15X,'allocation in DELQ routines.', 529 3 /' ',15X,'The temporary storage for PROPSI and PROPSC ', 530 4 /' ',15X,'must be modified.') 531C 532C Check the values of constants to make sure that the storage exists 533C 534 IF (NQ .GT. MAXQ ) THEN 535 WRITE (IOUT,1505) NQ, MAXQ 536 IERROR = IERROR + 1 537 ENDIF 538C 539 IF (NF .GT. MAXF ) THEN 540 WRITE (IOUT,1510) NF, MAXF 541 IERROR = IERROR + 1 542 ENDIF 543C 544 IF (NP .GT. MAXP ) THEN 545 WRITE (IOUT,1515) NP, MAXP 546 IERROR = IERROR + 1 547 ENDIF 548C 549 IF (NRANK .GT. MXRANK) THEN 550 WRITE (IOUT,1520) NRANK, MXRANK 551 IERROR = IERROR + 1 552 ENDIF 553C 554 IF (NTMLVL .GT. MXTMLV) THEN 555 WRITE (IOUT,1525) NTMLVL, MXTMLV 556 IERROR = IERROR + 1 557 ENDIF 558C 559 1505 FORMAT (/' ','MAIN : Number of variables required ',I5, 560 1 ' exceeds allocated ',I5, 561 2 /' ', 8X,'=> Increase MAXQ and re-submit.') 562 1510 FORMAT (/' ','MAIN : Number of fluxes required ',I5, 563 1 ' exceeds allocated ',I5, 564 2 /' ', 8X,'=> Increase MAXF and re-submit.') 565 1515 FORMAT (/' ','MAIN : Number of properties required ',I5, 566 1 ' exceeds allocated ',I5, 567 2 /' ', 8X,'=> Increase MAXP and re-submit.') 568 1520 FORMAT (/' ','MAIN : Size of matrices required ',I5, 569 1 ' exceeds allocated ',I5, 570 2 /' ', 8X,'=> Increase MXRANK and re-submit.') 571 1525 FORMAT (/' ','MAIN : Number of time levels required',I5, 572 1 ' exceeds allocated ',I5, 573 2 /' ', 8X,'=> Increase MXTMLV and re-submit.') 574C 575C Set up the identity matrix used in implicit terms 576C 577 CALL IDENT (AI) 578C 579C 3. CALCULATE POINTERS 580C a. Calculate dimensions for coarser grids and offset pointers. 581C b. Calculate offset pointers and boundary data. 582C c. Calculate ICUTS for coarser grids. 583C d. Verify coarse grid BC and CUT data 584C e. Set up IORDER for multigrid 585C f. Output grid summary 586C 587 NDSIZE = 0 588 NDCRS = 0 589 NDSECT = 0 590 NDRM = 0 591 NLVLS = 1 592 DO 10 ICFS = 1, NCFS 593 ILVL = ICFS + MGLVLS(ICFS) - 1 594 IF (ILVL .GT. NLVLS) NLVLS = ILVL 595 10 CONTINUE 596 IOFFT = 0 597 IOFFST = 0 598 IOFFQT = 0 599 IOFFPT = 0 600 IOFQNT = 0 601 IOFQCT = 0 602C 603C 3-a. Calculate dimensions for coarser grids and test allowable grid levels. 604C 605 MSTLVL = NLVLS 606 DO 20 ILVL = 2, NLVLS 607 DO 15 IBLK = 1, NBLKS 608 IDIM(IBLK,ILVL) = IDIM(IBLK,ILVL-1) / 2 + 1 609 JDIM(IBLK,ILVL) = JDIM(IBLK,ILVL-1) / 2 + 1 610 KDIM(IBLK,ILVL) = KDIM(IBLK,ILVL-1) / 2 + 1 611C 612 IDIMF = 2 * (IDIM(IBLK,ILVL) - 1) + 1 613 JDIMF = 2 * (JDIM(IBLK,ILVL) - 1) + 1 614 IF (THREED) THEN 615 KDIMF = 2 * (KDIM(IBLK,ILVL) - 1) + 1 616 ELSE 617 KDIMF = 2 618 ENDIF 619C 620 IF (IDIMF .NE. IDIM(IBLK,ILVL-1) .OR. 621 1 JDIMF .NE. JDIM(IBLK,ILVL-1) .OR. 622 2 KDIMF .NE. KDIM(IBLK,ILVL-1) ) THEN 623 MSTLVL = MIN (MSTLVL, (ILVL-1)) 624 ENDIF 625 15 CONTINUE 626 20 CONTINUE 627C 628 IF (NLVLS .GT. MSTLVL) THEN 629 WRITE (IOUT,1530) NLVLS, MSTLVL 630 NLVLS = MSTLVL 631 ENDIF 632 IF (NCFS .GT. NLVLS ) THEN 633 WRITE (IOUT,1535) NCFS, NLVLS 634 NCFS = NLVLS 635 ENDIF 636 1530 FORMAT (/' ','MAIN : Number of grid levels requested ',I3, 637 1 ' exceeds allowable ',I3, 638 2 /' ', 8X,'=> Coarser levels are being deleted and ', 639 3 'run continuing.') 640 1535 FORMAT (/' ','MAIN : Number of coarsenings requested ',I3, 641 1 ' exceeds allowable ',I3, 642 2 /' ', 8X,'=> Coarsenings are being deleted and ', 643 3 'run continuing.'//) 644C 645C Verify all coarsenings fall within NLVLS number of grid levels 646C 647 DO 25 ICFS = 1, NCFS 648 LVLTST = ICFS + MGLVLS(ICFS) - 1 649 IF (LVLTST .GT. NLVLS) THEN 650 MGLVLS(ICFS) = NLVLS - ICFS + 1 651 ENDIF 652 25 CONTINUE 653C 654C 3-b. Calculate offset pointers and boundary data. 655C 656 DO 50 ILVL = 1, NLVLS 657 NBRPTS(ILVL) = 0 658 DO 40 IBLK = 1, NBLKS 659 IF (ILVL .NE. 1) THEN 660C 661C Calculate storage required on coarser grids for forcing function 662C 663 NDCRS = NDCRS + (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 664 1 (KDIM(IBLK,ILVL)+3) 665C 666C Calculate IBCDAT for coarser grids. 667C 668 DO 35 IBC = 1, NBCS(IBLK) 669 IBCDAT( 1,IBC,IBLK,ILVL)=IBCDAT(1,IBC,IBLK,ILVL-1) 670 IBCDAT( 2,IBC,IBLK,ILVL)=IBCDAT(2,IBC,IBLK,ILVL-1) 671 IBCDAT( 3,IBC,IBLK,ILVL)=IBCDAT(3,IBC,IBLK,ILVL-1)/2+1 672 IBCDAT( 4,IBC,IBLK,ILVL)=IBCDAT(4,IBC,IBLK,ILVL-1)/2+1 673 IBCDAT( 5,IBC,IBLK,ILVL)=IBCDAT(5,IBC,IBLK,ILVL-1)/2+1 674 IBCDAT( 6,IBC,IBLK,ILVL)=IBCDAT(6,IBC,IBLK,ILVL-1)/2+1 675 IBCDAT( 7,IBC,IBLK,ILVL)=IBCDAT(7,IBC,IBLK,ILVL-1)/2+1 676 IBCDAT( 8,IBC,IBLK,ILVL)=IBCDAT(8,IBC,IBLK,ILVL-1)/2+1 677 IBCDAT( 9,IBC,IBLK,ILVL)=IBCDAT(9,IBC,IBLK,ILVL-1)/2+1 678 IBCDAT(10,IBC,IBLK,ILVL)=IBCDAT(10,IBC,IBLK,ILVL-1) 679 35 CONTINUE 680 ENDIF 681 NBRPTS(ILVL) = NBRPTS(ILVL) + (IDIM(IBLK,ILVL)-1)* 682 1 (JDIM(IBLK,ILVL)-1)*(KDIM(IBLK,ILVL)-1) 683 IOFF (IBLK,ILVL) = IOFFT 684 IOFFS (IBLK,ILVL) = IOFFST 685 IOFFQ (IBLK,ILVL) = IOFFQT 686 IOFFP (IBLK,ILVL) = IOFFPT 687 IOFFQN(IBLK,ILVL) = IOFQNT 688 IOFFQC(IBLK,ILVL) = IOFQCT 689 IOFFT = IOFFT + IDIM(IBLK,ILVL)*JDIM(IBLK,ILVL)* 690 1 KDIM(IBLK,ILVL) 691 IOFFST = IOFFST + (IDIM(IBLK,ILVL)+2)*(JDIM(IBLK,ILVL)+2)* 692 1 (KDIM(IBLK,ILVL)+2) 693 IOFFQT = IOFFQT + (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 694 1 (KDIM(IBLK,ILVL)+3) 695 IF (ILVL .NE. 1) THEN 696 IOFQCT = IOFQCT+(IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 697 1 (KDIM(IBLK,ILVL)+3) 698 ENDIF 699C 700 NDSIZE = MAX(NDSIZE,(IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 701 1 (KDIM(IBLK,ILVL)+3)) 702 NDSECT = MAX(NDSECT,(IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)) 703 NDSECT = MAX(NDSECT,(JDIM(IBLK,ILVL)+3)*(KDIM(IBLK,ILVL)+3)) 704C 705 IF (ISOLVR .EQ. ILGS .OR. ISOLVR .EQ. IMARCH) THEN 706 NDRM = MAX(NDRM,(JDIM(IBLK,ILVL)-1)*(KDIM(IBLK,ILVL)-1)) 707 ELSE 708 NDRM = 0 709 ENDIF 710C 711C If running inviscid case (ITURB < ITLMNR), then don't 712C increment the pointer to the PROPS array. This allows to not 713C store this data for the inviscid case to cut down on memory. 714C 715 INCP = (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 716 1 (KDIM(IBLK,ILVL)+3) 717 IF (ITURB .LT. ITLMNR) THEN 718 INCP = 0 719 ENDIF 720 IOFFPT = IOFFPT + INCP 721C 722C If running steady case, then don't increment pointer to previous 723C time level data to allow minimum memory configuration for steady cases. 724C 725 INCQN = (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 726 1 (KDIM(IBLK,ILVL)+3) 727 IF (ITIMED .EQ. ISTDY) THEN 728 INCQN = 0 729 ENDIF 730 IOFQNT = IOFQNT + INCQN 731 40 CONTINUE 732C 733C 3-c. Calculate ICUTS for coarser grids. 734C 735 IF (ILVL .NE. 1) THEN 736 DO 45 ICUT = 1, NCUTS 737 ICUTS( 1,ICUT,ILVL) = ICUTS( 1,ICUT,ILVL-1) 738 ICUTS( 2,ICUT,ILVL) = ICUTS( 2,ICUT,ILVL-1) 739 ICUTS( 3,ICUT,ILVL) = ICUTS( 3,ICUT,ILVL-1) / 2 + 1 740 ICUTS( 4,ICUT,ILVL) = ICUTS( 4,ICUT,ILVL-1) / 2 + 1 741 ICUTS( 5,ICUT,ILVL) = ICUTS( 5,ICUT,ILVL-1) / 2 + 1 742 ICUTS( 6,ICUT,ILVL) = ICUTS( 6,ICUT,ILVL-1) / 2 + 1 743 ICUTS( 7,ICUT,ILVL) = ICUTS( 7,ICUT,ILVL-1) / 2 + 1 744 ICUTS( 8,ICUT,ILVL) = ICUTS( 8,ICUT,ILVL-1) / 2 + 1 745 ICUTS( 9,ICUT,ILVL) = ICUTS( 9,ICUT,ILVL-1) / 2 + 1 746 ICUTS(10,ICUT,ILVL) = ICUTS(10,ICUT,ILVL-1) 747 ICUTS(11,ICUT,ILVL) = ICUTS(11,ICUT,ILVL-1) 748 ICUTS(12,ICUT,ILVL) = ICUTS(12,ICUT,ILVL-1) / 2 + 1 749 ICUTS(13,ICUT,ILVL) = ICUTS(13,ICUT,ILVL-1) / 2 + 1 750 ICUTS(14,ICUT,ILVL) = ICUTS(14,ICUT,ILVL-1) / 2 + 1 751 ICUTS(15,ICUT,ILVL) = ICUTS(15,ICUT,ILVL-1) / 2 + 1 752 ICUTS(16,ICUT,ILVL) = ICUTS(16,ICUT,ILVL-1) / 2 + 1 753 ICUTS(17,ICUT,ILVL) = ICUTS(17,ICUT,ILVL-1) / 2 + 1 754 ICUTS(18,ICUT,ILVL) = ICUTS(18,ICUT,ILVL-1) / 2 + 1 755 ICUTS(19,ICUT,ILVL) = ICUTS(19,ICUT,ILVL-1) 756 ICUTS(20,ICUT,ILVL) = ICUTS(20,ICUT,ILVL-1) 757 ICUTS(21,ICUT,ILVL) = ICUTS(21,ICUT,ILVL-1) 758 45 CONTINUE 759 ENDIF 760 50 CONTINUE 761C 762C 3-d. Verify coarse grid BC and CUT data 763C 764 IF (NLVLS .GT. 1) THEN 765 CALL CHKCRS (NLVLS, NXBLKS, NBLKS, IDIM, JDIM, KDIM, 766 1 NXBCS, NBCS, IBCDAT, 767 2 NXCUTS, NCUTS, ICUTS, CUTNAM, IERRCD) 768 IF (IERRCD .GT. INO) THEN 769 IERROR = IERROR + 1 770 ENDIF 771 ENDIF 772C 773C 3-e. Set up IORDER for multigrid 774C 775 IHIGH = 2 776 ILOW = 1 777 IF (FOURTH) THEN 778 IHIGH = 4 779 ILOW = 4 780 ENDIF 781C 782 DO IBLK = 1, NBLKS 783 DO ICFS = 1, NCFS 784 LVLONE = NCFS - ICFS + 1 785 DO ILVL = 1, NLVLS 786 IF (ILVL .EQ. LVLONE) THEN 787 IORDER(1,IBLK,ILVL,ICFS) = IHIGH 788 IORDER(2,IBLK,ILVL,ICFS) = IHIGH 789 IORDER(3,IBLK,ILVL,ICFS) = IHIGH 790 ELSE 791 IORDER(1,IBLK,ILVL,ICFS) = ILOW 792 IORDER(2,IBLK,ILVL,ICFS) = ILOW 793 IORDER(3,IBLK,ILVL,ICFS) = ILOW 794 ENDIF 795 END DO 796 END DO 797 END DO 798C 799C 3-f. Output grid summary 800C 801 WRITE (IOUT,1540) (ILVL, ILVL = 1, NLVLS) 802 DO ILVL = 1, NLVLS 803 NNODES(ILVL,2) = 0 804 END DO 805C 806 DO IBLK = 1, NBLKS 807 DO ILVL = 1, NLVLS 808 IF (THREED) THEN 809 K = KDIM(IBLK,ILVL) 810 ELSE 811 K = 1 812 ENDIF 813 NNODES(ILVL,1) = IDIM(IBLK,ILVL) * JDIM(IBLK,ILVL) * K 814 NNODES(ILVL,2) = NNODES(ILVL,2) + NNODES(ILVL,1) 815 END DO 816 WRITE (IOUT,1541) IBLK,IDIM(IBLK,1),JDIM(IBLK,1),KDIM(IBLK,1), 817 1 (NNODES(ILVL,1), ILVL = 1, NLVLS) 818 END DO 819 WRITE (IOUT,1542) (NNODES(ILVL,2), ILVL = 1, NLVLS) 820 WRITE (IOUT,1543) 821C 822 1540 FORMAT (///' ','Grid Summary', 823 1 /' ',27X,'Number of Grid Points by Mesh Level', 824 2 /' ',11X,'IDIM',1X,'JDIM',1X,'KDIM',6(5x,I2,2X)) 825 1541 FORMAT ( ' ','Block ',I3,':',3I5,6(1X,I8)) 826 1542 FORMAT ( /' ','Total ',19X,6(1X,I8)) 827 1543 FORMAT (///' ') 828C 829C 4. VERIFY STORAGE 830C Verify sufficient storage available or STOP 831C 832 IF (IOFFQT .GT. MXPTS ) THEN 833 WRITE (IOUT,1000) IOFFQT, MXPTS 834 IERROR = IERROR + 1 835 ENDIF 836C 837 IF (IOFFPT .GT. MXPPTS) THEN 838 WRITE (IOUT,1001) IOFFPT, MXPPTS 839 IERROR = IERROR + 1 840 ENDIF 841C 842 IF (IOFQNT .GT. MXQNPT) THEN 843 WRITE (IOUT,1002) IOFQNT, MXQNPT 844 IERROR = IERROR + 1 845 ENDIF 846C 847 IF (NDSIZE .GT. MXSIZE) THEN 848 WRITE (IOUT,1003) NDSIZE, MXSIZE 849 IERROR = IERROR + 1 850 ENDIF 851C 852 IF (NDCRS .GT. MXCRS ) THEN 853 WRITE (IOUT,1004) NDCRS , MXCRS 854 IERROR = IERROR + 1 855 ENDIF 856C 857 IF (NDSECT .GT. MXSECT) THEN 858 WRITE (IOUT,1005) NDSECT, MXSECT 859 IERROR = IERROR + 1 860 ENDIF 861C 862 IF (NDRM .GT. MXRM ) THEN 863 WRITE (IOUT,1006) NDRM, MXRM 864 IERROR = IERROR + 1 865 ENDIF 866C 867 IF (FOURTH) THEN 868 IF (NDSIZE .GT. MXSIZ4) THEN 869 WRITE (IOUT,1007) NDSIZE, MXSIZ4 870 IERROR = IERROR + 1 871 ENDIF 872 ENDIF 873C 874 IF (IERROR .GT. 0 ) THEN 875 WRITE (IOUT,1009) IERROR 876 STOP 877 ENDIF 878C 879 1000 FORMAT (/' ','MAIN : Storage required for Q variables ', 880 1 I9,' exceeds allocated ',I9, 881 2 /' ', 8X,'=> Increase MXPTS in MAIN and re-submit.') 882 1001 FORMAT (/' ','MAIN : Storage required for PROPS ', 883 1 I9,' exceeds allocated ',I9, 884 2 /' ', 8X,'=> Increase MXPPTS in MAIN and re-submit.') 885 1002 FORMAT (/' ','MAIN : Storage required for QN ', 886 1 I9,' exceeds allocated ',I9, 887 2 /' ', 8X,'=> Increase MXQNPT in MAIN and re-submit.') 888 1003 FORMAT (/' ','MAIN : Storage required for RES ', 889 1 I9,' exceeds allocated ',I9, 890 2 /' ', 8X,'=> Increase MXSIZE in MAIN and re-submit.') 891 1004 FORMAT (/' ','MAIN : Storage required for RESCRS ', 892 1 I9,' exceeds allocated ',I9, 893 2 /' ', 8X,'=> Increase MXCRS in MAIN and re-submit.') 894 1005 FORMAT (/' ','MAIN : Storage required for cross-section', 895 1 I9,' exceeds allocated ',I9, 896 2 /' ', 8X,'=> Increase MXSECT in MAIN and re-submit.') 897 1006 FORMAT (/' ','MAIN : Storage required for RM array ', 898 1 I9,' exceeds allocated ',I9, 899 2 /' ', 8X,'=> Increase MXRM in MAIN and re-submit.') 900 1007 FORMAT (/' ','MAIN : Storage required for FHAT 4th ordr', 901 1 I9,' exceeds allocated ',I9, 902 2 /' ', 8X,'=> Increase MXSIZ4 in MAIN and re-submit.') 903 1009 FORMAT (/' ','MAIN : Storage or input errors have occured = ',I3, 904 1 /' ', 8X,'Please fix and re-submit.', 905 2 /' ', 8X,'RUN ABORTING!') 906C 907C 5. INPUT GRID 908C a. Read in grid for all blocks. 909C b. Loop over cuts and check the validity of the grid cuts 910C c. Extend each block grid from 1:DIM,... to 0:DIM+1,... and calculate 911C coarser grids as every other grid point from next finer grid. 912C d. Loop over cuts on all levels blocks to set image points of the grid. 913C e. Calculate metrics and volumes, reset volumes to be non-zero for 914C boundary conditions that do not lie on a face 915C 916C 5-a. Read in grid for all blocks, finest grid level. 917C 918 CALL RDGRID (NBLKS, IDIM, JDIM, KDIM, IOFFS, GRDFIL, R, 919 1 IDIMIN, JDIMIN, KDIMIN, IGRDTP, GRDBIN, IEOF, IERRCD) 920 IF (IEOF .NE. INO .OR. IERRCD .NE. INO) THEN 921 STOP 922 ENDIF 923C 924C 5-b. Loop over cuts and check the validity of the grid cuts 925C 926 IERRCD = INO 927 ILVL = 1 928 DO 55 ICUT = 1, NCUTS 929 IBLK1 = ICUTS( 1,ICUT,ILVL) 930 IBLK2 = ICUTS(10,ICUT,ILVL) 931 IR1 = IOFFS(IBLK1,ILVL) * 3 + 1 932 IR2 = IOFFS(IBLK2,ILVL) * 3 + 1 933 CALL CHKCUT (IBLK1, IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 934 1 KDIM(IBLK1,ILVL), R(IR1), 935 2 IBLK2, IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 936 3 KDIM(IBLK2,ILVL), R(IR2), ICUTS(1,ICUT,ILVL), 937 4 CUTNAM(ICUT), IERRCD) 938 55 CONTINUE 939 IF (IERRCD .NE. INO) THEN 940 STOP 941 ENDIF 942C 943C 5-c. Extend each block grid and calculate coarser meshes. 944C 945 DO 65 IBLK = 1, NBLKS 946 ILVL = 1 947 IR = IOFFS(IBLK,ILVL) * 3 + 1 948 CALL GDXTND (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 949 1 KDIM(IBLK,ILVL), R(IR)) 950C 951 DO 60 ILVL = 2, NLVLS 952 IRF = IOFFS(IBLK,ILVL-1) * 3 + 1 953 IR = IOFFS(IBLK,ILVL ) * 3 + 1 954 CALL CRSGRD (IDIM(IBLK,ILVL-1), JDIM(IBLK,ILVL-1), 955 1 KDIM(IBLK,ILVL-1), R(IRF), 956 2 IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 957 3 KDIM(IBLK,ILVL), R(IR)) 958 CALL GDXTND (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 959 1 KDIM(IBLK,ILVL), R(IR)) 960 60 CONTINUE 961 65 CONTINUE 962C 963C 5-d. Loop over cuts setting grid image points 964C 965 DO 75 ILVL = 1, NLVLS 966 DO 70 ICUT = 1, NCUTS 967 IBLK1 = ICUTS( 1,ICUT,ILVL) 968 IBLK2 = ICUTS(10,ICUT,ILVL) 969 IR1 = IOFFS(IBLK1,ILVL) * 3 + 1 970 IR2 = IOFFS(IBLK2,ILVL) * 3 + 1 971 CALL GRDCUT (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 972 1 KDIM(IBLK1,ILVL), R(IR1), 973 2 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 974 3 KDIM(IBLK2,ILVL), R(IR2), ICUTS(1,ICUT,ILVL)) 975 70 CONTINUE 976 75 CONTINUE 977C 978C 5-e. Calculate metrics and volumes and reset volumes for some bcs. 979C Reset interior volumes on coarser meshes to be the collection of the 980C fine grid volumes to ensure that the coarse grid transfer operator 981C generates a smooth and accurate coarse grid representation of the fine 982C grid solution. 983C 984 IERRCD = INO 985 DO 85 IBLK = 1, NBLKS 986 DO 80 ILVL = 1, NLVLS 987 IR = IOFFS(IBLK,ILVL) * 3 + 1 988 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 989 IV = IOFFS(IBLK,ILVL) + 1 990 CALL METRIC (IBLK, ILVL, IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 991 1 KDIM(IBLK,ILVL), R(IR), S(IS), VOL(IV), IERRCD) 992 CALL METBC (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 993 1 KDIM(IBLK,ILVL), S(IS), VOL(IV), 994 2 NBCS(IBLK), IBCDAT(1,1,IBLK,ILVL)) 995 IF (ILVL .GT. 1) THEN 996 IVF = IOFFS(IBLK,ILVL-1) + 1 997 CALL COLV2H (IDIM(IBLK,ILVL-1), JDIM(IBLK,ILVL-1), 998 1 KDIM(IBLK,ILVL-1), VOL(IVF), 999 2 IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1000 3 KDIM(IBLK,ILVL), VOL(IV), QP) 1001 ENDIF 1002 80 CONTINUE 1003 85 CONTINUE 1004 IF (IERRCD .NE. INO) THEN 1005 STOP 1006 ENDIF 1007C 1008C 6. INITIALIZATION 1009C Initialize each block and each grid level to freestream. 1010C Initialize each block and each grid level properties. 1011C Initialize each block and each grid level at additional time levels. 1012C Initialize each block and each grid level wall functions. 1013C 1014 DO 100 ILVL = 1, NLVLS 1015 DO 90 IBLK = 1, NBLKS 1016 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1017 IV = IOFFS(IBLK,ILVL) + 1 1018 CALL INIT (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1019 1 KDIM(IBLK,ILVL), Q(IQ), ITURB) 1020C 1021 IF (IASM .EQ. IASMGS) THEN 1022 IP = IOFFP(IBLK,ILVL) * NP + 1 1023 CALL INITPR (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1024 1 KDIM(IBLK,ILVL), PROPS(IP), ITURB) 1025 ENDIF 1026C 1027 IF (ITIMED .EQ. IUNSTD) THEN 1028 IQN = IOFFQN(IBLK,ILVL) * NQ * NTMLVL + 1 1029 CALL INITQN (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1030 1 KDIM(IBLK,ILVL), Q(IQ), QN(IQN), 1031 2 QP, QM, IMODEL) 1032 ENDIF 1033C 1034 IF (IFWALF) THEN 1035 IP = IOFFP(IBLK,ILVL) * NP + 1 1036 CALL INITTQ (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1037 1 KDIM(IBLK,ILVL), PROPS(IP), 1038 2 NBCS(IBLK), IBCDAT(1,1,IBLK,ILVL)) 1039 ENDIF 1040 90 CONTINUE 1041 100 CONTINUE 1042C 1043C Loop over all boundary conditions to check for a PROFILE or 1044C SUBSONIC PROFILE bc, then read in the boundary data for the 1045C PROFILE or SUBSONIC PROFILE bcs on the fine grid 1046C and calculate the coarsened profile data for coarse grids. 1047C Set all offsets to zero in case there are no profile boundary segments. 1048C i. Re-order profile BC segment data based on profile number 1049C ii. Loop over all profile segments in order and calculate offsets, 1050C read in profile data and coarsen data 1051C 1052C i. Re-order profile BC segment data based on profile number 1053C 1054 DO 115 IBLK = 1, NBLKS 1055 DO 110 IBC = 1, NBCS(IBLK) 1056 IBCTYP = IBCDAT(1,IBC,IBLK,1) 1057 IF (IBCTYP .EQ. IPROFL .OR. IBCTYP .EQ. IPROSB) THEN 1058 INPROF = IBCDAT(10,IBC,IBLK,1) 1059 DO 105 ILVL = 1, NLVLS 1060 DO 104 L = 1, 10 1061 IBCPRF(L,INPROF,ILVL) = IBCDAT(L,IBC,IBLK,ILVL) 1062 104 CONTINUE 1063 105 CONTINUE 1064 ENDIF 1065 110 CONTINUE 1066 115 CONTINUE 1067C 1068C ii. Loop over all profile segments in order and calculate offsets, 1069C read in profile data and coarsen data 1070C 1071 NDBCPT = 0 1072 IERRCD = INO 1073 IERRC2 = INO 1074 DO 125 ILVL = 1, NLVLS 1075 IOFQBC(1,ILVL) = 0 1076 DO 120 INPROF = 1, NPROFL 1077 ISTRT = IBCPRF( 4,INPROF,ILVL) 1078 IEND = IBCPRF( 5,INPROF,ILVL) 1079 JSTRT = IBCPRF( 6,INPROF,ILVL) 1080 JEND = IBCPRF( 7,INPROF,ILVL) 1081 KSTRT = IBCPRF( 8,INPROF,ILVL) 1082 KEND = IBCPRF( 9,INPROF,ILVL) 1083C 1084 IOFQBC(INPROF,ILVL) = NDBCPT 1085 IQBC = IOFQBC(INPROF,ILVL) * NQ + 1 1086C 1087 NPTS = (IEND-ISTRT+1) * (JEND-JSTRT+1) * 1088 1 (KEND-KSTRT+1) 1089 NDBCPT = NDBCPT + NPTS 1090C 1091 IF (NDBCPT .GT. MXBCPT) THEN 1092 IERRCD = IYES 1093 ENDIF 1094C 1095 IF (IERRCD .EQ. INO) THEN 1096 IF (ILVL .EQ. 1) THEN 1097 WRITE (IOUT,1410) FPROFL(INPROF) 1098 FERROR = FPROFL(INPROF) 1099 OPEN (IRDBC,FILE=FPROFL(INPROF),STATUS='OLD', 1100 1 FORM='FORMATTED',ERR=990) 1101 CALL RDPRFL (NPTS, QBC(IQBC), 1102 2 IBCPRF(1,INPROF,ILVL), IERRC2) 1103 ELSE 1104 ISTRTF = IBCPRF( 4,INPROF,ILVL-1) 1105 IENDF = IBCPRF( 5,INPROF,ILVL-1) 1106 JSTRTF = IBCPRF( 6,INPROF,ILVL-1) 1107 JENDF = IBCPRF( 7,INPROF,ILVL-1) 1108 KSTRTF = IBCPRF( 8,INPROF,ILVL-1) 1109 KENDF = IBCPRF( 9,INPROF,ILVL-1) 1110 NPTSF = (IENDF-ISTRTF+1) * (JENDF-JSTRTF+1) * 1111 1 (KENDF-KSTRTF+1) 1112 IQBCF = IOFQBC(INPROF,ILVL-1) * NQ + 1 1113 CALL CRSPRF (NPTSF, QBC(IQBCF), 1114 1 IBCPRF(1,INPROF,ILVL-1), 1115 2 NPTS, QBC(IQBC), 1116 3 IBCPRF(1,INPROF,ILVL), IERRCD) 1117 ENDIF 1118 ENDIF 1119 120 CONTINUE 1120 125 CONTINUE 1121C 1122 IF (IERRCD .NE. INO) THEN 1123 WRITE (IOUT,1420) NDBCPT, MXBCPT 1124 ENDIF 1125C 1126 IF (IERRCD .NE. INO .OR. IERRC2 .NE. INO) THEN 1127 STOP 1128 ENDIF 1129C 1130 1410 FORMAT (//' ','MAIN : Reading input profile from the file:', 1131 1 /' ',10X,A80) 1132 1420 FORMAT (' ','MAIN : ERROR-> Number of boundary points ',I10, 1133 1 /' ',' exceeds allocated (MXBCPT) ',I10, 1134 2 /' ',' Increase MXBCPT in MAIN and resubmit.', 1135 3 /' ', 8X,'RUN ABORTING!') 1136C 1137C 7. RESTART 1138C If this is a restart solution, then read in restart data. 1139C The logic here is a little confusing so pay attention: 1140C (The restart/conversion process plays with NQ, NP. It must reset NQ, NP 1141C to be the correct value before leaving this section) 1142C 1143C Begin{Restart/Conversion} 1144C Open Files 1145C Set NQTO,NFTO, NPTO <- the number of variables of TO 1146C turbulence model 1147C (this is NQ,NP as set in RDDATA) 1148C Set NQFROM,NFFROM,NPFROM <- the number of variables of FROM 1149C turbulence model 1150C Set NQ <- NQFROM, NF <- NFFROM, NP <- NPFROM 1151C For Each Block 1152C If NORMAL Restart ( .not. contrb ) 1153C begin 1154C 1. Read restart data into Q, PROPS 1155C 2. If unsteady, convert QN to conserved variables 1156C N.B. storage for tke in convert is over entire block! 1157C end 1158C Else (doing a restart with turbulence model conversion) 1159C begin 1160C 1. Read restart data of NQFROM variables into RES, 1161C NPFROM properties into W 1162C (RES is dimensioned to hold 1 block of Q, 1163C W is large enough to hold 1 block of PROPS) 1164C 2. Call conversion routine to convert from one turbulence 1165C model to another. Returns the new model in Q, PROPS. 1166C Note: The PROPS will be updated after this section. 1167C 3. If unsteady, output error as this is not yet supported 1168C end 1169C Endif 1170C Next Block 1171C Set NQ <- NQTO, NF <- NFTO, NP <- NPTO 1172C Close Files 1173C End{Restart/Conversion} 1174C 1175 IF ( RESTRT ) THEN 1176 WRITE (IOUT,1500) 1177 1500 FORMAT (' ','MAIN : Reading Restart File'/) 1178 FERROR = FRDRST 1179 OPEN (IRDRST,FILE=FRDRST,STATUS='OLD',FORM='UNFORMATTED', 1180 1 ERR=990) 1181 REWIND (IRDRST) 1182C 1183 NQTO = NQ 1184 NFTO = NF 1185 NPTO = NP 1186C 1187 IF ( .NOT. CONTRB ) THEN 1188 NQFROM = NQ 1189 NFFROM = NF 1190 NPFROM = NP 1191 ELSE 1192 ITURBX = IMODLX( 1) 1193 IF (ITURBX .EQ. ITLMNR .OR. ITURBX .EQ. ITBLMX) THEN 1194 NQFROM = 5 1195 NFFROM = 5 1196 NPFROM = 5 1197 ELSE IF (ITURBX .EQ. ITKE .OR. ITURBX .EQ. ITKW) THEN 1198 NQFROM = 7 1199 NFFROM = 7 1200 NPFROM = 5 1201 ELSE IF (ITURBX .EQ. ITRS) THEN 1202 NQFROM = 12 1203 NFFROM = 12 1204 NPFROM = 5 1205 ELSE 1206 WRITE (IOUT,1600) ITURBX 1207 1600 FORMAT (' ','MAIN : ERROR-> Invalid FROM turbulence model = ',I5, 1208 1 /' ', 8X,'RUN ABORTING!') 1209 STOP 1210 ENDIF 1211 ENDIF 1212C 1213C Now reset NQ, NF, NP to be NQFROM, NFFROM, NPFROM 1214C 1215 NQ = NQFROM 1216 NF = NFFROM 1217 NP = NPFROM 1218C 1219C Loop over all grid levels that are written in the restart file 1220C 1221 ILVL = 1 1222 130 CONTINUE 1223 DO 135 IBLK = 1, NBLKS 1224 IQ = IOFFQ (IBLK,ILVL) * NQTO + 1 1225 IQN = IOFFQN(IBLK,ILVL) * NQTO * NTMLVL + 1 1226 IP = IOFFP (IBLK,ILVL) * NPTO + 1 1227 IS = IOFFS (IBLK,ILVL) * 4 * 3 + 1 1228 IV = IOFFS (IBLK,ILVL) + 1 1229 IF ( .NOT. CONTRB ) THEN 1230 CALL RDREST (ITURB , IBLK, ILVL, IDIM(IBLK,ILVL), 1231 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1232 2 Q(IQ), QN(IQN), PROPS(IP), IEOF, IERRCD) 1233C 1234 IF (IEOF .EQ. IYES) THEN 1235 GO TO 139 1236 ENDIF 1237 IF (IERRCD .NE. INO ) THEN 1238 WRITE (IOUT,1545) 1239 1545 FORMAT (' ','MAIN : ERROR-> Error in restart grid level ',I3, 1240 1 /' ', 8X,'RUN ABORTING!') 1241 STOP 1242 ENDIF 1243C 1244 IF (ITIMED .EQ. IUNSTD) THEN 1245 IF (INITTD) THEN 1246 IQ = IOFFQ (IBLK,ILVL) * NQ + 1 1247 IQN = IOFFQN(IBLK,ILVL) * NQ * NTMLVL + 1 1248 CALL INITQN (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1249 1 KDIM(IBLK,ILVL), Q(IQ), QN(IQN), 1250 2 QP, QM, IMODEL) 1251 ELSE 1252 NNPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1253 1 (KDIM(IBLK,ILVL)+3) * (NTIME-1) 1254 CALL CONP2C (ITURB, NNPTS, QN(IQN), RES) 1255 ENDIF 1256 ENDIF 1257C 1258 ELSE 1259 ITURBX = IMODLX( 1) 1260 CALL RDREST (ITURBX, IBLK, ILVL, IDIM(IBLK,ILVL), 1261 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1262 2 RES, QN(IQN), W, IEOF, IERRCD) 1263C 1264 IF (IEOF .EQ. IYES) THEN 1265 GO TO 139 1266 ENDIF 1267 IF (IERRCD .NE. INO ) THEN 1268 STOP 1269 ENDIF 1270 IF (ITIMED .EQ. IUNSTD) THEN 1271 WRITE (IOUT,1555) 1272 1555 FORMAT (' ','MAIN : ERROR-> Conversion of unsteady calculation ', 1273 1 'is not supported.', 1274 2 /' ', 8X,'RUN ABORTING!') 1275 STOP 1276 ENDIF 1277C 1278C Temporary variables passed to the conversion routine 1279C QP for QC (variables in one-dimensional array) 1280C FM for PROPSC (properties in one-dimensional array - this is 1281C passed in using FM. Must always be sure that 1282C NP <= NQ or else will need to fix this) 1283C QM for DQDX (derivatives of Q in x direction at cell centers) 1284C DQP for DQDY (derivatives of Q in y direction at cell centers) 1285C DQM for DQDZ (derivatives of Q in z direction at cell centers) 1286C 1287 WRITE (IOUT,1550) 1288 1550 FORMAT (/' ','MAIN : Converting Solution to New Turbulence ', 1289 1 'Model') 1290 CALL CONVRT (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1291 1 KDIM(IBLK,ILVL), 1292 2 IMODLX, RES, W, NQFROM, NPFROM, 1293 3 IMODEL, Q(IQ), PROPS(IP), NQTO, NPTO, 1294 4 S(IS), VOL(IV), 1295 5 QP, FM, QM, DQP, DQM, TAU, F, 1296 6 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), IFDELQ) 1297 ENDIF 1298 135 CONTINUE 1299C 1300C Output message, Increment grid level counter, Return to read next 1301C coarser grid. 1302C Do not read coarser grids if current run does not require it. 1303C 1304 WRITE (IOUT,1570) ILVL 1305 1570 FORMAT (' ','MAIN : Read Restart Grid Level ',I3) 1306 ILVL = ILVL + 1 1307 IF (ILVL .GT. NLVLS) GO TO 139 1308 GO TO 130 1309C 1310C Finished reading all grid levels 1311C 1312 139 CONTINUE 1313 WRITE (IOUT,*) 1314C 1315C Now reset NQ, NF, NP back to NQTO, NFTO, NPTO 1316C 1317 NQ = NQTO 1318 NF = NFTO 1319 NP = NPTO 1320 CLOSE (IRDRST) 1321C 1322C Check that ITTOT + NITS(all mesh sequences) is less than MXHIST 1323C 1324 NITTOT = ITTOT 1325 DO 138 ICFS = 1, NCFS 1326 NITTOT = NITTOT + NITS(ICFS) 1327 138 CONTINUE 1328 IF (NITTOT .GT. MXHIST) THEN 1329 WRITE (IOUT,1580) NITTOT, MXHIST 1330 STOP 1331 ENDIF 1332 ENDIF 1333 1580 FORMAT (' ','MAIN : ERROR-> Length of residual history ', 1334 1 ' (ITTOT + NITS) ',I6, 1335 2 /' ',15X,' exceeds maximum allocated (MXHIST ', 1336 3 'in hstry.h) ',I6, 1337 4 /' ',15X,' Increase MXHIST to be greater than ', 1338 5 'ITTOT + NITS and recompile.', 1339 6 /' ', 8X,'RUN ABORTING!') 1340C 1341C 8. CALCULATE YNORML 1342C a. Calculate y_normal for initial run or if user forces 1343C b. Calculate R(cell centers); 1344C Set y_normal = abs(y_normal) to allow for re-setting of transition 1345C locations 1346C c. Set y_normal = - (y_normal) to specify transition location 1347C d. If RESTART calculate forces if needed for the 2D farfield vortex bc 1348C 1349C 8-a. Calculate y_normal - Do not recalculate if restarting, unless 1350C user forces from input. 1351C Use DTJ as temporary storage for y_normal. 1352C 1353 IF ( (ITURB .GE. ITLMNR .AND. (.NOT. RESTRT)) 1354 1 .OR. YNCALC ) THEN 1355 NXNODE = MXNODE 1356 CALL YNORM (NXLVLS, NXBLKS, NXBCS, NXPPTS, NXNODE, 1357 1 IOFFS, IOFFP, 1358 2 NLVLS, NBLKS, NBCS, IBCDAT, 1359 3 IDIM, JDIM, KDIM, R, PROPS, DTJ) 1360 END IF 1361C 1362C 8-b. Calculate R(Cell Centers); Set y_normal = abs(y_normal) 1363C Set R to contain position vector of cell center for the normal dependent 1364C source terms of the Reynolds stress turbulent models. R(cell centers) 1365C is also needed in the output subroutine. 1366C Use RES as temporary storage of one block, then copy the R(cell centers) 1367C back into R. 1368C 1369 DO 88 ILVL = 1, NLVLS 1370 DO 87 IBLK = 1, NBLKS 1371 IR = IOFFS(IBLK,ILVL) * 3 + 1 1372 CALL RCENTR (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1373 1 KDIM(IBLK,ILVL), R(IR), RES, 1374 2 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK)) 1375 CALL RCCOPY (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1376 1 KDIM(IBLK,ILVL), RES, R(IR)) 1377C 1378 IP = IOFFP(IBLK,ILVL) * NP + 1 1379 NPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1380 1 (KDIM(IBLK,ILVL)+3) 1381 CALL YNABS (NPTS, PROPS(IP)) 1382 87 CONTINUE 1383 88 CONTINUE 1384C 1385C 8-c. Set y_normal = - (y_normal) to specify transition location 1386C Coarsen and set for all mesh levels 1387C N.B. This destroys the information in the ITRSEG array. 1388C 1389 DO ISEG = 1, NTRSEG 1390 IBLK = ITRSEG(1,ISEG) 1391 ILVL = 1 1392 IP = IOFFP(IBLK,ILVL) * NP + 1 1393 CALL YNTRAN (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1394 1 KDIM(IBLK,ILVL), PROPS(IP), ITRSEG(1,ISEG)) 1395 DO ILVL = 2, NLVLS 1396 DO L = 2, 7 1397 ITRSEG(L,ISEG) = ITRSEG(L,ISEG) / 2 + 1 1398 END DO 1399 IP = IOFFP(IBLK,ILVL) * NP + 1 1400 CALL YNTRAN (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1401 1 KDIM(IBLK,ILVL), PROPS(IP), ITRSEG(1,ISEG)) 1402 END DO 1403 END DO 1404C 1405C 8-d. If RESTART then calculate forces if needed for the 2D farfield vortex bc 1406C Calculate forces on the grid level to begin the restart on. 1407C 1408 IF (RESTRT) THEN 1409 ILVL = NCFS 1410 IF (IFFORC .EQ. IYES) THEN 1411 CALL FRCINT (ITURB,NBLKS,IDIM(1,ILVL), 1412 1 JDIM(1,ILVL), KDIM(1,ILVL), 1413 2 IOFF(1,ILVL), IOFFS(1,ILVL), 1414 3 IOFFQ(1,ILVL), IOFFP(1,ILVL), 1415 4 Q, PROPS, S, R, NXBCS, NBCS, 1416 5 IBCDAT(1,1,1,ILVL), FORCE, CL, CD) 1417 ENDIF 1418 ENDIF 1419C 1420C 9. INITIALIZE PROPERTIES 1421C Initialize the PROPERTIES array. 1422C Use QP, QM, F and DQP as temporary arrays. 1423C 1424 DO 150 ILVL = 1, NLVLS 1425 DO 140 IBLK = 1, NBLKS 1426 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1427 IP = IOFFP(IBLK,ILVL) * NP + 1 1428 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 1429 IV = IOFFS(IBLK,ILVL) + 1 1430 IPSTRT = 1 1431 IPEND = IDIM(IBLK,ILVL) + 1 1432 IF (ITURB .GE. ITLMNR) THEN 1433 CALL PRPRTY (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1434 1 KDIM(IBLK,ILVL), IPSTRT, IPEND, MXSECT, 1435 2 Q(IQ), PROPS(IP), S(IS), VOL(IV), 1436 3 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 1437 4 IMODEL, IFDELQ, 1438 5 STEMP, QP, F, FM, QM, DQP, DQM, WORKX) 1439 ENDIF 1440 140 CONTINUE 1441 150 CONTINUE 1442C 1443C 10. BOUNDARY CONDITION INITIALIZATION AND VERIFICATION 1444C Do boundary conditions for all blocks and all levels. Use QP for work. 1445C a. Initialize boundary data to infinity to check input bc data 1446C b. Update boundary conditions 1447C c. Update Cuts 1448C d. Update properties array (synchronize properties with mean flow) 1449C e. Update boundary conditions (ensure turbulent bcs use updated props) 1450C f. Check boundary data to verify input bc data 1451C 1452 DO 170 ILVL = 1, NLVLS 1453 DO 160 IBLK = 1, NBLKS 1454 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1455 IP = IOFFP(IBLK,ILVL) * NP + 1 1456 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 1457 IV = IOFFS(IBLK,ILVL) + 1 1458 IR = IOFFS(IBLK,ILVL) * 3 + 1 1459C 1460C Storage for PROPS 1461C 1462 IF (ITURB .GE. ITLMNR) THEN 1463 NPPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1464 1 (KDIM(IBLK,ILVL)+3) 1465 ELSE 1466 NPPTS = MXPPTS 1467 ENDIF 1468C 1469C 10-a. Initialize boundary data to infinity to check input bc data 1470C Use WORKX for QTEST. 1471C 1472 CALL INITBC (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1473 1 KDIM(IBLK,ILVL), Q(IQ), WORKX, ITURB) 1474C 1475C 10-b. Update boundary conditions 1476C 1477 CALL BC (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1478 1 1, IDIM(IBLK,ILVL)+1, Q(IQ), PROPS(IP), NPPTS, 1479 2 S(IS), VOL(IV), R(IR), 1480 3 NBCS(IBLK), IBCDAT(1,1,IBLK,ILVL), 1481 4 NXPROF, NPROFL, IOFQBC(1,ILVL), QBC, 1482 5 VISCOS(1,IBLK), IMODEL, QP, MXSECT, CL, RCNTR) 1483 160 CONTINUE 1484 170 CONTINUE 1485C 1486C 10-c. Update Cuts 1487C Do block boundary conditions (CUTS) for all levels 1488C 1489 DO 190 ILVL = 1, NLVLS 1490 DO 180 ICUT = 1, NCUTS 1491 IBLK1 = ICUTS( 1,ICUT,ILVL) 1492 IBLK2 = ICUTS(10,ICUT,ILVL) 1493 IQ1 = IOFFQ(IBLK1,ILVL) * NQ + 1 1494 IQ2 = IOFFQ(IBLK2,ILVL) * NQ + 1 1495 CALL CUT (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 1496 1 KDIM(IBLK1,ILVL), Q(IQ1), 1497 2 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 1498 3 KDIM(IBLK2,ILVL), Q(IQ2), ICUTS(1,ICUT,ILVL), NQ) 1499 IF (ITURB .GE. ITLMNR) THEN 1500 IP1 = IOFFP(IBLK1,ILVL) * NP + 1 1501 IP2 = IOFFP(IBLK2,ILVL) * NP + 1 1502 CALL CUT (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 1503 1 KDIM(IBLK1,ILVL), PROPS(IP1), 1504 2 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 1505 3 KDIM(IBLK2,ILVL), PROPS(IP2), 1506 4 ICUTS(1,ICUT,ILVL), NP) 1507 ENDIF 1508 180 CONTINUE 1509 190 CONTINUE 1510C 1511C 10-d,e. Recalculate properties and boundary conditions as properties 1512C are functions of quantities just set in BC and some BCs are 1513C functions of the properties. 1514C 1515 DO 200 ILVL = 1, NLVLS 1516 DO 195 IBLK = 1, NBLKS 1517 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1518 IP = IOFFP(IBLK,ILVL) * NP + 1 1519 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 1520 IV = IOFFS(IBLK,ILVL) + 1 1521 IR = IOFFS(IBLK,ILVL) * 3 + 1 1522C 1523C Storage for PROPS 1524C 1525 IF (ITURB .GE. ITLMNR) THEN 1526 NPPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1527 1 (KDIM(IBLK,ILVL)+3) 1528 ELSE 1529 NPPTS = MXPPTS 1530 ENDIF 1531C 1532 IF (ITURB .GE. ITLMNR) THEN 1533 IPSTRT = 1 1534 IPEND = IDIM(IBLK,ILVL) + 1 1535C 1536C 10-d. Update properties array (synchronize properties with mean flow) 1537C 1538 CALL PRPRTY (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1539 1 KDIM(IBLK,ILVL), IPSTRT, IPEND, MXSECT, 1540 2 Q(IQ), PROPS(IP), S(IS), VOL(IV), 1541 3 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 1542 4 IMODEL, IFDELQ, 1543 5 STEMP, QP, F, FM, QM, DQP, DQM, WORKX) 1544C 1545C 10-e. Update boundary conditions (ensure turbulent bcs use updated props) 1546C 1547 CALL BC (IDIM(IBLK,ILVL),JDIM(IBLK,ILVL),KDIM(IBLK,ILVL), 1548 1 1, IDIM(IBLK,ILVL)+1, Q(IQ), PROPS(IP), NPPTS, 1549 2 S(IS), VOL(IV), R(IR), 1550 3 NBCS(IBLK), IBCDAT(1,1,IBLK,ILVL), 1551 4 NXPROF, NPROFL, IOFQBC(1,ILVL), QBC, 1552 5 VISCOS(1,IBLK), IMODEL, QP, MXSECT, CL, RCNTR) 1553 ENDIF 1554 195 CONTINUE 1555 200 CONTINUE 1556C 1557C 10-f. Check boundary data to verify input bc data 1558C Check each block on each level to verify that all boundary conditions 1559C have been initialized. 1560C 1561 IBCERR = 0 1562 DO 220 ILVL = 1, NLVLS 1563 DO 210 IBLK = 1, NBLKS 1564 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1565 CALL CHCKBC (ILVL, IBLK, IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1566 1 KDIM(IBLK,ILVL), Q(IQ), IBCERR) 1567 210 CONTINUE 1568 220 CONTINUE 1569 IF (IBCERR .GT. 0) THEN 1570 WRITE (IOUT,1400) IBCERR 1571 STOP 1572 ENDIF 1573 1400 FORMAT (' ','MAIN : ERROR-> ',I8,' Unset Boundary Conditions.', 1574 1 /' ', 8X,'RUN ABORTING!') 1575C 1576C 11. MARCHING CALCULATION 1577C If Marching, then call marching procedure 1578C 1579 IF (ISOLVR .EQ. IMARCH) THEN 1580 IBLK = 1 1581 ILVL = 1 1582 ICFS = 1 1583 IQ = IOFFQ (IBLK,ILVL) * NQ + 1 1584 IQN = IOFFQN(IBLK,ILVL) * NQ * NTMLVL + 1 1585 IP = IOFFP (IBLK,ILVL) * NP + 1 1586 IS = IOFFS (IBLK,ILVL) * 4 * 3 + 1 1587 IV = IOFFS (IBLK,ILVL) + 1 1588 IR = IOFFS (IBLK,ILVL) * 3 + 1 1589 ID = IOFF (IBLK,ILVL) + 1 1590C 1591C Temporary(?) fix for the global storage of the LU when marching 1592C 1593 NABC = MXABC / 2 1594C 1595C Storage for QN - Marching shouldn't be unsteady! 1596C 1597 IF (ITIMED .EQ. IUNSTD) THEN 1598 NQNPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1599 1 (KDIM(IBLK,ILVL)+3) 1600 IFL2DQ = .TRUE. 1601 ELSE 1602 NQNPTS = MXQNPT 1603 IFL2DQ = .FALSE. 1604 ENDIF 1605C 1606C Storage for PROPS 1607C 1608 IF (ITURB .GE. ITLMNR) THEN 1609 NPPTS = (IDIM(IBLK,ILVL)+3) * (JDIM(IBLK,ILVL)+3) * 1610 1 (KDIM(IBLK,ILVL)+3) 1611 ELSE 1612 NPPTS = MXPPTS 1613 ENDIF 1614C 1615C Initialize error code 1616C 1617 IERRCD = INO 1618C 1619 CALL MARCH (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1620 1 MXSECT, NABC, Q(IQ), QN(IQN), NQNPTS, PROPS(IP), NPPTS, 1621 2 S(IS), VOL(IV), R(IR), 1622 3 DTJ(ID), RES, FHAT, DQP,DQM,QP,QM,TAU,F,FM,STEMP,WORKX, 1623 4 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 1624 5 NXPROF, NPROFL, IOFQBC(1,ILVL), QBC, 1625 6 ICUTS(1,1,ILVL), NCUTS, R2NORM, ITRSMX, RESMAX, IRESMX, 1626 7 IORDER(1,IBLK,1,1), LIMIT(1,IBLK), RKAP(1,IBLK), 1627 8 IFLUX(IBLK), VISCOS(1,IBLK), IMODEL, IFDELQ, 1628 9 ABC, RM, RTEMP, AT, AI, 1629 A IMSTRT, IMEND, NITS(ICFS), NITFO, NITALG, ITUPJ, 1630 B TOLER, VGNRN, SIGMA, ITER, QMIN, NFIXQ, NRELIZ, 1631 C NPRLIM, NPRNEG, IFL2DQ, R2DQDT, CL, RCNTR, IERRCD) 1632C 1633 IF (IERRCD .NE. INO) THEN 1634 WRITE (IOUT,1405) QMIN(1), NFIXQ(1), QMIN(5), NFIXQ(5) 1635 GO TO 910 1636 ENDIF 1637 1405 FORMAT (///' ','MAIN : FIXQ : ERROR-> Density < ',1PE9.2, 1638 1 ' at',I7,' locations' 1639 2 /' ',' Pressure < ',1PE9.2, 1640 3 ' at',I7,' locations' 1641 4 /' ','RUN ABORTING! Attempting to write output ', 1642 5 'and restart files.', 1643 6 ///' ') 1644 ELSE 1645C 1646C 12. TIME DEPENDENT CALCULATION 1647C 1648C Calculate execution time of run 1649C 1650 1651#ifdef CPU_TIME 1652C 1653 RUNTIME = DTIME (TARRAY) 1654#endif 1655 1656C 1657 MRCHNG = .FALSE. 1658C 1659C Initialize error code 1660C 1661 IERRCD = INO 1662C 1663C Loop over coarse to fine sequences 1664C 1665 ITRUN = 0 1666 DO 900 ICFS = 1, NCFS 1667 LVLONE = NCFS - ICFS + 1 1668C 1669C Loop over number of iterations 1670C 1671 DO 800 ITER = 1, NITS(LVLONE) 1672C 1673C Check for presence of STOP file; if it exists, write restart and quit 1674C 1675 OPEN (ISTOPFL,FILE='STOP',STATUS='OLD',ERR=701) 1676 WRITE (IOUT,1490) 1677 GO TO 905 1678 1490 FORMAT (//' ','MAIN : STOP file encountered! Run stopping.'//) 1679C 1680C STOP does not exist; continue run 1681C 1682 701 CONTINUE 1683C 1684C If unsteady calculation then 1685C i. Update TIME 1686C ii. Collect Q to coarser levels at time level n 1687C iii. Shift Q, QN for next time step 1688C 1689 IF (ITIMED .EQ. IUNSTD) THEN 1690C 1691C i. Update TIME 1692C 1693 TIME = TIME + TIMEDT 1694C 1695C ii. Collect Q to coarser levels at time level n 1696C 1697 DO 710 ILVL = LVLONE+1, LVLONE + MGLVLS(LVLONE) - 1 1698 DO 705 IBLK = 1, NBLKS 1699 IQF = IOFFQ(IBLK,ILVL-1) * NQ + 1 1700 IVF = IOFFS(IBLK,ILVL-1) + 1 1701 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1702 IV = IOFFS(IBLK,ILVL) + 1 1703C 1704C Temporary storage: QP for QS, RES for TKE 1705C 1706 CALL COLQ2H (IDIM(IBLK,ILVL-1), 1707 1 JDIM(IBLK,ILVL-1),KDIM(IBLK,ILVL-1), 1708 2 Q(IQF), VOL(IVF), 1709 3 IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1710 4 KDIM(IBLK,ILVL), Q(IQ), VOL(IV), 1711 5 QP, ITURB, RES) 1712 705 CONTINUE 1713 710 CONTINUE 1714C 1715C iii. Shift Q, QN for next time step 1716C 1717 DO 720 ILVL = LVLONE, LVLONE + MGLVLS(LVLONE) - 1 1718 DO 715 IBLK = 1, NBLKS 1719 IQ = IOFFQ (IBLK,ILVL) * NQ + 1 1720 IQN = IOFFQN(IBLK,ILVL) * NQ * NTMLVL + 1 1721C 1722 CALL QNSHFT (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1723 1 KDIM(IBLK,ILVL), Q(IQ), QN(IQN), 1724 2 QP, QM, IMODEL) 1725 715 CONTINUE 1726 720 CONTINUE 1727 ENDIF 1728C 1729C End of data shift for next time-accurate time step 1730C 1731C Q is Q^(p) and is initialized to Q^(n) to begin subiterations, 1732C QN(1) = Q^(n) and QN(2) = Q^(n-1) 1733C 1734C Perform sub-iterations for each time step 1735C 1736 DO 790 ITSUB = 1, NSUBIT 1737C 1738C Set starting multigrid level to be NCFS for coarsest mesh and 1739C starting level to be 1 for finest mesh and loop over 1740C the multigrid levels. 1741C 1742 DO 700 ILVL = LVLONE, LVLONE + MGLVLS(LVLONE) - 1 1743C 1744C Restrict from fine to coarse level if current level is not finest level 1745C Begin Restriction: 1746C 1747 IF (ILVL .NE. LVLONE) THEN 1748C 1749C Collect Q 1750C 1751 DO 300 IBLK = 1, NBLKS 1752 IQF = IOFFQ(IBLK,ILVL-1) * NQ + 1 1753 IVF = IOFFS(IBLK,ILVL-1) + 1 1754 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1755 IP = IOFFP(IBLK,ILVL) * NP + 1 1756 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 1757 IV = IOFFS(IBLK,ILVL) + 1 1758 IR = IOFFS(IBLK,ILVL) * 3 + 1 1759C 1760 IF (ITURB .GE. ITLMNR) THEN 1761 NPPTS = (IDIM(IBLK,ILVL)+3) * 1762 1 (JDIM(IBLK,ILVL)+3) * 1763 2 (KDIM(IBLK,ILVL)+3) 1764 ELSE 1765 NPPTS = MXPPTS 1766 ENDIF 1767C 1768C Temporary storage: QP for QS, RES for TKE 1769C 1770 CALL COLQ2H (IDIM(IBLK,ILVL-1), 1771 1 JDIM(IBLK,ILVL-1),KDIM(IBLK,ILVL-1), 1772 2 Q(IQF), VOL(IVF), 1773 3 IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1774 4 KDIM(IBLK,ILVL), Q(IQ), VOL(IV), 1775 5 QP, ITURB, RES) 1776C 1777 IPSTRT = 1 1778 IPEND = IDIM(IBLK,ILVL) + 1 1779C 1780 CALL BC (IDIM(IBLK,ILVL),JDIM(IBLK,ILVL), 1781 1 KDIM(IBLK,ILVL), IPSTRT, IPEND, 1782 2 Q(IQ), PROPS(IP), NPPTS, 1783 3 S(IS), VOL(IV), R(IR), 1784 4 NBCS(IBLK), IBCDAT(1,1,IBLK,ILVL), 1785 5 NXPROF,NPROFL,IOFQBC(1,ILVL),QBC, 1786 6 VISCOS(1,IBLK), IMODEL, QP, MXSECT, 1787 7 CL, RCNTR) 1788C 1789 IF (ITURB .GE. ITLMNR) THEN 1790 CALL PRPRTY (IDIM(IBLK,ILVL), 1791 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1792 2 IPSTRT, IPEND, MXSECT, 1793 3 Q(IQ), PROPS(IP), S(IS), VOL(IV), 1794 4 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 1795 5 IMODEL, IFDELQ, 1796 6 STEMP, QP, F, FM, QM, DQP,DQM,WORKX) 1797 ENDIF 1798 300 CONTINUE 1799C 1800C Update CUTS 1801C 1802 DO 320 ICUT = 1, NCUTS 1803 IBLK1 = ICUTS( 1,ICUT,ILVL) 1804 IBLK2 = ICUTS(10,ICUT,ILVL) 1805 IQ1 = IOFFQ(IBLK1,ILVL) * NQ + 1 1806 IQ2 = IOFFQ(IBLK2,ILVL) * NQ + 1 1807 CALL CUT (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 1808 2 KDIM(IBLK1,ILVL), Q(IQ1), 1809 3 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 1810 4 KDIM(IBLK2,ILVL), Q(IQ2), 1811 5 ICUTS(1,ICUT,ILVL), NQ) 1812 IF (ITURB .GE. ITLMNR) THEN 1813 IP1 = IOFFP(IBLK1,ILVL) * NP + 1 1814 IP2 = IOFFP(IBLK2,ILVL) * NP + 1 1815 CALL CUT 1816 1 (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 1817 2 KDIM(IBLK1,ILVL), PROPS(IP1), 1818 3 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 1819 4 KDIM(IBLK2,ILVL), PROPS(IP2), 1820 5 ICUTS(1,ICUT,ILVL), NP) 1821 ENDIF 1822 320 CONTINUE 1823C 1824C Copy QC and store 1825C 1826 DO 340 IBLK = 1, NBLKS 1827 IQ = IOFFQ (IBLK,ILVL) * NQ + 1 1828 IQCRS = IOFFQC(IBLK,ILVL) * NQ + 1 1829C 1830 NPTS = (IDIM(IBLK,ILVL)+3) * 1831 1 (JDIM(IBLK,ILVL)+3) * 1832 2 (KDIM(IBLK,ILVL)+3) * NQ 1833C 1834 CALL CPVCTR (NPTS, Q(IQ), QCRS(IQCRS)) 1835 340 CONTINUE 1836C 1837C Calculate coarse and fine grid residuals and coarse grid forcing function 1838C 1839 NPRLIM = 0 1840 NPRNEG = 0 1841 DO 360 IBLK = 1, NBLKS 1842 IQF = IOFFQ (IBLK,ILVL-1) * NQ + 1 1843 IQNF = IOFFQN(IBLK,ILVL-1) * NQ *NTMLVL + 1 1844 IPF = IOFFP (IBLK,ILVL-1) * NP + 1 1845 ISF = IOFFS (IBLK,ILVL-1) * 4 * 3 + 1 1846 IVF = IOFFS (IBLK,ILVL-1) + 1 1847 IDF = IOFF (IBLK,ILVL-1) + 1 1848 IRF = IOFFS (IBLK,ILVL-1) * 3 + 1 1849C 1850 IQC = IOFFQ (IBLK,ILVL) * NQ + 1 1851 IQNC = IOFFQN(IBLK,ILVL) * NQ *NTMLVL + 1 1852 IPC = IOFFP (IBLK,ILVL) * NP + 1 1853 ISC = IOFFS (IBLK,ILVL) * 4 * 3 + 1 1854 IVC = IOFFS (IBLK,ILVL) + 1 1855 IDC = IOFF (IBLK,ILVL) + 1 1856 IRC = IOFFS (IBLK,ILVL) * 3 + 1 1857C 1858 IRSCRS = IOFFQC(IBLK,ILVL) * NF + 1 1859C 1860 IF (ITIMED .EQ. IUNSTD) THEN 1861 NQNPTF = (IDIM(IBLK,ILVL-1)+3) * 1862 1 (JDIM(IBLK,ILVL-1)+3) * 1863 2 (KDIM(IBLK,ILVL-1)+3) 1864 NQNPTC = (IDIM(IBLK,ILVL) +3) * 1865 1 (JDIM(IBLK,ILVL) +3) * 1866 2 (KDIM(IBLK,ILVL) +3) 1867 ELSE 1868 NQNPTF = MXQNPT 1869 NQNPTC = MXQNPT 1870 ENDIF 1871C 1872 IF (ITURB .GE. ITLMNR) THEN 1873 NPPTSF = (IDIM(IBLK,ILVL-1)+3) * 1874 1 (JDIM(IBLK,ILVL-1)+3) * 1875 2 (KDIM(IBLK,ILVL-1)+3) 1876 NPPTSC = (IDIM(IBLK,ILVL) +3) * 1877 1 (JDIM(IBLK,ILVL) +3) * 1878 2 (KDIM(IBLK,ILVL) +3) 1879 ELSE 1880 NPPTSF = MXPPTS 1881 NPPTSC = MXPPTS 1882 ENDIF 1883C 1884 IF (ILVL .EQ. (LVLONE + 1) ) THEN 1885 INITRS = .TRUE. 1886 ELSE 1887 INITRS = .FALSE. 1888 IRSCRF = IOFFQC(IBLK,ILVL-1) * NF + 1 1889 NPTRSF = (IDIM(IBLK,ILVL-1)+1) * 1890 1 (JDIM(IBLK,ILVL-1)+1) * 1891 2 (KDIM(IBLK,ILVL-1)+1) * NF 1892 CALL CPVCTR (NPTRSF, RESCRS(IRSCRF), RES) 1893 ENDIF 1894 MRCHNG = .FALSE. 1895 IFL2DQ = .FALSE. 1896C 1897 IPSTRT = 2 1898 IPEND = IDIM(IBLK,ILVL-1) 1899 CALL RESID (IDIM(IBLK,ILVL-1), 1900 1 JDIM(IBLK,ILVL-1),KDIM(IBLK,ILVL-1), 1901 2 IPSTRT, IPEND, MXSECT, 1902 3 Q(IQF), QN(IQNF), NQNPTF, 1903 4 PROPS(IPF), NPPTSF, S(ISF), 1904 5 VOL(IVF), R(IRF), DTJ(IDF), RES, 1905 6 FHAT, DQP, DQM, QP, QM, TAU, F, FM, 1906 7 STEMP, WORKX, 1907 8 IBCDAT(1,1,IBLK,ILVL-1), NBCS(IBLK), 1908 9 IORDER(1,IBLK,ILVL-1,ICFS), 1909 A LIMIT(1,IBLK), RKAP(1,IBLK), 1910 B IFLUX(IBLK), VISCOS(1,IBLK), IMODEL, 1911 C INITRS, IFDELQ, MRCHNG, VGNRN,SIGMA, 1912 D NPRLIM, NPRNEG, IFL2DQ, R2DQDT) 1913C 1914C Set logical to calculate multigrid forcing function on first resid call 1915C 1916 IFMGFF = .TRUE. 1917C 1918 CALL COLR2H (IDIM(IBLK,ILVL-1), 1919 1 JDIM(IBLK,ILVL-1),KDIM(IBLK,ILVL-1), 1920 2 RES, 1921 3 IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 1922 4 KDIM(IBLK,ILVL), RESCRS(IRSCRS), QP) 1923 360 CONTINUE 1924 ENDIF 1925C 1926C End Restriction: 1927C 1928C Update time step. 1929C Use DQP as temporary storage for PFACTR. 1930C 1931 IF ( MOD((ITER-1), ITUPDT) .EQ. 0 .AND. 1932 1 ITSUB .EQ. 1 ) THEN 1933 DTMIN = 999.E0 1934 DTMAX = -999.E0 1935 DO 400 IBLK = 1, NBLKS 1936 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 1937 IP = IOFFP(IBLK,ILVL) * NP + 1 1938 ID = IOFF (IBLK,ILVL) + 1 1939 IV = IOFFS(IBLK,ILVL) + 1 1940 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 1941 CALL DELTAT (IBLK, IDIM(IBLK,ILVL), 1942 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 1943 2 2, IDIM(IBLK,ILVL), 1944 3 Q(IQ), PROPS(IP), DTJ(ID), VOL(IV), 1945 4 S(IS), DQP, 1946 5 VISCOS(1,IBLK), ITURB, ITER, CFLCUR, 1947 6 DTMIN,IBLKMN,IDTMIN,JDTMIN,KDTMIN, 1948 7 DTMAX,IBLKMX,IDTMAX,JDTMAX,KDTMAX) 1949 400 CONTINUE 1950C 1951C Output minimum and maximum time step/CFL 1952C 1953 IF (DT .GT. 0.0E0) THEN 1954 WRITE (IOUT,1440) CFLCUR, 1955 1 DTMIN,IBLKMN,IDTMIN,JDTMIN,KDTMIN, 1956 2 DTMAX,IBLKMX,IDTMAX,JDTMAX,KDTMAX 1957 ELSE 1958 WRITE (IOUT,1450) CFLCUR, 1959 1 DTMIN,IBLKMN,IDTMIN,JDTMIN,KDTMIN, 1960 2 DTMAX,IBLKMX,IDTMAX,JDTMAX,KDTMAX 1961 ENDIF 1962 ENDIF 1963 1440 FORMAT (' ','MAIN : DT =',1PE10.3,' CFLMIN =',1PE10.3, 1964 1 ' At BLOCK:I,J,K=(',I3,':',I4,',',I4,',',I4,')', 1965 2 /' ',22X,' CFLMAX =',1PE10.3, 1966 3 ' At BLOCK:I,J,K=(',I3,':',I4,',',I4,',',I4,')') 1967 1450 FORMAT (' ','MAIN : CFL =',1PE10.3,' DTMIN =',1PE10.3, 1968 1 ' At BLOCK:I,J,K=(',I3,':',I4,',',I4,',',I4,')', 1969 2 /' ',23X,' DTMAX =',1PE10.3, 1970 3 ' At BLOCK:I,J,K=(',I3,':',I4,',',I4,',',I4,')') 1971C 1972C Perform ITSLVL iterations (time step integrations) per level. 1973C 1974 DO 600 INRIT = 1, ITSLVL(ILVL,LVLONE) 1975C 1976C Reset norms to zero 1977C 1978 DO 420 L = 1, NF + 1 1979 R2NORM(L) = 0.E0 1980 R2DQDT(L) = 0.E0 1981 420 CONTINUE 1982C 1983C Perform NSTAGE stages of multistage scheme 1984C 1985 DO 550 ISTAGE = 1, NSTAGE 1986C 1987C Reset counters to zero 1988C 1989 DO 425 L = 1, NQ 1990 NFIXQ (L) = 0 1991 NRELIZ(L) = 0 1992 425 CONTINUE 1993C 1994C Counter for limiting of turbulence production 1995C 1996 NPRLIM = 0 1997 NPRNEG = 0 1998C 1999C Set up control parameter for calculating maximum residual 2000C 2001 RESMAX = 0.E0 2002 IF (INRIT .EQ. 1 .AND. 2003 1 ISTAGE .EQ. 1 .AND. 2004 2 MOD((ITER),ITRSMX).EQ.0) THEN 2005 IFRSMX = IYES 2006 ELSE 2007 IFRSMX = INO 2008 ENDIF 2009C 2010C Loop over all blocks. 2011C 2012 DO 500 IBLK = 1, NBLKS 2013 IQ = IOFFQ (IBLK,ILVL) * NQ + 1 2014 IQN = IOFFQN(IBLK,ILVL) * NQ * NTMLVL + 1 2015 IP = IOFFP (IBLK,ILVL) * NP + 1 2016 IS = IOFFS (IBLK,ILVL) * 4 * 3 + 1 2017 IV = IOFFS (IBLK,ILVL) + 1 2018 IR = IOFFS (IBLK,ILVL) * 3 + 1 2019 ID = IOFF (IBLK,ILVL) + 1 2020 IABC = IOFF (IBLK,ILVL)*NRANK*NRANK*3*3+1 2021C 2022C Storage for QN, calculate L2 norm of dQ/dt only if first stage of 2023C multistage scheme and on the finest mesh level 2024C 2025 IF (ITIMED .EQ. IUNSTD) THEN 2026 NQNPTS = (IDIM(IBLK,ILVL)+3) * 2027 1 (JDIM(IBLK,ILVL)+3) * 2028 2 (KDIM(IBLK,ILVL)+3) 2029 IF (ILVL .EQ. LVLONE .AND. 2030 1 ISTAGE .EQ. 1) THEN 2031 IFL2DQ = .TRUE. 2032 ELSE 2033 IFL2DQ = .FALSE. 2034 ENDIF 2035 ELSE 2036 NQNPTS = MXQNPT 2037 IFL2DQ = .FALSE. 2038 ENDIF 2039C 2040C Storage for PROPS 2041C 2042 IF (ITURB .GE. ITLMNR) THEN 2043 NPPTS = (IDIM(IBLK,ILVL)+3) * 2044 1 (JDIM(IBLK,ILVL)+3) * 2045 2 (KDIM(IBLK,ILVL)+3) 2046 ELSE 2047 NPPTS = MXPPTS 2048 ENDIF 2049C 2050C Set up switch from first order scheme to higher order 2051C 2052 ITORDR(1) = IORDER(1,IBLK,ILVL,ICFS) 2053 ITORDR(2) = IORDER(2,IBLK,ILVL,ICFS) 2054 ITORDR(3) = IORDER(3,IBLK,ILVL,ICFS) 2055 IF (ITER .LE. NITFO) THEN 2056 ITORDR(1) = 1 2057 ITORDR(2) = 1 2058 ITORDR(3) = 1 2059 ENDIF 2060C 2061C Set up use of algebraic turbulent eddy viscosity with higher order 2062C models for the initial transient 2063C 2064 IF (ITER .LE. NITALG) THEN 2065 MUTALG = .TRUE. 2066 ELSE 2067 MUTALG = .FALSE. 2068 ENDIF 2069C 2070C If running a viscous case, calculate the PROPS array 2071C Use QP, QM, F and DQP as temporary storage. 2072C 2073 IF (ITURB .GE. ITLMNR .AND. 2074 1 ISTAGE .EQ. 1) THEN 2075 CALL PRPRTY (IDIM(IBLK,ILVL), 2076 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), 2077 2 1, IDIM(IBLK,ILVL)+1, MXSECT, 2078 3 Q(IQ), PROPS(IP), S(IS), VOL(IV), 2079 4 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 2080 5 IMODEL, IFDELQ, 2081 6 STEMP, QP, F, FM, QM, DQP,DQM,WORKX) 2082 ENDIF 2083C 2084C Calculate the Residual for the current block and level 2085C 2086 INITRS = .TRUE. 2087 IPSTRT = 2 2088 IPEND = IDIM(IBLK,ILVL) 2089 CALL RESID (IDIM(IBLK,ILVL), 2090 1 JDIM(IBLK,ILVL),KDIM(IBLK,ILVL), 2091 2 IPSTRT, IPEND, MXSECT, 2092 3 Q(IQ), QN(IQN), NQNPTS, 2093 4 PROPS(IP), NPPTS, S(IS), 2094 5 VOL(IV), R(IR), DTJ(ID), RES, 2095 6 FHAT, DQP, DQM, QP, QM, TAU, F, FM, 2096 7 STEMP, WORKX, 2097 8 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 2098 9 ITORDR, 2099 A LIMIT(1,IBLK), RKAP(1,IBLK), 2100 B IFLUX(IBLK), VISCOS(1,IBLK), IMODEL, 2101 C INITRS, IFDELQ, MRCHNG, VGNRN,SIGMA, 2102 D NPRLIM, NPRNEG, IFL2DQ, R2DQDT) 2103C 2104C Calculate L2 norms only on the finest mesh 2105C If first stage, call l2norm to calculate L2 norm of residual, 2106C call l2max to calculate maximum residual 2107C 2108 IF (ILVL .EQ. LVLONE) THEN 2109 IF (ISTAGE .EQ. 1) THEN 2110 CALL L2NORM (IDIM(IBLK,ILVL), 2111 1 JDIM(IBLK,ILVL),KDIM(IBLK,ILVL), 2112 2 IPSTRT, IPEND, 2113 3 RES, VOL(IV), R2NORM) 2114 IF (IFRSMX .EQ. IYES) THEN 2115 CALL L2MAX (IBLK,IDIM(IBLK,ILVL), 2116 1 JDIM(IBLK,ILVL),KDIM(IBLK,ILVL), 2117 2 IPSTRT, IPEND, 2118 3 RES, VOL(IV), RESMAX, IRESMX) 2119 ENDIF 2120 ENDIF 2121 ENDIF 2122C 2123C For MG: i. For the first residual evaluation after coarsening, 2124C calculate the MG forcing function 2125C ii. Set the residual to RES + FF before time integration 2126C 2127 IF (ILVL .NE. LVLONE) THEN 2128 IRSCRS = IOFFQC(IBLK,ILVL) * NF + 1 2129 NPTS = (IDIM(IBLK,ILVL)+1) * 2130 1 (JDIM(IBLK,ILVL)+1) * 2131 2 (KDIM(IBLK,ILVL)+1) * NF 2132 IF (IFMGFF) THEN 2133c%%%%% CALL CMGFF (NPTS,RESCRS(IRSCRS),RES) 2134 CALL xCMGFF (idim(iblk,ilvl), 2135 1 jdim(iblk,ilvl), kdim(iblk,ilvl), 2136 2 RESCRS(IRSCRS),RES) 2137 ENDIF 2138c%%%%% CALL CRESFF (NPTS, RESCRS(IRSCRS), RES) 2139 CALL xCRESFF (idim(iblk,ilvl), 2140 1 jdim(iblk,ilvl), kdim(iblk,ilvl), 2141 2 RESCRS(IRSCRS),RES) 2142 ENDIF 2143C 2144 IF (ISOLVR .EQ. IRKN) THEN 2145 CALL RK (IBLK, IDIM(IBLK,ILVL), 2146 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL),MXSECT, 2147 2 NSTAGE, ISTAGE, RKALPH(ISTAGE), 2148 3 Q(IQ), QN(IQN), NQNPTS, 2149 4 PROPS(IP), NPPTS, S(IS), VOL(IV), 2150 5 R(IR),DTJ(ID), RES, FHAT, DQP,DQM,QP,QM, 2151 6 TAU, F, FM, STEMP, WORKX, 2152 7 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 2153 8 ITORDR, LIMIT(1,IBLK), RKAP(1,IBLK), 2154 9 IFLUX(IBLK), VISCOS(1,IBLK), 2155 A IMODEL, IFDELQ, W, 2156 B QMIN, NFIXQ, NRELIZ, NPRLIM, NPRNEG, 2157 C IERRCD) 2158C 2159 ELSE IF (ISOLVR .EQ. IAF3F) THEN 2160 CALL AF3F (IBLK, IDIM(IBLK,ILVL), 2161 1 JDIM(IBLK,ILVL), KDIM(IBLK,ILVL), MXSECT, 2162 2 MXABC, Q(IQ), QN(IQN), NQNPTS, 2163 3 PROPS(IP), NPPTS, S(IS), VOL(IV), 2164 4 R(IR), DTJ(ID), RES, FHAT, DQP, DQM, QP,QM, 2165 5 TAU, F, FM, STEMP, WORKX, 2166 6 IBCDAT(1,1,IBLK,ILVL), NBCS(IBLK), 2167 7 ITORDR, LIMIT(1,IBLK), RKAP(1,IBLK), 2168 8 IFLUX(IBLK), VISCOS(1,IBLK), DIAG(1,IBLK), 2169 9 IMODEL, IFDELQ, ABC, RTEMP, AT, AI, 2170 A QMIN, NFIXQ, NRELIZ, NPRLIM, NPRNEG,IERRCD) 2171C 2172 ENDIF 2173C 2174C Update Boundary Conditions 2175C 2176 CALL BC(IDIM(IBLK,ILVL),JDIM(IBLK,ILVL), 2177 1 KDIM(IBLK,ILVL), 1, 2178 2 IDIM(IBLK,ILVL)+1,Q(IQ),PROPS(IP), 2179 3 NPPTS, S(IS), VOL(IV), R(IR), 2180 4 NBCS(IBLK),IBCDAT(1,1,IBLK,ILVL), 2181 5 NXPROF,NPROFL,IOFQBC(1,ILVL),QBC, 2182 6 VISCOS(1,IBLK),IMODEL,QP, MXSECT, 2183 7 CL, RCNTR) 2184 500 CONTINUE 2185C 2186C Reset the logical controlling the calculation of the MG forcing function 2187C 2188 IFMGFF = .FALSE. 2189C 2190C Check error code and abort 2191C 2192 IF (IERRCD .NE. INO) THEN 2193 WRITE (IOUT,1405) QMIN(1), NFIXQ(1), 2194 1 QMIN(5), NFIXQ(5) 2195 GO TO 910 2196 ENDIF 2197C 2198C Do block boundary conditions (CUTS) 2199C 2200 DO 530 ICUT = 1, NCUTS 2201 IBLK1 = ICUTS( 1,ICUT,ILVL) 2202 IBLK2 = ICUTS(10,ICUT,ILVL) 2203 IQ1 = IOFFQ(IBLK1,ILVL) * NQ + 1 2204 IQ2 = IOFFQ(IBLK2,ILVL) * NQ + 1 2205 CALL CUT 2206 1 (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 2207 2 KDIM(IBLK1,ILVL), Q(IQ1), 2208 3 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 2209 4 KDIM(IBLK2,ILVL), Q(IQ2), 2210 5 ICUTS(1,ICUT,ILVL), NQ) 2211 IF (ITURB .GE. ITLMNR) THEN 2212 IP1 = IOFFP(IBLK1,ILVL) * NP + 1 2213 IP2 = IOFFP(IBLK2,ILVL) * NP + 1 2214 CALL CUT 2215 1 (IDIM(IBLK1,ILVL), JDIM(IBLK1,ILVL), 2216 2 KDIM(IBLK1,ILVL), PROPS(IP1), 2217 3 IDIM(IBLK2,ILVL), JDIM(IBLK2,ILVL), 2218 4 KDIM(IBLK2,ILVL), PROPS(IP2), 2219 5 ICUTS(1,ICUT,ILVL), NP) 2220 ENDIF 2221 530 CONTINUE 2222C 2223C End of stage loop 2224C 2225 550 CONTINUE 2226C 2227C Output monitoring data on finest mesh level only 2228C 2229 IF (ILVL .EQ. LVLONE) THEN 2230C 2231C Calculate iteration for output 2232C 2233 ITCUR = ITER + ITRUN 2234 IT = ITCUR + ITTOT 2235C 2236C i. Output monitoring data on subiterations 2237C 2238 IF (ITSUB .NE. 1) THEN 2239C 2240C Output residual data scaled by number of cell centers on current level 2241C 2242 DO 570 LL = 1, NF + 1 2243 R2NORM(LL) = SQRT(R2NORM(LL)) / 2244 1 NBRPTS(ILVL) 2245 R2DQDT(LL) = SQRT(R2DQDT(LL)) / 2246 1 NBRPTS(ILVL) 2247 570 CONTINUE 2248 WRITE (IOUT,2200) ITCUR, ITSUB, 2249 1 (R2NORM(LL),LL=1,NF+1) 2250 IF (ITIMED .EQ. IUNSTD) THEN 2251 WRITE (IOUT,2210) ITCUR, ITSUB, 2252 1 (R2DQDT(LL),LL=1,NF+1) 2253 ENDIF 2254 2200 FORMAT (' ','Iter =',I8,' Subiteration =',I4,8X, 2255 1 ' Res = ',13(1PE14.7,2X)) 2256 2210 FORMAT (' ','Iter =',I8,' Subiteration =',I4,2X, 2257 1 ' L2(dQ/dt) = ',13(1PE14.7,2X)) 2258C 2259C ii. Output monitoring data only on main iterations 2260C 2261 ELSE 2262C 2263C Output residual data scaled by number of cell centers on current level 2264C 2265 DO 580 LL = 1, NF + 1 2266 R2NORM(LL) = SQRT(R2NORM(LL)) / 2267 1 NBRPTS(ILVL) 2268 R2DQDT(LL) = SQRT(R2DQDT(LL)) / 2269 1 NBRPTS(ILVL) 2270 580 CONTINUE 2271 R2(IT) = R2NORM(NF+1) 2272 IF (ITIMED .EQ. ISTDY) THEN 2273 WRITE (IOUT,2300) ITCUR, IT, 2274 1 (R2NORM(LL),LL=1,NF+1) 2275 ELSE 2276 WRITE (IOUT,2305) ITCUR, IT, TIME, 2277 1 (R2NORM(LL),LL=1,NF+1) 2278 WRITE (IOUT,2210) ITCUR, ITSUB, 2279 1 (R2DQDT(LL),LL=1,NF+1) 2280 ENDIF 2281 2300 FORMAT (' ','ITER =',2I8,' RES = ',13(1PE14.7,2X)) 2282 2305 FORMAT (' ','ITER =',2I8,' TIME =',1PE12.5, 2283 1 ' RES = ',13(1PE14.7,2X)) 2284C 2285C Output maximum residual 2286C 2287 IF (IFRSMX .EQ. IYES) THEN 2288 RESMAX = SQRT (RESMAX) 2289 WRITE (IOUT,2350) ITCUR, IT, RESMAX, 2290 1 (IRESMX(L),L=0,3) 2291 ENDIF 2292 2350 FORMAT (' ','Iter =',2I8,' Max Res = ',1PE12.5, 2293 1 ' At BLOCK:I,J,K=(',I3,':',I4,',',I4,',',I4,')') 2294C 2295C Output if reset turbulence production 2296C 2297 IF (NPRNEG .GT. 0) THEN 2298 WRITE (IOUT,2400) NPRNEG 2299 ENDIF 2300 2400 FORMAT (' ','MAIN : SRCRES: Turb. Prod reset posit', 2301 1 'ive at ',I6,' locations') 2302C 2303 IF (NPRLIM .GT. 0) THEN 2304 WRITE (IOUT,2401) PRDLIM, NPRLIM 2305 ENDIF 2306 2401 FORMAT (' ','MAIN : SRCRES: Turb. Prod limited to ',1PE10.2, 2307 1 '*diss at ',I6,' locations') 2308C 2309C Output number of times positivity violated 2310C 2311 DO 585 L = 1, NQ 2312 IF (IFIXQ .EQ. IFIXMN) THEN 2313 FIXTYP = 'Minimum' 2314 ELSE IF (IFIXQ .EQ. IFIXAV) THEN 2315 FIXTYP = 'Average' 2316 ENDIF 2317 IF (NFIXQ(L) .GT. 0) THEN 2318 WRITE (IOUT,2410) L, QMIN(L), 2319 1 NFIXQ(L), FIXTYP 2320 ENDIF 2321 585 CONTINUE 2322 2410 FORMAT (' ','MAIN : FIXQ : Q(',I2,') < ',1PE9.2, 2323 1 ' at ',I7,' locations. Replaced with ',A7,'.') 2324C 2325C Output number of times realizability violated 2326C 2327 IF (ITURB. EQ. ITRS) THEN 2328 IF (NRELIZ( 9) .GT. 0) THEN 2329 WRITE (IOUT,2420) 'XY', NRELIZ( 9) 2330 ENDIF 2331C 2332 IF (NRELIZ(10) .GT. 0) THEN 2333 WRITE (IOUT,2420) 'XZ', NRELIZ(10) 2334 ENDIF 2335C 2336 IF (NRELIZ(11) .GT. 0) THEN 2337 WRITE (IOUT,2420) 'YZ', NRELIZ(11) 2338 ENDIF 2339 ENDIF 2340 2420 FORMAT (' ','MAIN : REALIZ: Tau_',A2,' violates Schwartz ineq. ', 2341 1 ' at ',I6,' locations. Reset.') 2342C 2343C Output Monitor data 2344C 2345 IF ( MOD(ITER, NITMON) .EQ. 0 ) THEN 2346 DO 595 IBLK = 1, NBLKS 2347 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 2348 CALL MONITR 2349 1 (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 2350 2 KDIM(IBLK,ILVL), Q(IQ), 2351 3 NPRMON(IBLK), IPRMON(1,1,IBLK)) 2352 595 CONTINUE 2353 ENDIF 2354 ENDIF 2355C 2356C Calculate forces: update on each subiteration for 2D farfield vortex bc 2357C 2358 IF (IFFORC .EQ. IYES) THEN 2359 CALL FRCINT (ITURB,NBLKS,IDIM(1,ILVL), 2360 1 JDIM(1,ILVL), KDIM(1,ILVL), 2361 2 IOFF(1,ILVL), IOFFS(1,ILVL), 2362 3 IOFFQ(1,ILVL), IOFFP(1,ILVL), 2363 4 Q, PROPS, S, R, NXBCS, NBCS, 2364 5 IBCDAT(1,1,1,ILVL), FORCE, CL, CD) 2365 ENDIF 2366C 2367C Output only on LAST sub-iteration 2368C 2369 IF (ITSUB .EQ. NSUBIT) THEN 2370 IF (IFRCPR .EQ. IYES) THEN 2371 IF (ITIMED .EQ. ISTDY) THEN 2372 WRITE (IOUT,2480) ITCUR, IT, CL, CD 2373 ELSE 2374 WRITE (IOUT,2485) ITCUR, IT, TIME, 2375 1 CL, CD 2376 ENDIF 2377 ENDIF 2378 2480 FORMAT (' ','Iter =',2I8,' CL = ',1PE13.6,' CD = ',1PE13.6) 2379 2485 FORMAT (' ','Iter =',2I8,' Time =',1PE12.5, 2380 1 ' CL = ',1PE13.6,' CD = ',1PE13.6) 2381 ENDIF 2382 ENDIF 2383C 2384C End of output 2385C 2386 600 CONTINUE 2387 700 CONTINUE 2388C 2389C Begin Prolongation from coarser meshes 2390C 1. For each block 2391C a. Interpolate coarse grid correction 2392C b. Update boundary conditions 2393C 2. Update Cuts 2394C 2395 DO 750 ILVL = LVLONE + MGLVLS(LVLONE) - 2, LVLONE, -1 2396 LVLC = ILVL + 1 2397 LVLF = ILVL 2398 DO 730 IBLK = 1, NBLKS 2399 IQC = IOFFQ (IBLK,LVLC) * NQ + 1 2400 IQCRS = IOFFQC(IBLK,LVLC) * NQ + 1 2401 IQF = IOFFQ (IBLK,LVLF) * NQ + 1 2402 IPF = IOFFP (IBLK,LVLF) * NP + 1 2403 ISF = IOFFS (IBLK,LVLF) * 4 * 3 + 1 2404 IVF = IOFFS (IBLK,LVLF) + 1 2405 IRF = IOFFS (IBLK,ILVL) * 3 + 1 2406C 2407C Use RES as temporary storage for DQF (RES should be big enough) 2408C 2409 CALL PROLNG ( 2410 1 IDIM(IBLK,LVLC), JDIM(IBLK,LVLC), 2411 2 KDIM(IBLK,LVLC), Q(IQC), QCRS(IQCRS), 2412 3 NBCS(IBLK), IBCDAT(1,1,IBLK,LVLC), 2413 4 IDIM(IBLK,LVLF), JDIM(IBLK,LVLF), 2414 5 KDIM(IBLK,LVLF), Q(IQF), RES, 2415 6 IBCDAT(1,1,IBLK,LVLF), 2416 7 ITURB, QMIN, NFIXQ, NRELIZ, 2417 8 MXSECT, ABC, RTEMP) 2418C 2419 IF (ITURB .GE. ITLMNR) THEN 2420 NPPTS = (IDIM(IBLK,LVLF)+3) * 2421 1 (JDIM(IBLK,LVLF)+3) * 2422 2 (KDIM(IBLK,LVLF)+3) 2423 ELSE 2424 NPPTS = MXPPTS 2425 ENDIF 2426 CALL BC ( 2427 1 IDIM(IBLK,LVLF),JDIM(IBLK,LVLF), 2428 2 KDIM(IBLK,LVLF), 1, IDIM(IBLK,LVLF)+1, 2429 3 Q(IQF),PROPS(IPF),NPPTS,S(ISF),VOL(IVF), 2430 4 R(IRF),NBCS(IBLK),IBCDAT(1,1,IBLK,LVLF), 2431 5 NXPROF, NPROFL, IOFQBC(1,LVLF), QBC, 2432 6 VISCOS(1,IBLK), IMODEL, QP, MXSECT, 2433 7 CL, RCNTR) 2434C 2435 IF (ITURB .GE. ITLMNR) THEN 2436 CALL PRPRTY (IDIM(IBLK,LVLF), 2437 1 JDIM(IBLK,LVLF), KDIM(IBLK,LVLF), 2438 2 1, IDIM(IBLK,LVLF)+1, MXSECT, 2439 3 Q(IQF), PROPS(IPF), S(ISF), VOL(IVF), 2440 4 IBCDAT(1,1,IBLK,LVLF), NBCS(IBLK), 2441 5 IMODEL, IFDELQ, 2442 6 STEMP, QP, F, FM, QM, DQP, DQM, WORKX) 2443 ENDIF 2444 730 CONTINUE 2445C 2446C Update CUTS 2447C 2448 DO 740 ICUT = 1, NCUTS 2449 IBLK1 = ICUTS( 1,ICUT,LVLF) 2450 IBLK2 = ICUTS(10,ICUT,LVLF) 2451 IQ1 = IOFFQ(IBLK1,LVLF) * NQ + 1 2452 IQ2 = IOFFQ(IBLK2,LVLF) * NQ + 1 2453 CALL CUT (IDIM(IBLK1,LVLF), JDIM(IBLK1,LVLF), 2454 1 KDIM(IBLK1,LVLF), Q(IQ1), 2455 2 IDIM(IBLK2,LVLF), JDIM(IBLK2,LVLF), 2456 3 KDIM(IBLK2,LVLF), Q(IQ2), 2457 4 ICUTS(1,ICUT,LVLF), NQ) 2458 IF (ITURB .GE. ITLMNR) THEN 2459 IP1 = IOFFP(IBLK1,LVLF) * NP + 1 2460 IP2 = IOFFP(IBLK2,LVLF) * NP + 1 2461 CALL CUT (IDIM(IBLK1,LVLF), JDIM(IBLK1,LVLF), 2462 1 KDIM(IBLK1,LVLF), PROPS(IP1), 2463 2 IDIM(IBLK2,LVLF), JDIM(IBLK2,LVLF), 2464 3 KDIM(IBLK2,LVLF), PROPS(IP2), 2465 4 ICUTS(1,ICUT,LVLF), NP) 2466 ENDIF 2467 740 CONTINUE 2468 750 CONTINUE 2469C 2470C End Prolongation 2471C 2472 790 CONTINUE 2473C 2474C End Sub-iterations; Q = Q^(n+1) 2475C 2476C Output intermediate restart file - but not on the last iteration 2477C 2478 IF (ITER/NITRST*NITRST .EQ. ITER .AND. 2479 1 ITER .NE. NITS(LVLONE)) THEN 2480 CALL WRREST (NXBLKS, NXLVLS, NXPTS, NXQNPT, NXPPTS, 2481 1 NXBCS, NXSIZE, NXSECT, 2482 2 IOFFS, IOFFQ, IOFFP, IOFFQN, 2483 3 NBLKS, NLVLS, IMODEL, IFDELQ, ITER, 2484 4 IDIM, JDIM, KDIM, 2485 5 Q, QN, PROPS, S, VOL, NBCS, IBCDAT, 2486 6 RES, QP, QM, DQP, DQM, F, FM,WORKX,STEMP, 2487 7 FWRRST, IERRCD) 2488 IF (IERRCD .NE. INO) GO TO 999 2489 ENDIF 2490 800 CONTINUE 2491C 2492C 12-c. If not finest grid, then interpolate to the next finer grid and 2493C update the properties and boundary conditions. 2494C 2495 IF (ICFS .NE. NCFS) THEN 2496 LVLC = NCFS - ICFS + 1 2497 LVLF = LVLC - 1 2498 WRITE (IOUT,8100) LVLC, LVLF 2499 8100 FORMAT (//' ','MAIN : Interpolating from grid level ',I3, 2500 1 ' to grid level ',I3//) 2501C 2502C 12-c-i. Interpolate coarse grid to fine grid for each block. 2503C Modify cell face bc data to be at ghost cells then call 2504C interpolation routine. 2505C 2506 DO 850 IBLK = 1, NBLKS 2507 IQC = IOFFQ(IBLK,LVLC) * NQ + 1 2508 IQF = IOFFQ(IBLK,LVLF) * NQ + 1 2509 CALL QBCMOD (IDIM(IBLK,LVLC), JDIM(IBLK,LVLC), 2510 1 KDIM(IBLK,LVLC), Q(IQC), 2511 2 NBCS(IBLK), IBCDAT(1,1,IBLK,LVLC)) 2512 CALL INT2HQ (IDIM(IBLK,LVLC), JDIM(IBLK,LVLC), 2513 1 KDIM(IBLK,LVLC), Q(IQC), 2514 2 IDIM(IBLK,LVLF), JDIM(IBLK,LVLF), 2515 3 KDIM(IBLK,LVLF), Q(IQF)) 2516 CALL QBCFIX (IDIM(IBLK,LVLC), JDIM(IBLK,LVLC), 2517 1 KDIM(IBLK,LVLC), Q(IQC), 2518 2 NBCS(IBLK), IBCDAT(1,1,IBLK,LVLC)) 2519 850 CONTINUE 2520C 2521C 12-c-ii. Update boundary conditions and properties on each block 2522C 2523 DO 870 IBLK = 1, NBLKS 2524 IQ = IOFFQ(IBLK,LVLF) * NQ + 1 2525 IP = IOFFP(IBLK,LVLF) * NP + 1 2526 IS = IOFFS(IBLK,LVLF) * 4 * 3 + 1 2527 IV = IOFFS(IBLK,LVLF) + 1 2528 IR = IOFFS(IBLK,LVLF) * 3 + 1 2529C 2530 IF (ITURB .GE. ITLMNR) THEN 2531 NPPTS = (IDIM(IBLK,LVLF)+3) * (JDIM(IBLK,LVLF)+3) * 2532 1 (KDIM(IBLK,LVLF)+3) 2533 IPSTRT = 1 2534 IPEND = IDIM(IBLK,LVLF) + 1 2535C 2536C Update properties array 2537C 2538 CALL PRPRTY (IDIM(IBLK,LVLF), JDIM(IBLK,LVLF), 2539 1 KDIM(IBLK,LVLF), IPSTRT, IPEND,MXSECT, 2540 2 Q(IQ), PROPS(IP), S(IS), VOL(IV), 2541 3 IBCDAT(1,1,IBLK,LVLF), NBCS(IBLK), 2542 4 IMODEL, IFDELQ, 2543 5 STEMP, QP, F, FM, QM, DQP, DQM, WORKX) 2544 ELSE 2545 NPPTS = MXPPTS 2546 ENDIF 2547C 2548C Update boundary conditions 2549C 2550 CALL BC (IDIM(IBLK,LVLF),JDIM(IBLK,LVLF), 2551 1 KDIM(IBLK,LVLF), 1, IDIM(IBLK,LVLF)+1, 2552 2 Q(IQ), PROPS(IP), NPPTS, S(IS), VOL(IV), R(IR), 2553 3 NBCS(IBLK), IBCDAT(1,1,IBLK,LVLF), 2554 4 NXPROF, NPROFL, IOFQBC(1,LVLF), QBC, 2555 5 VISCOS(1,IBLK), IMODEL, QP, MXSECT, CL, RCNTR) 2556 870 CONTINUE 2557C 2558C 12-c-iii. Do block boundary conditions (CUTS) 2559C 2560 DO 890 ICUT = 1, NCUTS 2561 IBLK1 = ICUTS( 1,ICUT,LVLF) 2562 IBLK2 = ICUTS(10,ICUT,LVLF) 2563 IQ1 = IOFFQ(IBLK1,LVLF) * NQ + 1 2564 IQ2 = IOFFQ(IBLK2,LVLF) * NQ + 1 2565 CALL CUT (IDIM(IBLK1,LVLF), JDIM(IBLK1,LVLF), 2566 1 KDIM(IBLK1,LVLF), Q(IQ1), 2567 2 IDIM(IBLK2,LVLF), JDIM(IBLK2,LVLF), 2568 3 KDIM(IBLK2,LVLF), Q(IQ2), 2569 4 ICUTS(1,ICUT,LVLF), NQ) 2570 IF (ITURB .GE. ITLMNR) THEN 2571 IP1 = IOFFP(IBLK1,LVLF) * NP + 1 2572 IP2 = IOFFP(IBLK2,LVLF) * NP + 1 2573 CALL CUT (IDIM(IBLK1,LVLF), JDIM(IBLK1,LVLF), 2574 1 KDIM(IBLK1,LVLF), PROPS(IP1), 2575 2 IDIM(IBLK2,LVLF), JDIM(IBLK2,LVLF), 2576 3 KDIM(IBLK2,LVLF), PROPS(IP2), 2577 4 ICUTS(1,ICUT,LVLF), NP) 2578 ENDIF 2579 890 CONTINUE 2580 ENDIF 2581 ITRUN = ITRUN + NITS(LVLONE) 2582 900 CONTINUE 2583C 2584C Jump here if STOP file is encountered 2585C 2586 905 CONTINUE 2587C 2588C Calculate execution time of run 2589C 2590 2591#ifdef CPU_TIME 2592C 2593 RUNTIME = DTIME (TARRAY) 2594C 2595 WRITE (IOUT,8200) RUNTIME 2596 8200 FORMAT (//' ','MAIN : Approximate solver execution time: ', 2597 1 1PE12.4,' seconds ') 2598C 2599#endif 2600 2601C 2602 ENDIF 2603C 2604C 13. OUTPUT RESTART DATA 2605C Update properties. 2606C If unsteady, convert QN to primitive (N.B. storage for tke is over nnpts) 2607C Write restart file. 2608C 2609 910 CONTINUE 2610C 2611 CALL WRREST (NXBLKS, NXLVLS, NXPTS, NXQNPT, NXPPTS, 2612 1 NXBCS, NXSIZE, NXSECT, 2613 2 IOFFS, IOFFQ, IOFFP, IOFFQN, 2614 3 NBLKS, NLVLS, IMODEL, IFDELQ, ITRUN, 2615 4 IDIM, JDIM, KDIM, 2616 5 Q, QN, PROPS, S, VOL, NBCS, IBCDAT, 2617 6 RES, QP, QM, DQP, DQM, F, FM, WORKX, STEMP, 2618 7 FWRRST, IERRCD) 2619 IF (IERRCD .NE. INO) GO TO 999 2620C 2621C 14. OUTPUT 2622C Write output and PLOT3D output 2623C (PLOT3D output in multiple grid format) 2624C Use RES array as temporary space to interpolate the grid 2625C nodal values to cell centers. 2626C 2627 ILVL = 1 2628 IF ( IFPLT3 ) THEN 2629C 2630C Binary PLOT3D output 2631C 2632 IF ( BINPLT ) THEN 2633C Grid 2634 FERROR = FPLT3G 2635 OPEN (IPLT3G,FILE=FPLT3G,FORM='UNFORMATTED',ERR=990) 2636 WRITE (IPLT3G) NBLKS 2637 IF ( THREED ) THEN 2638 WRITE (IPLT3G) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2639 1 KDIM(IBLK,ILVL)+1,IBLK=1,NBLKS) 2640 ELSE 2641 WRITE (IPLT3G) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2642 1 IBLK=1,NBLKS) 2643 ENDIF 2644C Solution 2645 FERROR = FPLT3Q 2646 OPEN (IPLT3Q,FILE=FPLT3Q,FORM='UNFORMATTED',ERR=990) 2647 WRITE (IPLT3Q) NBLKS 2648 IF ( THREED ) THEN 2649 WRITE (IPLT3Q) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2650 1 KDIM(IBLK,ILVL)+1,IBLK=1,NBLKS) 2651 ELSE 2652 WRITE (IPLT3Q) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2653 1 IBLK=1,NBLKS) 2654 ENDIF 2655C Function file 2656 IF (NQ .GT. 5) THEN 2657 FERROR = FPLT3F 2658 NVAR = (NQ - 5) + NP 2659 OPEN (IPLT3F,FILE=FPLT3F,FORM='UNFORMATTED',ERR=990) 2660 WRITE (IPLT3F) NBLKS 2661 IF ( THREED ) THEN 2662 WRITE (IPLT3F) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2663 1 KDIM(IBLK,ILVL)+1,NVAR,IBLK=1,NBLKS) 2664 ELSE 2665 WRITE (IPLT3F) (IDIM(IBLK,ILVL)+1,JDIM(IBLK,ILVL)+1, 2666 1 NVAR,IBLK=1,NBLKS) 2667 ENDIF 2668 ENDIF 2669C 2670C Formatted PLOT3D output 2671C 2672 ELSE 2673C Grid 2674 FERROR = FPLT3G 2675 OPEN (IPLT3G,FILE=FPLT3G,FORM='FORMATTED',ERR=990) 2676 WRITE (IPLT3G,'(I5)') NBLKS 2677 IF ( THREED ) THEN 2678 WRITE (IPLT3G,'(3I5)') (IDIM(IBLK,ILVL)+1, 2679 1 JDIM(IBLK,ILVL)+1, 2680 2 KDIM(IBLK,ILVL)+1,IBLK=1,NBLKS) 2681 ELSE 2682 WRITE (IPLT3G,'(3I5)') (IDIM(IBLK,ILVL)+1, 2683 1 JDIM(IBLK,ILVL)+1, 2684 2 IBLK=1,NBLKS) 2685 ENDIF 2686C Solution 2687 FERROR = FPLT3Q 2688 OPEN (IPLT3Q,FILE=FPLT3Q,FORM='FORMATTED',ERR=990) 2689 WRITE (IPLT3Q,'(I5)') NBLKS 2690 IF ( THREED ) THEN 2691 WRITE (IPLT3Q,'(3I5)') (IDIM(IBLK,ILVL)+1, 2692 1 JDIM(IBLK,ILVL)+1, 2693 2 KDIM(IBLK,ILVL)+1,IBLK=1,NBLKS) 2694 ELSE 2695 WRITE (IPLT3Q,'(3I5)') (IDIM(IBLK,ILVL)+1, 2696 1 JDIM(IBLK,ILVL)+1, 2697 2 IBLK=1,NBLKS) 2698 ENDIF 2699C Function file 2700 IF (NQ .GT. 5) THEN 2701 FERROR = FPLT3F 2702 NVAR = (NQ - 5) + NP 2703 OPEN (IPLT3F,FILE=FPLT3F,FORM='FORMATTED',ERR=990) 2704 WRITE (IPLT3F,'(I5)') NBLKS 2705 IF ( THREED ) THEN 2706 WRITE (IPLT3F,'(4I5)') (IDIM(IBLK,ILVL)+1, 2707 1 JDIM(IBLK,ILVL)+1, 2708 2 KDIM(IBLK,ILVL)+1, 2709 3 NVAR,IBLK=1,NBLKS) 2710 ELSE 2711 WRITE (IPLT3F,'(4I5)') (IDIM(IBLK,ILVL)+1, 2712 1 JDIM(IBLK,ILVL)+1, 2713 2 NVAR,IBLK=1,NBLKS) 2714 ENDIF 2715 ENDIF 2716 ENDIF 2717 ENDIF 2718C 2719C Plot3d function name file 2720C 2721 IF (NQ .GT. 5) THEN 2722 FERROR = FPLT3FN 2723 OPEN (IPLT3FN,FILE=FPLT3FN,FORM='FORMATTED',ERR=990) 2724 CALL PLOT3DFN (ITURB) 2725 CLOSE (IPLT3FN) 2726 ENDIF 2727C 2728C Finished writing PLOT3D headers 2729C 2730 DO 980 IBLK = 1, NBLKS 2731 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 2732 IP = IOFFP(IBLK,ILVL) * NP + 1 2733 IR = IOFFS(IBLK,ILVL) * 3 + 1 2734 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 2735 IV = IOFFS(IBLK,ILVL) + 1 2736C 2737 INCP = (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 2738 1 (KDIM(IBLK,ILVL)+3) 2739 IP2 = IP + INCP 2740 IP3 = IP2 + INCP 2741C 2742 CALL OUTPUT (IDIM(IBLK,ILVL),JDIM(IBLK,ILVL), 2743 1 KDIM(IBLK,ILVL),Q(IQ),PROPS(IP),PROPS(IP2), 2744 2 PROPS(IP3),ITURB, 2745 3 R(IR),IBLK,NPRINT(IBLK),IPRINT(1,1,IBLK)) 2746 IF ( IFPLT3 ) THEN 2747 IF ( BINPLT ) THEN 2748 WRITE (IPLT3Q) FSMACH, ALPHA, RE, TIME 2749 ELSE 2750 WRITE (IPLT3Q,'(4E16.8)') FSMACH, ALPHA, RE, TIME 2751 ENDIF 2752 DO 959 ICUT = 1, NCUTS 2753 IF (ICUTS(1,ICUT,ILVL) .EQ. IBLK) THEN 2754 CALL CUT3D (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 2755 1 KDIM(IBLK,ILVL), Q(IQ), 2756 2 ICUTS(1,ICUT,ILVL)) 2757 ENDIF 2758 959 CONTINUE 2759 CALL PLOT3D (IDIM(IBLK,ILVL),JDIM(IBLK,ILVL), 2760 1 KDIM(IBLK,ILVL),Q(IQ),PROPS(IP),R(IR),BINPLT) 2761 ENDIF 2762 980 CONTINUE 2763C 2764C Output skin friction and heat transfer data 2765C 2766 DO 982 IBLK = 1, NBLKS 2767 IQ = IOFFQ(IBLK,ILVL) * NQ + 1 2768 IP = IOFFP(IBLK,ILVL) * NP + 1 2769 IR = IOFFS(IBLK,ILVL) * 3 + 1 2770 IS = IOFFS(IBLK,ILVL) * 4 * 3 + 1 2771 IV = IOFFS(IBLK,ILVL) + 1 2772C 2773 INCP = (IDIM(IBLK,ILVL)+3)*(JDIM(IBLK,ILVL)+3)* 2774 1 (KDIM(IBLK,ILVL)+3) 2775 IP2 = IP + INCP 2776 IP3 = IP2 + INCP 2777C 2778 DO 981 IBC = 1, NBCS(IBLK) 2779 IF ( (IBCDAT(1,IBC,IBLK,ILVL) .EQ. IWALL .OR. 2780 1 IBCDAT(1,IBC,IBLK,ILVL) .EQ. IWALFN ) ) THEN 2781 CALL SFHT (IDIM(IBLK,ILVL), JDIM(IBLK,ILVL), 2782 1 KDIM(IBLK,ILVL), 2, IDIM(IBLK,ILVL),ITURB, 2783 2 Q(IQ), PROPS(IP), S(IS), VOL(IV), 2784 3 IBCDAT(1,IBC,IBLK,ILVL), R(IR), IBLK) 2785 ENDIF 2786 981 CONTINUE 2787 982 CONTINUE 2788C 2789 GO TO 999 2790C 2791C 15. Error Handling Section 2792C 15-a. Bad file name 2793C 2794 990 CONTINUE 2795 WRITE (IOUT,9990) FERROR 2796 9990 FORMAT (/' ','MAIN : ERROR-> Error opening file ',A80, 2797 1 /' ', 8X,'RUN ABORTING!') 2798 GO TO 999 2799C 2800C 16. Close output and finish 2801C 2802 999 CONTINUE 2803 CLOSE (IOUT) 2804C 2805C Finished with ISAAC run. 2806C 2807 STOP 2808 END 2809