ADflow  v1.0
ADflow is a finite volume RANS solver tailored for gradient-based aerodynamic design optimization.
residuals.F90
Go to the documentation of this file.
1 module residuals
2 contains
3 
4  subroutine residual_block
5  !
6  ! residual computes the residual of the mean flow equations on
7  ! the current MG level.
8  !
9  use blockpointers
10  use cgnsgrid
11  use flowvarrefstate
12  use inputiteration
13  use inputdiscretization
14  use inputtimespectral
15  use inputunsteady ! Added by HDN
16  use iteration
17  use inputadjoint
18  use flowutils, only: computespeedofsoundsquared, allnodalgradients
19  use fluxes
20 #ifndef USE_TAPENADE
21  use aleutils, only: interplevelale_block, recoverlevelale_block
22 #endif
23  implicit none
24  !
25  ! Local variables.
26  !
27  integer(kind=intType) :: discr
28  integer(kind=intType) :: i, j, k, l
29  integer(kind=intType) :: iale, jale, kale, lale, male ! For loops of ALE
30  real(kind=realtype), parameter :: k1 = 1.05_realtype
31  ! The line below is only used for the low-speed preconditioner part of this routine
32  real(kind=realtype), parameter :: k2 = 0.6_realtype ! Random given number
33 
34  real(kind=realtype), parameter :: m0 = 0.2_realtype ! Mach number preconditioner activation
35  real(kind=realtype), parameter :: alpha = 0_realtype
36  real(kind=realtype), parameter :: delta = 0_realtype
37  !real(kind=realType), parameter :: hinf = 2_realType ! Test phase
38  real(kind=realtype), parameter :: cpres = 4.18_realtype ! Test phase
39  real(kind=realtype), parameter :: temp = 297.15_realtype
40 
41  !
42  ! Local variables
43  !
44  real(kind=realtype) :: k3, h, velxrho, velyrho, velzrho, sos, hinf
45  real(kind=realtype) :: resm, a11, a12, a13, a14, a15, a21, a22, a23, a24, a25, a31, a32, a33, a34, a35
46  real(kind=realtype) :: a41, a42, a43, a44, a45, a51, a52, a53, a54, a55, b11, b12, b13, b14, b15
47  real(kind=realtype) :: b21, b22, b23, b24, b25, b31, b32, b33, b34, b35
48  real(kind=realtype) :: b41, b42, b43, b44, b45, b51, b52, b53, b54, b55
49  real(kind=realtype) :: rhohdash, betamr2
50  real(kind=realtype) :: g, q
51  real(kind=realtype) :: b1, b2, b3, b4, b5
52  real(kind=realtype) :: dwo(nwf)
53  logical :: fineGrid
54 
55  ! Set the value of rFil, which controls the fraction of the old
56  ! dissipation residual to be used. This is only for the runge-kutta
57  ! schemes; for other smoothers rFil is simply set to 1.0.
58  ! Note the index rkStage+1 for cdisRK. The reason is that the
59  ! residual computation is performed before rkStage is incremented.
60 
61  if (smoother == rungekutta) then
62  rfil = cdisrk(rkstage + 1)
63  else
64  rfil = one
65  end if
66 
67  ! Set the value of the discretization, depending on the grid level,
68  ! and the logical fineGrid, which indicates whether or not this
69  ! is the finest grid level of the current mg cycle.
70 
71  discr = spacediscrcoarse
72  if (currentlevel == 1) discr = spacediscr
73 
74  finegrid = .false.
75  if (currentlevel == groundlevel) finegrid = .true.
76 
77  ! ===========================================================
78  !
79  ! Assuming ALE has nothing to do with MG
80  ! The geometric data will be interpolated if in MD mode
81  !
82  ! ===========================================================
83 #ifndef USE_TAPENADE
84  call interplevelale_block
85 #endif
86  ! ===========================================================
87  !
88  ! The fluxes are calculated as usual
89  !
90  ! ===========================================================
91 
92  call inviscidcentralflux
93 
94  select case (discr)
95 
96  case (dissscalar) ! Standard scalar dissipation scheme.
97 
98  if (finegrid) then
99  if (.not. lumpeddiss) then
100  call invisciddissfluxscalar
101  else
102  call invisciddissfluxscalarapprox
103  end if
104  else
105 #ifndef USE_TAPENADE
106  call invisciddissfluxscalarcoarse
107 #endif
108  end if
109 
110  !===========================================================
111 
112  case (dissmatrix) ! Matrix dissipation scheme.
113 
114  if (finegrid) then
115  if (.not. lumpeddiss) then
116  call invisciddissfluxmatrix
117  else
118  call invisciddissfluxmatrixapprox
119  end if
120  else
121 #ifndef USE_TAPENADE
122  call invisciddissfluxmatrixcoarse
123 #endif
124  end if
125 
126  !===========================================================
127 
128  case (upwind) ! Dissipation via an upwind scheme.
129 
130  call inviscidupwindflux(finegrid)
131 
132  end select
133 
134  !-------------------------------------------------------
135  ! Lastly, recover the old s[I,J,K], sFace[I,J,K]
136  ! This shall be done before difussive and source terms
137  ! are computed.
138  !-------------------------------------------------------
139 #ifndef USE_TAPENADE
140  call recoverlevelale_block
141 #endif
142 
143  if (viscous) then
144  ! Only compute viscous fluxes if rFil > 0
145  if (abs(rfil) > thresholdreal) then
146  ! not lumpedDiss means it isn't the PC...call the vicousFlux
147  if (.not. lumpeddiss) then
148  call computespeedofsoundsquared
149  call allnodalgradients
150  call viscousflux
151  else
152  ! This is a PC calc...only include viscous fluxes if viscPC
153  ! is used
154  ! if full visc is true, also need full viscous terms, even if
155  ! lumpedDiss is true
156  call computespeedofsoundsquared
157  if (viscpc) then
158  call allnodalgradients
159  call viscousflux
160  else
161  call viscousfluxapprox
162  end if
163  end if
164  end if
165  end if
166 
167  !===========================================================
168 
169  ! Add the dissipative and possibly viscous fluxes to the
170  ! Euler fluxes. Loop over the owned cells and add fw to dw.
171  ! Also multiply by iblank so that no updates occur in holes
172  if (lowspeedpreconditioner) then
173  do k = 2, kl
174  do j = 2, jl
175  do i = 2, il
176  ! Compute speed of sound
177  sos = sqrt(gamma(i, j, k) * p(i, j, k) / w(i, j, k, irho))
178 
179  ! Compute velocities without rho from state vector
180  ! (w is pointer.. see type blockType setup in block.F90)
181  ! w(0:ib,0:jb,0:kb,1:nw) is allocated in block.F90
182  ! these are per definition nw=[rho,u,v,w,rhoeE]
183  ! so the velocity is simply just taken out below...
184  ! we do not have to divide with rho since it is already
185  ! without rho...
186  velxrho = w(i, j, k, ivx) ! ivx: l. 60 in constants.F90
187  velyrho = w(i, j, k, ivy)
188  velzrho = w(i, j, k, ivz)
189 
190  q = (velxrho**2 + velyrho**2 + velzrho**2)
191 
192  resm = sqrt(q) / sos
193  ! resM above is used as M_a (thesis) and M (paper 2015)
194  ! and is the Free Stream Mach number
195 
196  ! see routine setup above:
197  ! l. 30: real(kind=realType), parameter :: K1 = 1.05_realType
198  ! Random given number for K2:
199  ! l. 31: real(kind=realType), parameter :: K2 = 0.6_realType
200  ! Mach number preconditioner activation for K3:
201  ! l. 32: real(kind=realType), parameter :: M0 = 0.2_realType
202  !
203  ! Compute K3
204  ! eq. 2.7 in Garg 2015. K1, M0 and resM are scalars
205  !
206  ! unfortunately, Garg has switched the K1 and K3 here in the
207  ! code. In both paper and thesis it is K3 that is used to det-
208  ! ermine K1 below
209 
210  !
211  ! Compute K3
212 
213  k3 = k1 * (1 + ((1 - k1 * m0**2) * resm**2) / (k1 * m0**4))
214  ! Compute BetaMr2
215  ! betaMr2 -> eq. 7 in Garg 2015
216  ! (use eq. 2.6 in thesis thesis since paper has an error)
217  ! where a==SoS
218  !
219  ! again, K1 and K3 are switched compared with paper/thesis
220  ! Compute BetaMr2
221  betamr2 = min(max(k3 * (velxrho**2 + velyrho**2 &
222  + velzrho**2), ((k2) * (winf(ivx)**2 &
223  + winf(ivy)**2 + winf(ivz)**2))), sos**2)
224 
225  ! above, the wInf is the free stream velocity
226  !
227  ! Should this first line's first element have SoS^4 or SoS^2
228 
229  a11 = (betamr2) * (1 / sos**4)
230  a12 = zero
231  a13 = zero
232  a14 = zero
233  a15 = (-betamr2) / sos**4
234 
235  a21 = one * velxrho / sos**2
236  a22 = one * w(i, j, k, irho)
237  a23 = zero
238  a24 = zero
239  a25 = one * (-velxrho) / sos**2
240 
241  a31 = one * velyrho / sos**2
242  a32 = zero
243  a33 = one * w(i, j, k, irho)
244  a34 = zero
245  a35 = one * (-velyrho) / sos**2
246 
247  a41 = one * velzrho / sos**2
248  a42 = zero
249  a43 = zero
250  a44 = one * w(i, j, k, irho)
251  a45 = zero + one * (-velzrho) / sos**2
252 
253  ! mham: seems he fixed the above line an irregular way?
254 
255  a51 = one * ((1 / (gamma(i, j, k) - 1)) + (resm**2) / 2)
256  a52 = one * w(i, j, k, irho) * velxrho
257  a53 = one * w(i, j, k, irho) * velyrho
258  a54 = one * w(i, j, k, irho) * velzrho
259  a55 = one * ((-(resm**2)) / 2)
260 
261  b11 = a11 * (gamma(i, j, k) - 1) * q / 2 + a12 * (-velxrho) &
262  / w(i, j, k, irho) + a13 * (-velyrho) / w(i, j, k, irho) + &
263  a14 * (-velzrho) / w(i, j, k, irho) &
264  + a15 * (((gamma(i, j, k) - 1) * q / 2) - sos**2)
265  b12 = a11 * (1 - gamma(i, j, k)) * velxrho + a12 * 1 / w(i, j, k, irho) &
266  + a15 * (1 - gamma(i, j, k)) * velxrho
267  b13 = a11 * (1 - gamma(i, j, k)) * velyrho + a13 &
268  / w(i, j, k, irho) + a15 * (1 - gamma(i, j, k)) * velyrho
269  b14 = a11 * (1 - gamma(i, j, k)) * velzrho &
270  + a14 / w(i, j, k, irho) + a15 * (1 - gamma(i, j, k)) * velzrho
271  b15 = a11 * (gamma(i, j, k) - 1) + a15 * (gamma(i, j, k) - 1)
272 
273  b21 = a21 * (gamma(i, j, k) - 1) * q / 2 + a22 * (-velxrho) &
274  / w(i, j, k, irho) + a23 * (-velyrho) / w(i, j, k, irho) + a24 * (-velzrho) &
275  / w(i, j, k, irho) + a25 * (((gamma(i, j, k) - 1) * q / 2) - sos**2)
276  b22 = a21 * (1 - gamma(i, j, k)) * velxrho + a22 &
277  / w(i, j, k, irho) + a25 * (1 - gamma(i, j, k)) * velxrho
278  b23 = a21 * (1 - gamma(i, j, k)) * velyrho &
279  + a23 * 1 / w(i, j, k, irho) + a25 * (1 - gamma(i, j, k)) * velyrho
280  b24 = a21 * (1 - gamma(i, j, k)) * velzrho &
281  + a24 * 1 / w(i, j, k, irho) + a25 * (1 - gamma(i, j, k)) * velzrho
282  b25 = a21 * (gamma(i, j, k) - 1) + a25 * (gamma(i, j, k) - 1)
283 
284  b31 = a31 * (gamma(i, j, k) - 1) * q / 2 + a32 * (-velxrho) &
285  / w(i, j, k, irho) + a33 * (-velyrho) / w(i, j, k, irho) + &
286  a34 * (-velzrho) / w(i, j, k, irho) &
287  + a35 * (((gamma(i, j, k) - 1) * q / 2) - sos**2)
288  b32 = a31 * (1 - gamma(i, j, k)) * velxrho + a32 &
289  / w(i, j, k, irho) + a35 * (1 - gamma(i, j, k)) * velxrho
290  b33 = a31 * (1 - gamma(i, j, k)) * velyrho &
291  + a33 * 1 / w(i, j, k, irho) + a35 * (1 - gamma(i, j, k)) * velyrho
292  b34 = a31 * (1 - gamma(i, j, k)) * velzrho &
293  + a34 * 1 / w(i, j, k, irho) + a35 * (1 - gamma(i, j, k)) * velzrho
294  b35 = a31 * (gamma(i, j, k) - 1) + a35 * (gamma(i, j, k) - 1)
295 
296  b41 = a41 * (gamma(i, j, k) - 1) * q / 2 + a42 * (-velxrho) &
297  / w(i, j, k, irho) + a43 * (-velyrho) / w(i, j, k, irho) + a44 * (-velzrho) &
298  / w(i, j, k, irho) + a45 * (((gamma(i, j, k) - 1) * q / 2) - sos**2)
299  b42 = a41 * (1 - gamma(i, j, k)) * velxrho + a42 &
300  / w(i, j, k, irho) + a45 * (1 - gamma(i, j, k)) * velxrho
301  b43 = a41 * (1 - gamma(i, j, k)) * velyrho &
302  + a43 * 1 / w(i, j, k, irho) + a45 * (1 - gamma(i, j, k)) * velyrho
303  b44 = a41 * (1 - gamma(i, j, k)) * velzrho &
304  + a44 * 1 / w(i, j, k, irho) + a45 * (1 - gamma(i, j, k)) * velzrho
305  b45 = a41 * (gamma(i, j, k) - 1) + a45 * (gamma(i, j, k) - 1)
306 
307  b51 = a51 * (gamma(i, j, k) - 1) * q / 2 + a52 * (-velxrho) &
308  / w(i, j, k, irho) + a53 * (-velyrho) / w(i, j, k, irho) + a54 * (-velzrho) &
309  / w(i, j, k, irho) + a55 * (((gamma(i, j, k) - 1) * q / 2) - sos**2)
310  b52 = a51 * (1 - gamma(i, j, k)) * velxrho + a52 &
311  / w(i, j, k, irho) + a55 * (1 - gamma(i, j, k)) * velxrho
312  b53 = a51 * (1 - gamma(i, j, k)) * velyrho &
313  + a53 * 1 / w(i, j, k, irho) + a55 * (1 - gamma(i, j, k)) * velyrho
314  b54 = a51 * (1 - gamma(i, j, k)) * velzrho &
315  + a54 * 1 / w(i, j, k, irho) + a55 * (1 - gamma(i, j, k)) * velzrho
316  b55 = a51 * (gamma(i, j, k) - 1) + a55 * (gamma(i, j, k) - 1)
317 
318  ! dwo is the orginal redisual
319  do l = 1, nwf
320  dwo(l) = (dw(i, j, k, l) + fw(i, j, k, l)) * max(real(iblank(i, j, k), realtype), zero)
321  end do
322 
323  dw(i, j, k, 1) = b11 * dwo(1) + b12 * dwo(2) + b13 * dwo(3) + b14 * dwo(4) + b15 * dwo(5)
324  dw(i, j, k, 2) = b21 * dwo(1) + b22 * dwo(2) + b23 * dwo(3) + b24 * dwo(4) + b25 * dwo(5)
325  dw(i, j, k, 3) = b31 * dwo(1) + b32 * dwo(2) + b33 * dwo(3) + b34 * dwo(4) + b35 * dwo(5)
326  dw(i, j, k, 4) = b41 * dwo(1) + b42 * dwo(2) + b43 * dwo(3) + b44 * dwo(4) + b45 * dwo(5)
327  dw(i, j, k, 5) = b51 * dwo(1) + b52 * dwo(2) + b53 * dwo(3) + b54 * dwo(4) + b55 * dwo(5)
328 
329  end do
330  end do
331  end do ! end of lowspeedpreconditioners three cells loops
332 
333  else ! else.. i.e. if we do not have preconditioner turned on...
334  do l = 1, nwf
335  do k = 2, kl
336  do j = 2, jl
337  do i = 2, il
338  dw(i, j, k, l) = (dw(i, j, k, l) + fw(i, j, k, l)) &
339  * max(real(iblank(i, j, k), realtype), zero)
340  end do
341  end do
342  end do
343  end do
344  end if
345 
346  end subroutine residual_block
347 
348  subroutine sourceterms_block(nn, res, iRegion, pLocal)
349 
350  ! Apply the source terms for the given block. Assume that the
351  ! block pointers are already set.
352  use constants
353  use actuatorregiondata
354  use blockpointers, only: vol, dw, w
355  use flowvarrefstate, only: pref, uref, lref
356  use communication
357  use iteration, only: ordersconverged
358  implicit none
359 
360  ! Input
361  integer(kind=intType), intent(in) :: nn, iRegion
362  logical, intent(in) :: res
363  real(kind=realtype), intent(inout) :: plocal
364 
365  ! Working
366  integer(kind=intType) :: i, j, k, ii, iStart, iEnd
367  real(kind=realtype) :: ftmp(3), vx, vy, vz, f_fact(3), q_fact, qtmp, redim, factor, ostart, oend
368 
369  redim = pref * uref
370 
371  ! Compute the relaxation factor based on the ordersConverged
372 
373  ! How far we are into the ramp:
374  if (ordersconverged < actuatorregions(iregion)%relaxStart) then
375  factor = zero
376  else if (ordersconverged > actuatorregions(iregion)%relaxEnd) then
377  factor = one
378  else ! In between
379  ostart = actuatorregions(iregion)%relaxStart
380  oend = actuatorregions(iregion)%relaxEnd
381  factor = (ordersconverged - ostart) / (oend - ostart)
382  end if
383 
384  ! Compute the constant force factor
385  f_fact = factor * actuatorregions(iregion)%force / actuatorregions(iregion)%volume / pref
386 
387  ! Heat factor. This is heat added per unit volume per unit time
388  q_fact = factor * actuatorregions(iregion)%heat / actuatorregions(iregion)%volume / (pref * uref * lref * lref)
389 
390  ! Loop over the ranges for this block
391  istart = actuatorregions(iregion)%blkPtr(nn - 1) + 1
392  iend = actuatorregions(iregion)%blkPtr(nn)
393 
394  !$AD II-LOOP
395  do ii = istart, iend
396 
397  ! Extract the cell ID.
398  i = actuatorregions(iregion)%cellIDs(1, ii)
399  j = actuatorregions(iregion)%cellIDs(2, ii)
400  k = actuatorregions(iregion)%cellIDs(3, ii)
401 
402  ! This actually gets the force
403  ftmp = vol(i, j, k) * f_fact
404 
405  vx = w(i, j, k, ivx)
406  vy = w(i, j, k, ivy)
407  vz = w(i, j, k, ivz)
408 
409  ! this gets the heat addition rate
410  qtmp = vol(i, j, k) * q_fact
411 
412  if (res) then
413  ! Momentum residuals
414  dw(i, j, k, imx:imz) = dw(i, j, k, imx:imz) - ftmp
415 
416  ! energy residuals
417  dw(i, j, k, irhoe) = dw(i, j, k, irhoe) - &
418  ftmp(1) * vx - ftmp(2) * vy - ftmp(3) * vz - qtmp
419  else
420  ! Add in the local power contribution:
421  plocal = plocal + (vx * ftmp(1) + vy * ftmp(2) + vz * ftmp(3)) * redim
422  end if
423  end do
424 
425  end subroutine sourceterms_block
426 
427  subroutine initres_block(varStart, varEnd, nn, sps)
428  !
429  ! initres initializes the given range of the residual. Either to
430  ! zero, steady computation, or to an unsteady term for the time
431  ! spectral and unsteady modes. For the coarser grid levels the
432  ! residual forcing term is taken into account.
433  !
434  use blockpointers
435  use flowvarrefstate
436  use inputiteration
437  use inputphysics
438  use inputtimespectral
439  use inputunsteady
440  use iteration
441 
442  implicit none
443  !
444  ! Subroutine arguments.
445  !
446  integer(kind=intType), intent(in) :: varStart, varEnd, nn, sps
447  !
448  ! Local variables.
449  !
450  integer(kind=intType) :: mm, ll, ii, jj, i, j, k, l, m
451  real(kind=realtype) :: oneoverdt, tmp
452 
453  real(kind=realtype), dimension(:, :, :, :), pointer :: ww, wsp, wsp1
454  real(kind=realtype), dimension(:, :, :), pointer :: volsp
455 
456  ! Return immediately of no variables are in the range.
457 
458  if (varend < varstart) return
459 
460  ! Determine the equation mode and act accordingly.
461 
462  select case (equationmode)
463  case (steady)
464 
465  ! Steady state computation.
466  ! Determine the currently active multigrid level.
467 
468  steadyleveltest: if (currentlevel == groundlevel) then
469 
470  ! Ground level of the multigrid cycle. Initialize the
471  ! owned residuals to zero.
472 
473  do l = varstart, varend
474  do k = 2, kl
475  do j = 2, jl
476  do i = 2, il
477  dw(i, j, k, l) = zero
478  end do
479  end do
480  end do
481  end do
482 
483  else steadyleveltest
484 
485  ! Coarse grid level. Initialize the owned cells to the
486  ! residual forcing terms.
487 
488  do l = varstart, varend
489  do k = 2, kl
490  do j = 2, jl
491  do i = 2, il
492  dw(i, j, k, l) = wr(i, j, k, l)
493  end do
494  end do
495  end do
496  end do
497  end if steadyleveltest
498 
499 #ifndef USE_TAPENADE
500 
501  case (unsteady)
502 
503  ! Unsteady computation.
504  ! A further distinction must be made.
505 
506  select case (timeintegrationscheme)
507 
508  case (explicitrk)
509 
510  ! We are always on the finest grid.
511  ! Initialize the residual to zero.
512 
513  do l = varstart, varend
514  do k = 2, kl
515  do j = 2, jl
516  do i = 2, il
517  dw(i, j, k, l) = zero
518  end do
519  end do
520  end do
521  end do
522 
523  !=======================================================
524 
525  case (md, bdf) ! Modified by HDN
526 
527  ! Store the inverse of the physical nonDimensional
528  ! time step a bit easier.
529 
530  oneoverdt = timeref / deltat
531 
532  ! Store the pointer for the variable to be used to compute
533  ! the unsteady source term. For a runge-kutta smoother this
534  ! is the solution of the zeroth runge-kutta stage. As for
535  ! rkStage == 0 this variable is not yet set w is used.
536  ! For other smoothers w is to be used as well.
537 
538  if (smoother == rungekutta .and. rkstage > 0) then
539  ww => wn
540  else
541  ww => w
542  end if
543 
544  ! Determine the currently active multigrid level.
545 
546  unsteadyleveltest: if (currentlevel == groundlevel) then
547 
548  ! Ground level of the multigrid cycle. Initialize the
549  ! owned cells to the unsteady source term. First the
550  ! term for the current time level. Note that in w the
551  ! velocities are stored and not the momentum variables.
552  ! Therefore the if-statement is present to correct this.
553 
554  do l = varstart, varend
555 
556  if (l == ivx .or. l == ivy .or. l == ivz) then
557 
558  ! Momentum variables.
559 
560  do k = 2, kl
561  do j = 2, jl
562  do i = 2, il
563  dw(i, j, k, l) = coeftime(0) * vol(i, j, k) &
564  * ww(i, j, k, l) * ww(i, j, k, irho)
565  end do
566  end do
567  end do
568 
569  else
570 
571  ! Non-momentum variables, for which the variable
572  ! to be solved is stored; for the flow equations this
573  ! is the conservative variable, for the turbulent
574  ! equations the primitive variable.
575 
576  do k = 2, kl
577  do j = 2, jl
578  do i = 2, il
579  dw(i, j, k, l) = coeftime(0) * vol(i, j, k) &
580  * ww(i, j, k, l)
581  end do
582  end do
583  end do
584 
585  end if
586 
587  end do
588 
589  ! The terms from the older time levels. Here the
590  ! conservative variables are stored. In case of a
591  ! deforming mesh, also the old volumes must be taken.
592 
593  deformingtest: if (deforming_grid) then
594 
595  ! Mesh is deforming and thus the volumes can change.
596  ! Use the old volumes as well.
597 
598  do m = 1, noldlevels
599  do l = varstart, varend
600  do k = 2, kl
601  do j = 2, jl
602  do i = 2, il
603  dw(i, j, k, l) = dw(i, j, k, l) &
604  + coeftime(m) * volold(m, i, j, k) &
605  * wold(m, i, j, k, l)
606  end do
607  end do
608  end do
609  end do
610  end do
611 
612  else deformingtest
613 
614  ! Rigid mesh. The volumes remain constant.
615 
616  do m = 1, noldlevels
617  do l = varstart, varend
618  do k = 2, kl
619  do j = 2, jl
620  do i = 2, il
621  dw(i, j, k, l) = dw(i, j, k, l) &
622  + coeftime(m) * vol(i, j, k) &
623  * wold(m, i, j, k, l)
624  end do
625  end do
626  end do
627  end do
628  end do
629 
630  end if deformingtest
631 
632  ! Multiply the time derivative by the inverse of the
633  ! time step to obtain the true time derivative.
634  ! This is done after the summation has been done, because
635  ! otherwise you run into finite accuracy problems for
636  ! very small time steps.
637 
638  do l = varstart, varend
639  do k = 2, kl
640  do j = 2, jl
641  do i = 2, il
642  dw(i, j, k, l) = oneoverdt * dw(i, j, k, l)
643  end do
644  end do
645  end do
646  end do
647 
648  else unsteadyleveltest
649 
650  ! Coarse grid level. Initialize the owned cells to the
651  ! residual forcing term plus a correction for the
652  ! multigrid treatment of the time derivative term.
653  ! As the velocities are stored instead of the momentum,
654  ! these terms must be multiplied by the density.
655 
656  tmp = oneoverdt * coeftime(0)
657 
658  do l = varstart, varend
659 
660  if (l == ivx .or. l == ivy .or. l == ivz) then
661 
662  ! Momentum variables.
663 
664  do k = 2, kl
665  do j = 2, jl
666  do i = 2, il
667  dw(i, j, k, l) = tmp * vol(i, j, k) &
668  * (ww(i, j, k, l) * ww(i, j, k, irho) &
669  - w1(i, j, k, l) * w1(i, j, k, irho))
670  dw(i, j, k, l) = dw(i, j, k, l) + wr(i, j, k, l)
671  end do
672  end do
673  end do
674 
675  else
676 
677  ! Non-momentum variables.
678 
679  do k = 2, kl
680  do j = 2, jl
681  do i = 2, il
682  dw(i, j, k, l) = tmp * vol(i, j, k) &
683  * (ww(i, j, k, l) - w1(i, j, k, l))
684  dw(i, j, k, l) = dw(i, j, k, l) + wr(i, j, k, l)
685  end do
686  end do
687  end do
688 
689  end if
690 
691  end do
692 
693  end if unsteadyleveltest
694 
695  end select
696 
697  !===========================================================
698 
699  case (timespectral)
700 
701  ! Time spectral computation. The time derivative of the
702  ! current solution is given by a linear combination of
703  ! all other solutions, i.e. a matrix vector product.
704 
705  ! First store the section to which this block belongs
706  ! in jj.
707 
708  jj = sectionid
709 
710  ! Determine the currently active multigrid level.
711 
712  spectralleveltest: if (currentlevel == groundlevel) then
713 
714  ! Finest multigrid level. The residual must be
715  ! initialized to the time derivative.
716 
717  ! Initialize it to zero.
718 
719  do l = varstart, varend
720  do k = 2, kl
721  do j = 2, jl
722  do i = 2, il
723  dw(i, j, k, l) = zero
724  end do
725  end do
726  end do
727  end do
728 
729  ! Loop over the number of terms which contribute
730  ! to the time derivative.
731 
732  timeloopfine: do mm = 1, ntimeintervalsspectral
733 
734  ! Store the pointer for the variable to be used to
735  ! compute the unsteady source term and the volume.
736  ! Also store in ii the offset needed for vector
737  ! quantities.
738 
739  wsp => flowdoms(nn, currentlevel, mm)%w
740  volsp => flowdoms(nn, currentlevel, mm)%vol
741  ii = 3 * (mm - 1)
742 
743  ! Loop over the number of variables to be set.
744 
745  varloopfine: do l = varstart, varend
746 
747  ! Test for a momentum variable.
748 
749  if (l == ivx .or. l == ivy .or. l == ivz) then
750 
751  ! Momentum variable. A special treatment is
752  ! needed because it is a vector and the velocities
753  ! are stored instead of the momentum. Set the
754  ! coefficient ll, which defines the row of the
755  ! matrix used later on.
756 
757  if (l == ivx) ll = 3 * sps - 2
758  if (l == ivy) ll = 3 * sps - 1
759  if (l == ivz) ll = 3 * sps
760 
761  ! Loop over the owned cell centers to add the
762  ! contribution from wsp.
763 
764  do k = 2, kl
765  do j = 2, jl
766  do i = 2, il
767 
768  ! Store the matrix vector product with the
769  ! velocity in tmp.
770 
771  tmp = dvector(jj, ll, ii + 1) * wsp(i, j, k, ivx) &
772  + dvector(jj, ll, ii + 2) * wsp(i, j, k, ivy) &
773  + dvector(jj, ll, ii + 3) * wsp(i, j, k, ivz)
774 
775  ! Update the residual. Note the
776  ! multiplication with the density to obtain
777  ! the correct time derivative for the
778  ! momentum variable.
779 
780  dw(i, j, k, l) = dw(i, j, k, l) &
781  + tmp * volsp(i, j, k) * wsp(i, j, k, irho)
782 
783  end do
784  end do
785  end do
786 
787  else
788 
789  ! Scalar variable. Loop over the owned cells to
790  ! add the contribution of wsp to the time
791  ! derivative.
792 
793  do k = 2, kl
794  do j = 2, jl
795  do i = 2, il
796  dw(i, j, k, l) = dw(i, j, k, l) &
797  + dscalar(jj, sps, mm) &
798  * volsp(i, j, k) * wsp(i, j, k, l)
799 
800  end do
801  end do
802  end do
803 
804  end if
805 
806  end do varloopfine
807 
808  end do timeloopfine
809  else spectralleveltest
810 
811  ! Coarse grid level. Initialize the owned cells to the
812  ! residual forcing term plus a correction for the
813  ! multigrid treatment of the time derivative term.
814 
815  ! Initialization to the residual forcing term.
816 
817  do l = varstart, varend
818  do k = 2, kl
819  do j = 2, jl
820  do i = 2, il
821  dw(i, j, k, l) = wr(i, j, k, l)
822  end do
823  end do
824  end do
825  end do
826 
827  ! Loop over the number of terms which contribute
828  ! to the time derivative.
829 
830  timeloopcoarse: do mm = 1, ntimeintervalsspectral
831 
832  ! Store the pointer for the variable to be used to
833  ! compute the unsteady source term and the pointers
834  ! for wsp1, the solution when entering this MG level
835  ! and for the volume.
836  ! Furthermore store in ii the offset needed for
837  ! vector quantities.
838 
839  wsp => flowdoms(nn, currentlevel, mm)%w
840  wsp1 => flowdoms(nn, currentlevel, mm)%w1
841  volsp => flowdoms(nn, currentlevel, mm)%vol
842  ii = 3 * (mm - 1)
843 
844  ! Loop over the number of variables to be set.
845 
846  varloopcoarse: do l = varstart, varend
847 
848  ! Test for a momentum variable.
849 
850  if (l == ivx .or. l == ivy .or. l == ivz) then
851 
852  ! Momentum variable. A special treatment is
853  ! needed because it is a vector and the velocities
854  ! are stored instead of the momentum. Set the
855  ! coefficient ll, which defines the row of the
856  ! matrix used later on.
857 
858  if (l == ivx) ll = 3 * sps - 2
859  if (l == ivy) ll = 3 * sps - 1
860  if (l == ivz) ll = 3 * sps
861 
862  ! Add the contribution of wps to the correction
863  ! of the time derivative. The difference between
864  ! the current time derivative and the one when
865  ! entering this grid level must be added, because
866  ! the residual forcing term only takes the spatial
867  ! part of the coarse grid residual into account.
868 
869  do k = 2, kl
870  do j = 2, jl
871  do i = 2, il
872 
873  ! Store the matrix vector product with the
874  ! momentum in tmp.
875 
876  tmp = dvector(jj, ll, ii + 1) &
877  * (wsp(i, j, k, irho) * wsp(i, j, k, ivx) &
878  - wsp1(i, j, k, irho) * wsp1(i, j, k, ivx)) &
879  + dvector(jj, ll, ii + 2) &
880  * (wsp(i, j, k, irho) * wsp(i, j, k, ivy) &
881  - wsp1(i, j, k, irho) * wsp1(i, j, k, ivy)) &
882  + dvector(jj, ll, ii + 3) &
883  * (wsp(i, j, k, irho) * wsp(i, j, k, ivz) &
884  - wsp1(i, j, k, irho) * wsp1(i, j, k, ivz))
885 
886  ! Add tmp to the residual. Multiply it by
887  ! the volume to obtain the finite volume
888  ! formulation of the derivative of the
889  ! momentum.
890 
891  dw(i, j, k, l) = dw(i, j, k, l) + tmp * volsp(i, j, k)
892 
893  end do
894  end do
895  end do
896 
897  else
898 
899  ! Scalar variable. Loop over the owned cells
900  ! to add the contribution of wsp to the correction
901  ! of the time derivative.
902 
903  do k = 2, kl
904  do j = 2, jl
905  do i = 2, il
906  dw(i, j, k, l) = dw(i, j, k, l) &
907  + dscalar(jj, sps, mm) &
908  * volsp(i, j, k) &
909  * (wsp(i, j, k, l) - wsp1(i, j, k, l))
910  end do
911  end do
912  end do
913 
914  end if
915 
916  end do varloopcoarse
917 
918  end do timeloopcoarse
919  end if spectralleveltest
920 #endif
921  end select
922 
923  ! Set the residual in the halo cells to zero. This is just
924  ! to avoid possible problems. Their values do not matter.
925 
926  do l = varstart, varend
927  do k = 0, kb
928  do j = 0, jb
929  dw(0, j, k, l) = zero
930  dw(1, j, k, l) = zero
931  dw(ie, j, k, l) = zero
932  dw(ib, j, k, l) = zero
933  end do
934  end do
935 
936  do k = 0, kb
937  do i = 2, il
938  dw(i, 0, k, l) = zero
939  dw(i, 1, k, l) = zero
940  dw(i, je, k, l) = zero
941  dw(i, jb, k, l) = zero
942  end do
943  end do
944 
945  do j = 2, jl
946  do i = 2, il
947  dw(i, j, 0, l) = zero
948  dw(i, j, 1, l) = zero
949  dw(i, j, ke, l) = zero
950  dw(i, j, kb, l) = zero
951  end do
952  end do
953  end do
954 
955  end subroutine initres_block
956 
957  ! ----------------------------------------------------------------------
958  ! |
959  ! No Tapenade Routine below this line |
960  ! |
961  ! ----------------------------------------------------------------------
962 
963 #ifndef USE_TAPENADE
964  subroutine initres(varStart, varEnd)
965  !
966  ! Shell function to call initRes_block on all blocks
967  !
968  use blockpointers
969  use constants
970  use inputtimespectral
971  use iteration
972  use section
973  use utils, only: setpointers
974  !
975  ! Subroutine argument.
976  !
977  integer(kind=intType), intent(in) :: varStart, varEnd
978  !
979  ! Local variables.
980  !
981  integer(kind=intType) :: sps, nn
982 
983  ! Loop over the number of spectral solutions.
984 
985  spectralloop: do sps = 1, ntimeintervalsspectral
986 
987  ! Loop over the number of blocks.
988 
989  domains: do nn = 1, ndom
990 
991  ! Set the pointers for this block.
992 
993  call setpointers(nn, currentlevel, sps)
994 
995  call initres_block(varstart, varend, nn, sps)
996 
997  end do domains
998 
999  end do spectralloop
1000 
1001  end subroutine initres
1002 
1003  subroutine sourceterms
1004 
1005  ! Shell function to call sourceTerms_block on all blocks
1006 
1007  use constants
1008  use blockpointers
1009  use utils, only: setpointers
1010  use actuatorregiondata
1011  implicit none
1012 
1013  integer(kind=intType) :: nn, iRegion
1014  real(kind=realtype) :: dummy
1015 
1016  ! Loop over the number of domains.
1017  domains: do nn = 1, ndom
1018 
1019  ! Set the pointers for this block.
1020  call setpointers(nn, 1, 1)
1021  do iregion = 1, nactuatorregions
1022  call sourceterms_block(nn, .true., iregion, dummy)
1023  end do
1024  end do domains
1025 
1026  end subroutine sourceterms
1027 
1028  subroutine residual
1029  !
1030  ! Shell function to call residual_block on all blocks
1031  !
1032  use blockpointers
1033  use constants
1034  use inputtimespectral
1035  use iteration
1036  use utils, only: setpointers
1037  implicit none
1038  !
1039  ! Local variables.
1040  !
1041  integer(kind=intType) :: sps, nn
1042 
1043  ! Loop over the number of spectral solutions.
1044 
1045  spectralloop: do sps = 1, ntimeintervalsspectral
1046 
1047  ! Loop over the number of blocks.
1048 
1049  domains: do nn = 1, ndom
1050 
1051  ! Set the pointers for this block.
1052 
1053  call setpointers(nn, currentlevel, sps)
1054 
1055  call residual_block
1056 
1057  end do domains
1058 
1059  end do spectralloop
1060  end subroutine residual
1061 
1062  subroutine computedwdadi
1063  !
1064  ! executeRkStage executes one runge kutta stage. The stage
1065  ! number, rkStage, is defined in the local module iteration.
1066  !
1067  use blockpointers
1068  use constants
1069  use flowvarrefstate
1070  use inputiteration
1071  use inputphysics
1072  use inputtimespectral
1073  use inputunsteady
1074  use iteration
1075  implicit none
1076  !
1077  ! Local parameter.
1078  !
1079  ! Local variables.
1080  !
1081  integer(kind=intType) :: i, j, k, n
1082  real(kind=realtype) :: gm1, epsval, fac, eps2
1083  real(kind=realtype) :: uvel, vvel, wvel, cijk, cijkinv, c2inv, uvw
1084  real(kind=realtype) :: ri1, ri2, ri3, rj1, rj2, rj3, rk1, rk2, rk3, uu
1085  real(kind=realtype) :: ri, rj, rk, qsi, qsj, qsk, currentcfl
1086  real(kind=realtype) :: xfact, cinf, cinf2
1087  real(kind=realtype) :: dw1, dw2, dw3, dw4, dw5
1088  real(kind=realtype) :: a1, a2, a3, a4, a5, a6, a7, mut, ge
1089  real(kind=realtype) :: viscterm1, viscterm2, viscterm3
1090  real(kind=realtype) :: metterm
1091  real(kind=realtype) :: unsteadyimpl, mult
1092  real(kind=realtype) :: sqrt2, sqrt2inv, alph, alphinv
1093  real(kind=realtype) :: volhalf, volhalfrho, volfact
1094  real(kind=realtype) :: mutpi, mutmi, mutpj, mutmj, mutpk, mutmk, &
1095  mettermp, mettermm, &
1096  visctermi, visctermj, visctermk
1097 
1098  real(kind=realtype), dimension(:, :, :), pointer :: qq_i, qq_j, qq_k, &
1099  cc_i, cc_j, cc_k, spectral_i, &
1100  spectral_j, spectral_k, dual_dt
1101  real(kind=realtype), dimension(ie, 5) :: bbi, cci, ddi, ffi
1102  real(kind=realtype), dimension(je, 5) :: bbj, ccj, ddj, ffj
1103  real(kind=realtype), dimension(ke, 5) :: bbk, cck, ddk, ffk
1104  real(kind=realtype), dimension(ie) :: mettermi
1105  real(kind=realtype), dimension(je) :: mettermj
1106  real(kind=realtype), dimension(ke) :: mettermk
1107  real(kind=realtype), dimension(5) :: diagplus, diagminus
1108 
1109  ! Set the value of the current cfl number,
1110  ! depending on the situation. On the finest grid in the mg cycle
1111  ! the second halo is computed, otherwise not.
1112 
1113  currentcfl = cflcoarse
1114  if (currentlevel == 1) then
1115  currentcfl = cfl
1116  end if
1117  qq_i => scratch(:, :, :, 1)
1118  qq_j => scratch(:, :, :, 2)
1119  qq_k => scratch(:, :, :, 3)
1120  cc_i => scratch(:, :, :, 4)
1121  cc_j => scratch(:, :, :, 5)
1122  cc_k => scratch(:, :, :, 6)
1123  spectral_i => scratch(:, :, :, 7)
1124  spectral_j => scratch(:, :, :, 8)
1125  spectral_k => scratch(:, :, :, 9)
1126  dual_dt => scratch(:, :, :, 10)
1127 
1128  ! havent thought about iblank
1129 
1130  ! these are factors for robustness.
1131 
1132  epsval = 0.08
1133  fac = 1.05
1134  mut = zero
1135  cinf2 = gammainf * pinf / rhoinf
1136  cinf = sqrt(cinf2)
1137 
1138  qsi = zero
1139  qsj = zero
1140  qsk = zero
1141  visctermi = zero
1142  visctermj = zero
1143  visctermk = zero
1144  sqrt2 = sqrt(two)
1145  sqrt2inv = one / sqrt2
1146 
1147  if (equationmode .eq. steady) then
1148  do k = 2, kl
1149  do j = 2, jl
1150  do i = 2, il
1151  dual_dt(i, j, k) = currentcfl * dtl(i, j, k) * vol(i, j, k)
1152  end do
1153  end do
1154  end do
1155  else
1156  do k = 2, kl
1157  do j = 2, jl
1158  do i = 2, il
1159  unsteadyimpl = coeftime(0) * timeref / deltat
1160  mult = currentcfl * dtl(i, j, k)
1161  mult = mult / (mult * unsteadyimpl * vol(i, j, k) + one)
1162  dual_dt(i, j, k) = mult * vol(i, j, k)
1163  end do
1164  end do
1165  end do
1166  end if
1167 
1168  ! Set up some arrays
1169  do k = 2, kl
1170  do j = 2, jl
1171  do i = 2, il
1172  volhalf = half / vol(i, j, k)
1173  volhalfrho = volhalf / w(i, j, k, irho)
1174  uvel = w(i, j, k, ivx)
1175  vvel = w(i, j, k, ivy)
1176  wvel = w(i, j, k, ivz)
1177 
1178  cijk = sqrt(gamma(i, j, k) * p(i, j, k) / w(i, j, k, irho))
1179 
1180  ri1 = volhalf * (si(i, j, k, 1) + si(i - 1, j, k, 1))
1181  ri2 = volhalf * (si(i, j, k, 2) + si(i - 1, j, k, 2))
1182  ri3 = volhalf * (si(i, j, k, 3) + si(i - 1, j, k, 3))
1183 
1184  rj1 = volhalf * (sj(i, j, k, 1) + sj(i, j - 1, k, 1))
1185  rj2 = volhalf * (sj(i, j, k, 2) + sj(i, j - 1, k, 2))
1186  rj3 = volhalf * (sj(i, j, k, 3) + sj(i, j - 1, k, 3))
1187 
1188  rk1 = volhalf * (sk(i, j, k, 1) + sk(i, j, k - 1, 1))
1189  rk2 = volhalf * (sk(i, j, k, 2) + sk(i, j, k - 1, 2))
1190  rk3 = volhalf * (sk(i, j, k, 3) + sk(i, j, k - 1, 3))
1191 
1192  if (addgridvelocities) then
1193  qsi = (sfacei(i - 1, j, k) + sfacei(i, j, k)) * volhalf
1194  qsj = (sfacej(i, j - 1, k) + sfacej(i, j, k)) * volhalf
1195  qsk = (sfacek(i, j, k - 1) + sfacek(i, j, k)) * volhalf
1196  end if
1197 
1198  qq_i(i, j, k) = ri1 * w(i, j, k, ivx) + ri2 * w(i, j, k, ivy) + &
1199  ri3 * w(i, j, k, ivz) - qsi
1200  qq_j(i, j, k) = rj1 * w(i, j, k, ivx) + rj2 * w(i, j, k, ivy) + &
1201  rj3 * w(i, j, k, ivz) - qsj
1202  qq_k(i, j, k) = rk1 * w(i, j, k, ivx) + rk2 * w(i, j, k, ivy) + &
1203  rk3 * w(i, j, k, ivz) - qsk
1204 
1205  ri = sqrt(ri1 * ri1 + ri2 * ri2 + ri3 * ri3)
1206  rj = sqrt(rj1 * rj1 + rj2 * rj2 + rj3 * rj3)
1207  rk = sqrt(rk1 * rk1 + rk2 * rk2 + rk3 * rk3)
1208 
1209  cc_i(i, j, k) = cijk * ri
1210  cc_j(i, j, k) = cijk * rj
1211  cc_k(i, j, k) = cijk * rk
1212  if (viscous) then
1213  mutpi = rlv(i, j, k) + rlv(i + 1, j, k)
1214  mutmi = rlv(i, j, k) + rlv(i - 1, j, k)
1215  mutpj = rlv(i, j, k) + rlv(i, j + 1, k)
1216  mutmj = rlv(i, j, k) + rlv(i, j - 1, k)
1217  mutpk = rlv(i, j, k) + rlv(i, j, k + 1)
1218  mutmk = rlv(i, j, k) + rlv(i, j, k - 1)
1219  if (eddymodel) then
1220  mutpi = mutpi + rev(i, j, k) + rev(i + 1, j, k)
1221  mutmi = mutpi + rev(i, j, k) + rev(i - 1, j, k)
1222  mutpj = mutpj + rev(i, j, k) + rev(i, j + 1, k)
1223  mutmj = mutpj + rev(i, j, k) + rev(i, j - 1, k)
1224  mutpk = mutpk + rev(i, j, k) + rev(i, j, k + 1)
1225  mutmk = mutpk + rev(i, j, k) + rev(i, j, k - 1)
1226  end if
1227 
1228  volfact = two / (vol(i, j, k) + vol(i + 1, j, k))
1229  mettermp = (si(i, j, k, 1) * si(i, j, k, 1) &
1230  + si(i, j, k, 2) * si(i, j, k, 2) &
1231  + si(i, j, k, 3) * si(i, j, k, 3)) * mutpi * volfact
1232  volfact = two / (vol(i, j, k) + vol(i - 1, j, k))
1233  mettermm = (si(i - 1, j, k, 1) * si(i - 1, j, k, 1) &
1234  + si(i - 1, j, k, 2) * si(i - 1, j, k, 2) &
1235  + si(i - 1, j, k, 3) * si(i - 1, j, k, 3)) * mutmi * volfact
1236  visctermi = (mettermp + mettermm) * volhalfrho
1237 
1238  volfact = two / (vol(i, j, k) + vol(i, j + 1, k))
1239  mettermp = (sj(i, j, k, 1) * sj(i, j, k, 1) &
1240  + sj(i, j, k, 2) * sj(i, j, k, 2) &
1241  + sj(i, j, k, 3) * sj(i, j, k, 3)) * mutpj * volfact
1242  volfact = two / (vol(i, j, k) + vol(i, j - 1, k))
1243  mettermm = (sj(i, j - 1, k, 1) * sj(i, j - 1, k, 1) &
1244  + sj(i, j - 1, k, 2) * sj(i, j - 1, k, 2) &
1245  + sj(i, j - 1, k, 3) * sj(i, j - 1, k, 3)) * mutmj * volfact
1246  visctermj = (mettermp + mettermm) * volhalfrho
1247 
1248  volfact = two / (vol(i, j, k) + vol(i, j, k + 1))
1249  mettermp = (sk(i, j, k, 1) * sk(i, j, k, 1) &
1250  + sk(i, j, k, 2) * sk(i, j, k, 2) &
1251  + sk(i, j, k, 3) * sk(i, j, k, 3)) * mutpk * volfact
1252  volfact = two / (vol(i, j, k) + vol(i, j, k - 1))
1253  mettermm = (sk(i, j, k - 1, 1) * sk(i, j, k - 1, 1) &
1254  + sk(i, j, k - 1, 2) * sk(i, j, k - 1, 2) &
1255  + sk(i, j, k - 1, 3) * sk(i, j, k - 1, 3)) * mutmk * volfact
1256  visctermk = (mettermp + mettermm) * volhalfrho
1257 
1258  end if
1259 
1260  eps2 = ri * cinf * epsval
1261  spectral_i(i, j, k) = (fac * sqrt((abs(qq_i(i, j, k)) + cc_i(i, j, k))**2 &
1262  + eps2**2) + visctermi) * dual_dt(i, j, k)
1263  eps2 = rj * cinf * epsval
1264  spectral_j(i, j, k) = (fac * sqrt((abs(qq_j(i, j, k)) + cc_j(i, j, k))**2 &
1265  + eps2**2) + visctermj) * dual_dt(i, j, k)
1266  eps2 = rk * cinf * epsval
1267  spectral_k(i, j, k) = (fac * sqrt((abs(qq_k(i, j, k)) + cc_k(i, j, k))**2 &
1268  + eps2**2) + visctermk) * dual_dt(i, j, k)
1269  spectral_i(i, j, k) = spectral_i(i, j, k) * zero
1270  spectral_j(i, j, k) = spectral_j(i, j, k) * zero
1271  spectral_k(i, j, k) = spectral_k(i, j, k) * zero
1272  end do
1273  end do
1274  end do
1275 
1276  ! Multiply by T_eta^inv
1277  do k = 2, kl
1278  do j = 2, jl
1279  do i = 2, il
1280 
1281  gm1 = gamma(i, j, k) - one
1282  cijk = sqrt(gamma(i, j, k) * p(i, j, k) / w(i, j, k, irho))
1283  c2inv = one / (cijk * cijk)
1284  xfact = two * cijk
1285  alphinv = sqrt2 * cijk / w(i, j, k, irho)
1286 
1287  uvel = w(i, j, k, ivx)
1288  vvel = w(i, j, k, ivy)
1289  wvel = w(i, j, k, ivz)
1290  uvw = half * (uvel * uvel + vvel * vvel + wvel * wvel)
1291 
1292  rj1 = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1))
1293  rj2 = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2))
1294  rj3 = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3))
1295  rj = sqrt(rj1 * rj1 + rj2 * rj2 + rj3 * rj3)
1296  uu = uvel * rj1 + vvel * rj2 + wvel * rj3
1297  rj1 = rj1 / rj
1298  rj2 = rj2 / rj
1299  rj3 = rj3 / rj
1300 
1301  dw1 = dw(i, j, k, 1)
1302  dw2 = dw(i, j, k, 2)
1303  dw3 = dw(i, j, k, 3)
1304  dw4 = dw(i, j, k, 4)
1305  dw5 = dw(i, j, k, 5)
1306 
1307  a1 = dw2 * uvel + dw3 * vvel + dw4 * wvel - dw5
1308  a1 = a1 * gm1 * c2inv + dw1 * (one - uvw * gm1 * c2inv)
1309 
1310  a2 = (rj2 * wvel - rj3 * vvel) * dw1 + rj3 * dw3 - rj2 * dw4
1311  a3 = (rj3 * uvel - rj1 * wvel) * dw1 + rj1 * dw4 - rj3 * dw2
1312  a4 = (rj1 * vvel - rj2 * uvel) * dw1 + rj2 * dw2 - rj1 * dw3
1313 
1314  a5 = uvw * dw1 - uvel * dw2 - vvel * dw3 - wvel * dw4 + dw5
1315  a5 = a5 * gm1 * c2inv
1316 
1317  a6 = uu * dw1 / rj - rj1 * dw2 - rj2 * dw3 - rj3 * dw4
1318 
1319  dw(i, j, k, 1) = a1 * rj1 + a2 / w(i, j, k, irho)
1320  dw(i, j, k, 2) = a1 * rj2 + a3 / w(i, j, k, irho)
1321  dw(i, j, k, 3) = a1 * rj3 + a4 / w(i, j, k, irho)
1322  dw(i, j, k, 4) = (half * a5 - a6 / xfact) * alphinv
1323  dw(i, j, k, 5) = (half * a5 + a6 / xfact) * alphinv
1324 
1325  end do
1326  end do
1327  end do
1328 
1329  if (jl .gt. 2) then
1330 
1331  ! Inversion in j
1332 
1333  do k = 2, kl
1334  do i = 2, il
1335  do j = 1, jl
1336  if (viscous) mut = rlv(i, j, k) + rlv(i, j + 1, k)
1337  if (eddymodel) mut = mut + rev(i, j, k) + rev(i, j + 1, k)
1338 
1339  volfact = one / (vol(i, j, k) + vol(i, j + 1, k))
1340  metterm = (sj(i, j, k, 1) * sj(i, j, k, 1) &
1341  + sj(i, j, k, 2) * sj(i, j, k, 2) &
1342  + sj(i, j, k, 3) * sj(i, j, k, 3))
1343  mettermj(j) = metterm * mut * volfact
1344  end do
1345 
1346  do j = 2, jl
1347 
1348  viscterm1 = mettermj(j) / vol(i, j, k) / w(i, j, k, irho)
1349  viscterm3 = mettermj(j - 1) / vol(i, j, k) / w(i, j, k, irho)
1350  viscterm2 = viscterm1 + viscterm3
1351 
1352  volhalf = half / vol(i, j, k)
1353  rj1 = volhalf * (sj(i, j, k, 1) + sj(i, j - 1, k, 1))
1354  rj2 = volhalf * (sj(i, j, k, 2) + sj(i, j - 1, k, 2))
1355  rj3 = volhalf * (sj(i, j, k, 3) + sj(i, j - 1, k, 3))
1356  metterm = rj1 * rj1 + rj2 * rj2 + rj3 * rj3
1357  eps2 = epsval * epsval * cinf2 * metterm
1358  diagplus(1) = half * (qq_j(i, j, k) + fac * sqrt(qq_j(i, j, k)**2 + eps2))
1359  diagplus(2) = diagplus(1)
1360  diagplus(3) = diagplus(1)
1361  diagplus(4) = half * (qq_j(i, j, k) + cc_j(i, j, k) &
1362  + fac * sqrt((qq_j(i, j, k) + cc_j(i, j, k))**2 + eps2))
1363  diagplus(5) = half * (qq_j(i, j, k) - cc_j(i, j, k) &
1364  + fac * sqrt((qq_j(i, j, k) - cc_j(i, j, k))**2 + eps2))
1365  diagminus(1) = half * (qq_j(i, j, k) - fac * sqrt(qq_j(i, j, k)**2 + eps2))
1366  diagminus(2) = diagminus(1)
1367  diagminus(3) = diagminus(1)
1368  diagminus(4) = half * (qq_j(i, j, k) + cc_j(i, j, k) &
1369  - fac * sqrt((qq_j(i, j, k) + cc_j(i, j, k))**2 + eps2))
1370  diagminus(5) = half * (qq_j(i, j, k) - cc_j(i, j, k) &
1371  - fac * sqrt((qq_j(i, j, k) - cc_j(i, j, k))**2 + eps2))
1372 
1373  do n = 1, 5
1374  bbj(j + 1, n) = -viscterm1 - diagplus(n)
1375  ddj(j - 1, n) = -viscterm3 + diagminus(n)
1376  ccj(j, n) = viscterm2 + diagplus(n) - diagminus(n)
1377  end do
1378  end do
1379 
1380  do n = 1, 5
1381  bbj(je, n) = zero
1382  ddj(1, n) = zero
1383  do j = 2, jl
1384  bbj(j, n) = bbj(j, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1385  ddj(j, n) = ddj(j, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1386  ccj(j, n) = one + ccj(j, n) * dual_dt(i, j, k) * &
1387  max(real(iblank(i, j, k), realtype), zero) + &
1388  spectral_i(i, j, k) + spectral_k(i, j, k)
1389  ffj(j, n) = dw(i, j, k, n)
1390  end do
1391  end do
1392 
1393  call tridiagsolve(bbj, ccj, ddj, ffj, jl)
1394 
1395  do n = 1, 5
1396  do j = 2, jl
1397  dw(i, j, k, n) = ffj(j, n)
1398  end do
1399  end do
1400 
1401  end do
1402  end do
1403 
1404  end if
1405  ! Multiply by T_xi^inv T_eta
1406 
1407  do k = 2, kl
1408  do j = 2, jl
1409  do i = 2, il
1410  ri1 = half * (si(i, j, k, 1) + si(i - 1, j, k, 1))
1411  ri2 = half * (si(i, j, k, 2) + si(i - 1, j, k, 2))
1412  ri3 = half * (si(i, j, k, 3) + si(i - 1, j, k, 3))
1413  ri = sqrt(ri1 * ri1 + ri2 * ri2 + ri3 * ri3)
1414  ri1 = ri1 / ri
1415  ri2 = ri2 / ri
1416  ri3 = ri3 / ri
1417 
1418  rj1 = half * (sj(i, j, k, 1) + sj(i, j - 1, k, 1))
1419  rj2 = half * (sj(i, j, k, 2) + sj(i, j - 1, k, 2))
1420  rj3 = half * (sj(i, j, k, 3) + sj(i, j - 1, k, 3))
1421  rj = sqrt(rj1 * rj1 + rj2 * rj2 + rj3 * rj3)
1422  rj1 = rj1 / rj
1423  rj2 = rj2 / rj
1424  rj3 = rj3 / rj
1425 
1426  dw1 = dw(i, j, k, 1)
1427  dw2 = dw(i, j, k, 2)
1428  dw3 = dw(i, j, k, 3)
1429  dw4 = dw(i, j, k, 4)
1430  dw5 = dw(i, j, k, 5)
1431 
1432  a1 = (ri1 * rj1 + ri2 * rj2 + ri3 * rj3)
1433  a2 = (ri1 * rj2 - rj1 * ri2)
1434  a3 = (ri3 * rj2 - rj3 * ri2)
1435  a4 = (ri1 * rj3 - rj1 * ri3)
1436  a5 = (dw4 - dw5) * sqrt2inv
1437  a6 = (dw4 + dw5) * half
1438  a7 = (a3 * dw1 + a4 * dw2 - a2 * dw3 - a5 * a1) * sqrt2inv
1439 
1440  dw(i, j, k, 1) = a1 * dw1 + a2 * dw2 + a4 * dw3 + a5 * a3
1441  dw(i, j, k, 2) = -a2 * dw1 + a1 * dw2 - a3 * dw3 + a5 * a4
1442  dw(i, j, k, 3) = -a4 * dw1 + a3 * dw2 + a1 * dw3 - a5 * a2
1443  dw(i, j, k, 4) = -a7 + a6
1444  dw(i, j, k, 5) = a7 + a6
1445 
1446  end do
1447  end do
1448  end do
1449 
1450  ! Multiply by diagonal
1451 
1452  do k = 2, kl
1453  do j = 2, jl
1454  do i = 2, il
1455  xfact = one + spectral_i(i, j, k) + spectral_j(i, j, k) &
1456  + spectral_k(i, j, k)
1457  dw(i, j, k, 1) = dw(i, j, k, 1) * xfact
1458  dw(i, j, k, 2) = dw(i, j, k, 2) * xfact
1459  dw(i, j, k, 3) = dw(i, j, k, 3) * xfact
1460  dw(i, j, k, 4) = dw(i, j, k, 4) * xfact
1461  dw(i, j, k, 5) = dw(i, j, k, 5) * xfact
1462  end do
1463  end do
1464  end do
1465 
1466  if (il .gt. 2) then
1467 
1468  ! Inversion in i
1469 
1470  do k = 2, kl
1471  do j = 2, jl
1472  do i = 1, il
1473  if (viscous) mut = rlv(i, j, k) + rlv(i + 1, j, k)
1474  if (eddymodel) mut = mut + rev(i, j, k) + rev(i + 1, j, k)
1475 
1476  volfact = one / (vol(i, j, k) + vol(i + 1, j, k))
1477  metterm = (si(i, j, k, 1) * si(i, j, k, 1) &
1478  + si(i, j, k, 2) * si(i, j, k, 2) &
1479  + si(i, j, k, 3) * si(i, j, k, 3))
1480  mettermi(i) = metterm * mut * volfact
1481  end do
1482 
1483  do i = 2, il
1484 
1485  viscterm1 = mettermi(i) / vol(i, j, k) / w(i, j, k, irho)
1486  viscterm3 = mettermi(i - 1) / vol(i, j, k) / w(i, j, k, irho)
1487  viscterm2 = viscterm1 + viscterm3
1488 
1489  volhalf = half / vol(i, j, k)
1490  ri1 = volhalf * (si(i, j, k, 1) + si(i - 1, j, k, 1))
1491  ri2 = volhalf * (si(i, j, k, 2) + si(i - 1, j, k, 2))
1492  ri3 = volhalf * (si(i, j, k, 3) + si(i - 1, j, k, 3))
1493  metterm = ri1 * ri1 + ri2 * ri2 + ri3 * ri3
1494  eps2 = epsval * epsval * cinf2 * metterm
1495  diagplus(1) = half * (qq_i(i, j, k) + fac * sqrt(qq_i(i, j, k)**2 + eps2))
1496  diagplus(2) = diagplus(1)
1497  diagplus(3) = diagplus(1)
1498  diagplus(4) = half * (qq_i(i, j, k) + cc_i(i, j, k) &
1499  + fac * sqrt((qq_i(i, j, k) + cc_i(i, j, k))**2 + eps2))
1500  diagplus(5) = half * (qq_i(i, j, k) - cc_i(i, j, k) &
1501  + fac * sqrt((qq_i(i, j, k) - cc_i(i, j, k))**2 + eps2))
1502  diagminus(1) = half * (qq_i(i, j, k) - fac * sqrt(qq_i(i, j, k)**2 + eps2))
1503  diagminus(2) = diagminus(1)
1504  diagminus(3) = diagminus(1)
1505  diagminus(4) = half * (qq_i(i, j, k) + cc_i(i, j, k) &
1506  - fac * sqrt((qq_i(i, j, k) + cc_i(i, j, k))**2 + eps2))
1507  diagminus(5) = half * (qq_i(i, j, k) - cc_i(i, j, k) &
1508  - fac * sqrt((qq_i(i, j, k) - cc_i(i, j, k))**2 + eps2))
1509  do n = 1, 5
1510  bbi(i + 1, n) = -viscterm1 - diagplus(n)
1511  ddi(i - 1, n) = -viscterm3 + diagminus(n)
1512  cci(i, n) = viscterm2 + diagplus(n) - diagminus(n)
1513  end do
1514  end do
1515 
1516  do n = 1, 5
1517  bbi(ie, n) = zero
1518  ddi(1, n) = zero
1519  do i = 2, il
1520  bbi(i, n) = bbi(i, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1521  ddi(i, n) = ddi(i, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1522  cci(i, n) = one + cci(i, n) * dual_dt(i, j, k) * &
1523  max(real(iblank(i, j, k), realtype), zero) + &
1524  spectral_j(i, j, k) + spectral_k(i, j, k)
1525  ffi(i, n) = dw(i, j, k, n)
1526  end do
1527  end do
1528 
1529  call tridiagsolve(bbi, cci, ddi, ffi, il)
1530 
1531  do n = 1, 5
1532  do i = 2, il
1533  dw(i, j, k, n) = ffi(i, n)
1534  end do
1535  end do
1536 
1537  end do
1538  end do
1539 
1540  end if
1541  ! Multiply by T_zeta^inv T_xi
1542 
1543  do k = 2, kl
1544  do j = 2, jl
1545  do i = 2, il
1546  ri1 = half * (si(i, j, k, 1) + si(i - 1, j, k, 1))
1547  ri2 = half * (si(i, j, k, 2) + si(i - 1, j, k, 2))
1548  ri3 = half * (si(i, j, k, 3) + si(i - 1, j, k, 3))
1549  ri = sqrt(ri1 * ri1 + ri2 * ri2 + ri3 * ri3)
1550  ri1 = ri1 / ri
1551  ri2 = ri2 / ri
1552  ri3 = ri3 / ri
1553 
1554  rk1 = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1))
1555  rk2 = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2))
1556  rk3 = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3))
1557  rk = sqrt(rk1 * rk1 + rk2 * rk2 + rk3 * rk3)
1558  rk1 = rk1 / rk
1559  rk2 = rk2 / rk
1560  rk3 = rk3 / rk
1561 
1562  dw1 = dw(i, j, k, 1)
1563  dw2 = dw(i, j, k, 2)
1564  dw3 = dw(i, j, k, 3)
1565  dw4 = dw(i, j, k, 4)
1566  dw5 = dw(i, j, k, 5)
1567 
1568  a1 = ri1 * rk1 + ri2 * rk2 + ri3 * rk3
1569  a2 = rk1 * ri2 - ri1 * rk2
1570  a3 = rk3 * ri2 - ri3 * rk2
1571  a4 = rk1 * ri3 - ri1 * rk3
1572  a5 = (dw4 - dw5) * sqrt2inv
1573  a6 = (dw4 + dw5) * half
1574  a7 = (a3 * dw1 + a4 * dw2 - a2 * dw3 - a5 * a1) * sqrt2inv
1575 
1576  dw(i, j, k, 1) = a1 * dw1 + a2 * dw2 + a4 * dw3 + a5 * a3
1577  dw(i, j, k, 2) = -a2 * dw1 + a1 * dw2 - a3 * dw3 + a5 * a4
1578  dw(i, j, k, 3) = -a4 * dw1 + a3 * dw2 + a1 * dw3 - a5 * a2
1579  dw(i, j, k, 4) = -a7 + a6
1580  dw(i, j, k, 5) = a7 + a6
1581  end do
1582  end do
1583  end do
1584 
1585  ! Multiply by diagonal
1586 
1587  do k = 2, kl
1588  do j = 2, jl
1589  do i = 2, il
1590  xfact = one + spectral_i(i, j, k) + spectral_j(i, j, k) &
1591  + spectral_k(i, j, k)
1592  dw(i, j, k, 1) = dw(i, j, k, 1) * xfact
1593  dw(i, j, k, 2) = dw(i, j, k, 2) * xfact
1594  dw(i, j, k, 3) = dw(i, j, k, 3) * xfact
1595  dw(i, j, k, 4) = dw(i, j, k, 4) * xfact
1596  dw(i, j, k, 5) = dw(i, j, k, 5) * xfact
1597  end do
1598  end do
1599  end do
1600 
1601  if (kl .gt. 2) then
1602 
1603  ! Inversion in k
1604 
1605  do j = 2, jl
1606  do i = 2, il
1607  do k = 1, kl
1608  if (viscous) mut = rlv(i, j, k) + rlv(i, j, k + 1)
1609  if (eddymodel) mut = mut + rev(i, j, k) + rev(i, j, k + 1)
1610 
1611  volfact = 1./(vol(i, j, k) + vol(i, j, k + 1))
1612  metterm = sk(i, j, k, 1) * sk(i, j, k, 1) &
1613  + sk(i, j, k, 2) * sk(i, j, k, 2) &
1614  + sk(i, j, k, 3) * sk(i, j, k, 3)
1615  mettermk(k) = metterm * mut * volfact
1616  end do
1617 
1618  do k = 2, kl
1619 
1620  viscterm1 = mettermk(k) / vol(i, j, k) / w(i, j, k, irho)
1621  viscterm3 = mettermk(k - 1) / vol(i, j, k) / w(i, j, k, irho)
1622  viscterm2 = viscterm1 + viscterm3
1623 
1624  volhalf = half / vol(i, j, k)
1625  rk1 = volhalf * (sk(i, j, k, 1) + sj(i, j, k - 1, 1))
1626  rk2 = volhalf * (sk(i, j, k, 2) + sj(i, j, k - 1, 2))
1627  rk3 = volhalf * (sk(i, j, k, 3) + sj(i, j, k - 1, 3))
1628  metterm = rk1 * rk1 + rk2 * rk2 + rk3 * rk3
1629  eps2 = epsval * epsval * cinf2 * metterm
1630  diagplus(1) = half * (qq_k(i, j, k) + fac * sqrt(qq_k(i, j, k)**2 + eps2))
1631  diagplus(2) = diagplus(1)
1632  diagplus(3) = diagplus(1)
1633  diagplus(4) = half * (qq_k(i, j, k) + cc_k(i, j, k) &
1634  + fac * sqrt((qq_k(i, j, k) + cc_k(i, j, k))**2 + eps2))
1635  diagplus(5) = half * (qq_k(i, j, k) - cc_k(i, j, k) &
1636  + fac * sqrt((qq_k(i, j, k) - cc_k(i, j, k))**2 + eps2))
1637  diagminus(1) = half * (qq_k(i, j, k) - fac * sqrt(qq_k(i, j, k)**2 + eps2))
1638  diagminus(2) = diagminus(1)
1639  diagminus(3) = diagminus(1)
1640  diagminus(4) = half * (qq_k(i, j, k) + cc_k(i, j, k) &
1641  - fac * sqrt((qq_k(i, j, k) + cc_k(i, j, k))**2 + eps2))
1642  diagminus(5) = half * (qq_k(i, j, k) - cc_k(i, j, k) &
1643  - fac * sqrt((qq_k(i, j, k) - cc_k(i, j, k))**2 + eps2))
1644 
1645  do n = 1, 5
1646  bbk(k + 1, n) = -viscterm1 - diagplus(n)
1647  ddk(k - 1, n) = -viscterm3 + diagminus(n)
1648  cck(k, n) = viscterm2 + diagplus(n) - diagminus(n)
1649  end do
1650  end do
1651 
1652  do n = 1, 5
1653  bbk(ke, n) = zero
1654  ddk(1, n) = zero
1655  do k = 2, kl
1656  bbk(k, n) = bbk(k, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1657  ddk(k, n) = ddk(k, n) * dual_dt(i, j, k) * max(real(iblank(i, j, k), realtype), zero)
1658  cck(k, n) = one + cck(k, n) * dual_dt(i, j, k) * &
1659  max(real(iblank(i, j, k), realtype), zero) + &
1660  spectral_i(i, j, k) + spectral_j(i, j, k)
1661  ffk(k, n) = dw(i, j, k, n)
1662  end do
1663  end do
1664 
1665  call tridiagsolve(bbk, cck, ddk, ffk, kl)
1666 
1667  do n = 1, 5
1668  do k = 2, kl
1669  dw(i, j, k, n) = ffk(k, n)
1670  end do
1671  end do
1672 
1673  end do
1674  end do
1675  end if
1676 
1677  ! Multiply by T_zeta
1678 
1679  do k = 2, kl
1680  do j = 2, jl
1681  do i = 2, il
1682  uvel = w(i, j, k, ivx)
1683  vvel = w(i, j, k, ivy)
1684  wvel = w(i, j, k, ivz)
1685 
1686  rk1 = half * (sk(i, j, k, 1) + sk(i, j, k - 1, 1))
1687  rk2 = half * (sk(i, j, k, 2) + sk(i, j, k - 1, 2))
1688  rk3 = half * (sk(i, j, k, 3) + sk(i, j, k - 1, 3))
1689 
1690  rk = sqrt(rk1 * rk1 + rk2 * rk2 + rk3 * rk3)
1691  uu = uvel * rk1 + vvel * rk2 + wvel * rk3
1692 
1693  rk1 = rk1 / rk
1694  rk2 = rk2 / rk
1695  rk3 = rk3 / rk
1696 
1697  uvw = half * (uvel * uvel + vvel * vvel + wvel * wvel)
1698  cijkinv = sqrt(w(i, j, k, irho) / gamma(i, j, k) / p(i, j, k))
1699  alph = w(i, j, k, irho) * cijkinv * sqrt2inv
1700  xfact = two / cijkinv
1701 
1702  ge = gamma(i, j, k) * w(i, j, k, irhoe) / w(i, j, k, irho) - &
1703  (gamma(i, j, k) - one) * uvw
1704 
1705  dw1 = dw(i, j, k, 1)
1706  dw2 = dw(i, j, k, 2)
1707  dw3 = dw(i, j, k, 3)
1708  dw4 = dw(i, j, k, 4) * alph
1709  dw5 = dw(i, j, k, 5) * alph
1710 
1711  a1 = dw1 * rk1 + dw2 * rk2 + dw3 * rk3 + dw4 + dw5
1712  a2 = half * xfact * (dw4 - dw5)
1713  a3 = uvw * (rk1 * dw1 + rk2 * dw2 + rk3 * dw3)
1714 
1715  dw(i, j, k, 1) = a1
1716  dw(i, j, k, 2) = a1 * uvel - w(i, j, k, irho) * (rk3 * dw2 - rk2 * dw3) &
1717  + a2 * rk1
1718  dw(i, j, k, 3) = a1 * vvel - w(i, j, k, irho) * (rk1 * dw3 - rk3 * dw1) &
1719  + a2 * rk2
1720  dw(i, j, k, 4) = a1 * wvel - w(i, j, k, irho) * (rk2 * dw1 - rk1 * dw2) &
1721  + a2 * rk3
1722  dw(i, j, k, 5) = a3 + w(i, j, k, irho) * &
1723  ((vvel * rk3 - wvel * rk2) * dw1 &
1724  + (wvel * rk1 - uvel * rk3) * dw2 &
1725  + (uvel * rk2 - vvel * rk1) * dw3) &
1726  + (ge + half * xfact * uu / rk) * dw4 &
1727  + (ge - half * xfact * uu / rk) * dw5
1728 
1729  end do
1730  end do
1731  end do
1732 
1733  ! For consistency with update.
1734 
1735  do k = 2, kl
1736  do j = 2, jl
1737  do i = 2, il
1738  volfact = -one / vol(i, j, k)
1739  dw(i, j, k, 1) = dw(i, j, k, 1) * volfact
1740  dw(i, j, k, 2) = dw(i, j, k, 2) * volfact
1741  dw(i, j, k, 3) = dw(i, j, k, 3) * volfact
1742  dw(i, j, k, 4) = dw(i, j, k, 4) * volfact
1743  dw(i, j, k, 5) = dw(i, j, k, 5) * volfact
1744  end do
1745  end do
1746  end do
1747 
1748  end subroutine computedwdadi
1749 
1750  subroutine tridiagsolve(bb, cc, dd, ff, nn)
1751  use precision
1752  implicit none
1753  !
1754  ! Subroutine arguments
1755  !
1756  integer(kind=intType) :: nn
1757  real(kind=realtype), dimension(nn + 1, 5) :: bb, cc, dd, ff
1758 
1759  ! local variables
1760 
1761  integer(kind=intType) :: m, n
1762  real(kind=realtype) :: d0, d2
1763 
1764  do n = 1, 5
1765  m = 2
1766  d0 = 1./cc(m, n)
1767  dd(m, n) = dd(m, n) * d0
1768  ff(m, n) = ff(m, n) * d0
1769 
1770  do m = 3, nn
1771  d2 = bb(m, n)
1772  d0 = 1./(cc(m, n) - d2 * dd(m - 1, n))
1773  ff(m, n) = (ff(m, n) - d2 * ff(m - 1, n)) * d0
1774  dd(m, n) = dd(m, n) * d0
1775  end do
1776 
1777  do m = nn - 1, 2, -1
1778  ff(m, n) = ff(m, n) - dd(m, n) * ff(m + 1, n)
1779  end do
1780 
1781  end do
1782 
1783  end subroutine tridiagsolve
1784 
1785  subroutine residualaveraging
1786  !
1787  ! Implicit residual smoothing is a simple procedure that
1788  ! replaces the residual at each point by a weighted sum of all
1789  ! of the residuals in the block (although the residuals that are
1790  ! closer to the cell under consideration are weighted more
1791  ! heavily). This smoothing can be applied explicitly, but may
1792  ! result in zero smoothed residuals for non-zero initial
1793  ! residual modes. For this reason, the smoothing is applied
1794  ! implicitly in the following form:
1795  ! -epz R{i+1} + (1 + 2 epz) R{i} -epz R{i-1} = r{i}
1796  ! Where r{i} is the original residual at point i, and R{i} is
1797  ! the implicitly smoothed residual at point i. The analysis for
1798  ! the 1-D scalar convection-diffusion equation shows that if
1799  ! Epz >= (1/4)*((lambda/lambda*)^2 -1),
1800  ! where lambda is the cfl number desired to be used, and
1801  ! lambda* is the CFL limit of the unsmoothed scheme, the scheme
1802  ! can be made unconditionally stable (arbitrarily large lambda).
1803  ! In practice, lambda = 6-8 is common for Euler solutions. For
1804  ! RANS solutions lambda = 3-4 is what we have used in practice
1805  ! resulting in a slight improvement in convergence rate, but,
1806  ! more importantly, an increase in the robustness of the solver.
1807  ! Note that this theoretical result can be shown for infinite
1808  ! 1-D problems, but also for finite-periodic 1-D problems and
1809  ! finite-size 1-D problems (i.e. with block boundaries). Such
1810  ! theoretical results are not available for 3-D cases, where the
1811  ! scheme is applied in an ADI fashion:
1812  ! (1 -epzI d_ii)(1 -epzJ d_jj)(1 -epzK d_kk) r = r
1813  ! Where d_ii, d_jj, d_kk are second difference operators in each
1814  ! of the mesh coordinate directions.
1815  ! For each of the coordinate direction solves, the initial
1816  ! matrix problem is of the form:
1817  ! - - - - - -
1818  ! | (1+2 epz) -epz | |r| |r|
1819  ! | -epz (1 +2 epz) -epz | |r| |r|
1820  ! | -epz (1 + 2 epz) -epz | |r| = |r|
1821  ! | . . . | |.| |.|
1822  ! | . . . | |.| |.|
1823  ! - - - - - -
1824  ! And after the forward elimination phase a normalization is
1825  ! applied so the result looks like:
1826  ! - - - - - -
1827  ! | 1 -d | |r| |r|
1828  ! | 0 1 -d | |r| |r|
1829  ! | 0 1 -d | |r| = |r|
1830  ! | . . . | |.| |.|
1831  ! | . . . | |.| |.|
1832  ! - - - - - -
1833  ! Which can then be used with a straightforward backsolve to
1834  ! obtain the answer.
1835  ! It is assumed that the pointers in blockPointers already
1836  ! point to the correct block.
1837  !
1838  use constants
1839  use blockpointers, only: nx, ny, nz, il, jl, kl, dw, p, iblank
1840  use inputiteration, only: cfl, cflcoarse, cfllimit, smoop
1841  use flowvarrefstate, only: pinfcorr, nwf
1842  use iteration, only: currentlevel
1843 
1844  implicit none
1845  !
1846  ! Local variables.
1847  !
1848  integer(kind=intType) :: i, j, k, l
1849 
1850  real(kind=realtype), parameter :: b = 2.0_realtype
1851 
1852  real(kind=realtype) :: currentcfl, rfl0, plim
1853  real(kind=realtype) :: dpi, dpj, dpk, r
1854  real(kind=realtype), dimension(il, max(jl, kl)) :: epz, d, t, rfl
1855 
1856  ! rfl0 is a measure of the ratio lambda/lambda*
1857  !
1858  currentcfl = cflcoarse
1859  if (currentlevel == 1) then
1860  currentcfl = cfl
1861  end if
1862 
1863  rfl0 = half * currentcfl / cfllimit
1864 
1865  plim = 0.001_realtype * pinfcorr
1866 
1867  ! Smoothing in the i-direction. Only done when enough cells are
1868  ! present in the i-direction.
1869  !
1870  if (nx > 1) then
1871 
1872  do k = 2, kl
1873 
1874  ! Compute smoothing coeficients
1875 
1876  do j = 2, jl
1877  do i = 2, il
1878  dpi = abs(p(i + 1, j, k) - two * p(i, j, k) + p(i - 1, j, k)) &
1879  / (p(i + 1, j, k) + two * p(i, j, k) + p(i - 1, j, k) + plim)
1880  dpj = abs(p(i, j + 1, k) - two * p(i, j, k) + p(i, j - 1, k)) &
1881  / (p(i, j + 1, k) + two * p(i, j, k) + p(i, j - 1, k) + plim)
1882  dpk = abs(p(i, j, k + 1) - two * p(i, j, k) + p(i, j, k - 1)) &
1883  / (p(i, j, k + 1) + two * p(i, j, k) + p(i, j, k - 1) + plim)
1884  rfl(i, j) = one / (one + b * (dpi + dpj + dpk))
1885  end do
1886  end do
1887 
1888  do j = 2, jl
1889  do i = 2, nx
1890  r = rfl0 * (rfl(i, j) + rfl(i + 1, j))
1891  epz(i, j) = fourth * smoop * dim(r * r, one) * max(real(iblank(i, j, k), realtype), zero)
1892  end do
1893  end do
1894 
1895  ! Zero out coefficients for boundary condition treatment
1896 
1897  do j = 2, jl
1898  epz(1, j) = zero
1899  epz(il, j) = zero
1900  d(1, j) = zero
1901  end do
1902 
1903  ! Compute coefficients for forward elimination process
1904 
1905  do i = 2, il
1906  do j = 2, jl
1907  t(i, j) = one &
1908  / (one + epz(i, j) + epz(i - 1, j) - epz(i - 1, j) * d(i - 1, j))
1909  d(i, j) = t(i, j) * epz(i, j)
1910  end do
1911  end do
1912 
1913  ! Apply same transformation to the rhs vector of residuals
1914 
1915  do i = 2, il
1916  do j = 2, jl
1917  do l = 1, nwf
1918  dw(i, j, k, l) = t(i, j) &
1919  * (dw(i, j, k, l) + epz(i - 1, j) * dw(i - 1, j, k, l))
1920  end do
1921  end do
1922  end do
1923 
1924  ! Backsolve operation. Smoothed residuals are left
1925  ! in dw(i,j,k,l)
1926 
1927  do i = nx, 2, -1
1928  do j = 2, jl
1929  do l = 1, nwf
1930  dw(i, j, k, l) = dw(i, j, k, l) + d(i, j) * dw(i + 1, j, k, l)
1931  end do
1932  end do
1933  end do
1934 
1935  end do
1936  end if
1937  !
1938  ! Smoothing in the j-direction. Only done when enough cells are
1939  ! present in the j-direction.
1940  !
1941  if (ny > 1) then
1942 
1943  do k = 2, kl
1944 
1945  ! Compute smoothing coeficients
1946 
1947  do j = 2, jl
1948  do i = 2, il
1949  dpi = abs(p(i + 1, j, k) - two * p(i, j, k) + p(i - 1, j, k)) &
1950  / (p(i + 1, j, k) + two * p(i, j, k) + p(i - 1, j, k) + plim)
1951  dpj = abs(p(i, j + 1, k) - two * p(i, j, k) + p(i, j - 1, k)) &
1952  / (p(i, j + 1, k) + two * p(i, j, k) + p(i, j - 1, k) + plim)
1953  dpk = abs(p(i, j, k + 1) - two * p(i, j, k) + p(i, j, k - 1)) &
1954  / (p(i, j, k + 1) + two * p(i, j, k) + p(i, j, k - 1) + plim)
1955  rfl(i, j) = one / (one + b * (dpi + dpj + dpk))
1956  end do
1957  end do
1958 
1959  do j = 2, ny
1960  do i = 2, il
1961  r = rfl0 * (rfl(i, j) + rfl(i, j + 1))
1962  epz(i, j) = fourth * smoop * dim(r * r, one) * max(real(iblank(i, j, k), realtype), zero)
1963  end do
1964  end do
1965 
1966  ! Zero out coefficients for boundary condition treatment
1967 
1968  do i = 2, il
1969  epz(i, 1) = zero
1970  epz(i, jl) = zero
1971  d(i, 1) = zero
1972  end do
1973 
1974  ! Compute coefficients for forward eliMination process
1975 
1976  do j = 2, jl
1977  do i = 2, il
1978  t(i, j) = one &
1979  / (one + epz(i, j) + epz(i, j - 1) - epz(i, j - 1) * d(i, j - 1))
1980  d(i, j) = t(i, j) * epz(i, j)
1981  end do
1982  end do
1983 
1984  ! Apply same transformation to the rhs vector of residuals
1985 
1986  do j = 2, jl
1987  do i = 2, il
1988  do l = 1, nwf
1989  dw(i, j, k, l) = t(i, j) &
1990  * (dw(i, j, k, l) + epz(i, j - 1) * dw(i, j - 1, k, l))
1991  end do
1992  end do
1993  end do
1994 
1995  ! Backsolve operation. Smoothed residuals are left
1996  ! in dw(i,j,k,l)
1997 
1998  do j = ny, 2, -1
1999  do i = 2, il
2000  do l = 1, nwf
2001  dw(i, j, k, l) = dw(i, j, k, l) + d(i, j) * dw(i, j + 1, k, l)
2002  end do
2003  end do
2004  end do
2005 
2006  end do
2007  end if
2008  !
2009  ! Smoothing in the k-direction. Only done when enough cells are
2010  ! present in the k-direction.
2011  !
2012  if (nz > 1) then
2013 
2014  do j = 2, jl
2015 
2016  ! Compute smoothing coeficients
2017 
2018  do k = 2, kl
2019  do i = 2, il
2020  dpi = abs(p(i + 1, j, k) - two * p(i, j, k) + p(i - 1, j, k)) &
2021  / (p(i + 1, j, k) + two * p(i, j, k) + p(i - 1, j, k) + plim)
2022  dpj = abs(p(i, j + 1, k) - two * p(i, j, k) + p(i, j - 1, k)) &
2023  / (p(i, j + 1, k) + two * p(i, j, k) + p(i, j - 1, k) + plim)
2024  dpk = abs(p(i, j, k + 1) - two * p(i, j, k) + p(i, j, k - 1)) &
2025  / (p(i, j, k + 1) + two * p(i, j, k) + p(i, j, k - 1) + plim)
2026  rfl(i, k) = one / (one + b * (dpi + dpj + dpk))
2027  end do
2028  end do
2029 
2030  do k = 2, nz
2031  do i = 2, il
2032  r = rfl0 * (rfl(i, k) + rfl(i, k + 1))
2033  epz(i, k) = fourth * smoop * dim(r * r, one) * max(real(iblank(i, j, k), realtype), zero)
2034  end do
2035  end do
2036 
2037  ! Zero out coefficients for boundary condition treatment
2038 
2039  do i = 2, il
2040  epz(i, 1) = zero
2041  epz(i, kl) = zero
2042  d(i, 1) = zero
2043  end do
2044 
2045  ! Compute coefficients for forward eliMination process
2046 
2047  do k = 2, kl
2048  do i = 2, il
2049  t(i, k) = one &
2050  / (one + epz(i, k) + epz(i, k - 1) - epz(i, k - 1) * d(i, k - 1))
2051  d(i, k) = t(i, k) * epz(i, k)
2052  end do
2053  end do
2054 
2055  ! Apply same transformation to the rhs vector of residuals
2056 
2057  do k = 2, kl
2058  do i = 2, il
2059  do l = 1, nwf
2060  dw(i, j, k, l) = t(i, k) &
2061  * (dw(i, j, k, l) + epz(i, k - 1) * dw(i, j, k - 1, l))
2062  end do
2063  end do
2064  end do
2065 
2066  ! Backsolve operation. Smoothed residuals are left
2067  ! in dw(i,j,k,l)
2068 
2069  do k = nz, 2, -1
2070  do i = 2, il
2071  do l = 1, nwf
2072  dw(i, j, k, l) = dw(i, j, k, l) + d(i, k) * dw(i, j, k + 1, l)
2073  end do
2074  end do
2075  end do
2076 
2077  end do
2078  end if
2079 
2080  end subroutine residualaveraging
2081 #endif
2082 end module residuals
actuatorregiondata::actuatorregions
type(actuatorregiontype), dimension(nactuatorregionsmax), target actuatorregions
Definition: actuatorRegionData.F90:39
actuatorregiondata::nactuatorregions
integer(kind=inttype) nactuatorregions
Definition: actuatorRegionData.F90:40
aleutils::recoverlevelale_block
subroutine recoverlevelale_block
Definition: ALEUtils.F90:358
aleutils::interplevelale_block
subroutine interplevelale_block
Definition: ALEUtils.F90:245
blockpointers::sfacek
real(kind=realtype), dimension(:, :, :), pointer sfacek
Definition: blockPointers.F90:119
blockpointers::w1
real(kind=realtype), dimension(:, :, :, :), pointer w1
Definition: blockPointers.F90:131
blockpointers::gamma
real(kind=realtype), dimension(:, :, :), pointer gamma
Definition: blockPointers.F90:123
blockpointers::volold
real(kind=realtype), dimension(:, :, :, :), pointer volold
Definition: blockPointers.F90:99
blockpointers::addgridvelocities
logical addgridvelocities
Definition: blockPointers.F90:117
blockpointers::jl
integer(kind=inttype) jl
Definition: blockPointers.F90:38
blockpointers::wr
real(kind=realtype), dimension(:, :, :, :), pointer wr
Definition: blockPointers.F90:131
blockpointers::nx
integer(kind=inttype) nx
Definition: blockPointers.F90:38
blockpointers::p
real(kind=realtype), dimension(:, :, :), pointer p
Definition: blockPointers.F90:123
blockpointers::ny
integer(kind=inttype) ny
Definition: blockPointers.F90:38
blockpointers::ie
integer(kind=inttype) ie
Definition: blockPointers.F90:39
blockpointers::w
real(kind=realtype), dimension(:, :, :, :), pointer w
Definition: blockPointers.F90:120
blockpointers::scratch
real(kind=realtype), dimension(:, :, :, :), pointer scratch
Definition: blockPointers.F90:129
blockpointers::sfacei
real(kind=realtype), dimension(:, :, :), pointer sfacei
Definition: blockPointers.F90:119
blockpointers::iblank
integer(kind=inttype), dimension(:, :, :), pointer iblank
Definition: blockPointers.F90:71
blockpointers::wn
real(kind=realtype), dimension(:, :, :, :), pointer wn
Definition: blockPointers.F90:149
blockpointers::jb
integer(kind=inttype) jb
Definition: blockPointers.F90:39
blockpointers::kb
integer(kind=inttype) kb
Definition: blockPointers.F90:39
blockpointers::rlv
real(kind=realtype), dimension(:, :, :), pointer rlv
Definition: blockPointers.F90:125
blockpointers::si
real(kind=realtype), dimension(:, :, :, :), pointer si
Definition: blockPointers.F90:96
blockpointers::sj
real(kind=realtype), dimension(:, :, :, :), pointer sj
Definition: blockPointers.F90:96
blockpointers::sectionid
integer(kind=inttype) sectionid
Definition: blockPointers.F90:30
blockpointers::rev
real(kind=realtype), dimension(:, :, :), pointer rev
Definition: blockPointers.F90:125
blockpointers::dw
real(kind=realtype), dimension(:, :, :, :), pointer dw
Definition: blockPointers.F90:128
blockpointers::sk
real(kind=realtype), dimension(:, :, :, :), pointer sk
Definition: blockPointers.F90:96
blockpointers::vol
real(kind=realtype), dimension(:, :, :), pointer vol
Definition: blockPointers.F90:97
blockpointers::ib
integer(kind=inttype) ib
Definition: blockPointers.F90:39
blockpointers::dtl
real(kind=realtype), dimension(:, :, :), pointer dtl
Definition: blockPointers.F90:151
blockpointers::sfacej
real(kind=realtype), dimension(:, :, :), pointer sfacej
Definition: blockPointers.F90:119
blockpointers::fw
real(kind=realtype), dimension(:, :, :, :), pointer fw
Definition: blockPointers.F90:128
blockpointers::nz
integer(kind=inttype) nz
Definition: blockPointers.F90:38
blockpointers::je
integer(kind=inttype) je
Definition: blockPointers.F90:39
blockpointers::ke
integer(kind=inttype) ke
Definition: blockPointers.F90:39
blockpointers::wold
real(kind=realtype), dimension(:, :, :, :, :), pointer wold
Definition: blockPointers.F90:121
blockpointers::kl
integer(kind=inttype) kl
Definition: blockPointers.F90:38
blockpointers::il
integer(kind=inttype) il
Definition: blockPointers.F90:38
constants::zero
real(kind=realtype), parameter zero
Definition: constants.F90:71
constants::irho
integer, parameter irho
Definition: constants.F90:34
constants::imx
integer, parameter imx
Definition: constants.F90:65
constants::ivx
integer, parameter ivx
Definition: constants.F90:35
constants::irhoe
integer, parameter irhoe
Definition: constants.F90:38
constants::one
real(kind=realtype), parameter one
Definition: constants.F90:72
constants::half
real(kind=realtype), parameter half
Definition: constants.F90:80
constants::steady
integer(kind=inttype), parameter steady
Definition: constants.F90:115
constants::ivz
integer, parameter ivz
Definition: constants.F90:37
constants::two
real(kind=realtype), parameter two
Definition: constants.F90:73
constants::fourth
real(kind=realtype), parameter fourth
Definition: constants.F90:82
constants::imz
integer, parameter imz
Definition: constants.F90:67
constants::ivy
integer, parameter ivy
Definition: constants.F90:36
flowutils::allnodalgradients
subroutine allnodalgradients
Definition: flowUtils.F90:1677
flowutils::computespeedofsoundsquared
subroutine computespeedofsoundsquared
Definition: flowUtils.F90:489
flowvarrefstate::gammainf
real(kind=realtype) gammainf
Definition: flowVarRefState.F90:70
flowvarrefstate::pinfcorr
real(kind=realtype) pinfcorr
Definition: flowVarRefState.F90:69
flowvarrefstate::pinf
real(kind=realtype) pinf
Definition: flowVarRefState.F90:69
flowvarrefstate::uref
real(kind=realtype) uref
Definition: flowVarRefState.F90:40
flowvarrefstate::nwf
integer(kind=inttype) nwf
Definition: flowVarRefState.F90:17
flowvarrefstate::winf
real(kind=realtype), dimension(:), allocatable winf
Definition: flowVarRefState.F90:74
flowvarrefstate::lref
real(kind=realtype) lref
Definition: flowVarRefState.F90:48
flowvarrefstate::rhoinf
real(kind=realtype) rhoinf
Definition: flowVarRefState.F90:69
flowvarrefstate::pref
real(kind=realtype) pref
Definition: flowVarRefState.F90:39
flowvarrefstate::timeref
real(kind=realtype) timeref
Definition: flowVarRefState.F90:40
fluxes
Definition: fluxes.F90:1
fluxes::invisciddissfluxmatrixcoarse
subroutine invisciddissfluxmatrixcoarse
Definition: fluxes.F90:5206
fluxes::invisciddissfluxmatrixapprox
subroutine invisciddissfluxmatrixapprox
Definition: fluxes.F90:4345
fluxes::invisciddissfluxscalarcoarse
subroutine invisciddissfluxscalarcoarse
Definition: fluxes.F90:4978
fluxes::invisciddissfluxscalar
subroutine invisciddissfluxscalar
Definition: fluxes.F90:1050
fluxes::invisciddissfluxmatrix
subroutine invisciddissfluxmatrix
Definition: fluxes.F90:404
fluxes::inviscidcentralflux
subroutine inviscidcentralflux
Definition: fluxes.F90:5
fluxes::viscousfluxapprox
subroutine viscousfluxapprox
Definition: fluxes.F90:3488
fluxes::inviscidupwindflux
subroutine inviscidupwindflux(fineGrid)
Definition: fluxes.F90:1439
fluxes::viscousflux
subroutine viscousflux
Definition: fluxes.F90:2535
fluxes::invisciddissfluxscalarapprox
subroutine invisciddissfluxscalarapprox
Definition: fluxes.F90:3862
inputadjoint::viscpc
logical viscpc
Definition: inputParam.F90:793
inputdiscretization::lowspeedpreconditioner
logical lowspeedpreconditioner
Definition: inputParam.F90:96
inputdiscretization::spacediscrcoarse
integer(kind=inttype) spacediscrcoarse
Definition: inputParam.F90:72
inputdiscretization::spacediscr
integer(kind=inttype) spacediscr
Definition: inputParam.F90:72
inputiteration::smoother
integer(kind=inttype) smoother
Definition: inputParam.F90:264
inputiteration::cflcoarse
real(kind=realtype) cflcoarse
Definition: inputParam.F90:275
inputiteration::cdisrk
real(kind=realtype), dimension(:), allocatable cdisrk
Definition: inputParam.F90:283
inputiteration::smoop
real(kind=realtype) smoop
Definition: inputParam.F90:275
inputiteration::cfllimit
real(kind=realtype) cfllimit
Definition: inputParam.F90:268
inputiteration::cfl
real(kind=realtype) cfl
Definition: inputParam.F90:275
inputphysics::equationmode
integer(kind=inttype) equationmode
Definition: inputParam.F90:583
inputtimespectral::dvector
real(kind=realtype), dimension(:, :, :), allocatable dvector
Definition: inputParam.F90:660
inputtimespectral::ntimeintervalsspectral
integer(kind=inttype) ntimeintervalsspectral
Definition: inputParam.F90:645
inputtimespectral::dscalar
real(kind=realtype), dimension(:, :, :), allocatable dscalar
Definition: inputParam.F90:659
inputunsteady::deltat
real(kind=realtype) deltat
Definition: inputParam.F90:733
inputunsteady::timeintegrationscheme
integer(kind=inttype) timeintegrationscheme
Definition: inputParam.F90:719
iteration::noldlevels
integer(kind=inttype) noldlevels
Definition: iteration.f90:79
iteration::currentlevel
integer(kind=inttype) currentlevel
Definition: iteration.f90:18
iteration::groundlevel
integer(kind=inttype) groundlevel
Definition: iteration.f90:18
iteration::rkstage
integer(kind=inttype) rkstage
Definition: iteration.f90:19
iteration::coeftime
real(kind=realtype), dimension(:), allocatable coeftime
Definition: iteration.f90:80
iteration::rfil
real(kind=realtype) rfil
Definition: iteration.f90:36
iteration::deforming_grid
logical deforming_grid
Definition: iteration.f90:69
iteration::ordersconverged
real(kind=realtype) ordersconverged
Definition: iteration.f90:128
precision::realtype
integer, parameter realtype
Definition: precision.F90:112
residuals::residualaveraging
subroutine residualaveraging
Definition: residuals.F90:1786
residuals::sourceterms
subroutine sourceterms
Definition: residuals.F90:1004
residuals::computedwdadi
subroutine computedwdadi
Definition: residuals.F90:1063
residuals::sourceterms_block
subroutine sourceterms_block(nn, res, iRegion, pLocal)
Definition: residuals.F90:349
residuals::initres_block
subroutine initres_block(varStart, varEnd, nn, sps)
Definition: residuals.F90:428
residuals::residual
subroutine residual
Definition: residuals.F90:1029
residuals::initres
subroutine initres(varStart, varEnd)
Definition: residuals.F90:965
residuals::residual_block
subroutine residual_block
Definition: residuals.F90:5
residuals::tridiagsolve
subroutine tridiagsolve(bb, cc, dd, ff, nn)
Definition: residuals.F90:1751
utils
Definition: utils.F90:1
utils::setpointers
subroutine setpointers(nn, mm, ll)
Definition: utils.F90:3237
Generated by doxygen 1.9.1