!----------------------------------------------------------------------
!
! Compressible/Incompressible stretched vortex subgrid model (LES)
!
!----------------------------------------------------------------------
!
! Copyright (C) California Institute of Technology 2000-2005
! Authors: D.I. Pullin (dale@galcit.caltech.edu)
! D.J. Hill
! C. Pantano
!
!----------------------------------------------------------------------
! USAGE
!
! x,y,z dimensions are (ixlo:ixhi,iylo:iyhi,izlo:izhi)
!
! ux(ncomps,,,) conservative vector of state with components
! rho, rho u, rho v, rho w, E, rho Y1, rho Y2, ...
! vx(nvars,,,) primitive vector of state with components
! rho, u, v, w, p, Y1, Y2,...
!
! Subroutines -------------------------------------------------------
!
! SetUpLES(version)
! needs to be called only once during the whole program
! version: 0 incompressible, 1 compressible
!
! AllocateLES/DeallocateLES() allocates/deallocates working storage
! call after previous function
!
! InitializeLES(ux,vx)
! initializes the sgs model. Should be called at the begining of each time
! step or stage if multistage method
!
! AddSgsFlux(fxi,ux,vx,dir)
! computes and adds the fluxes on each cell. This subroutine is specifc to
! our AMR fluxed based method. One should alter the end of each direction
! switch to reflect your method.
! fxi fluxes
! dir directions of the flux: 1,2,3
!-----------------------------------------------------------------------
MODULE LES_Spiral
SAVE
include 'lessvm.i'
! compressible version of SGS model
INTEGER :: lessvm_version = LESSVM_COMPRESSIBLE
! alignment with vorticity component allows backscatter
INTEGER :: lessvm_align_mode = LESSVM_ANSATZ2
INTEGER :: lessvm_bscatter = LESSVM_BACKSCATTER_ON
INTEGER :: lessvm_cutoff_safe = LESSVM_CUTOFF_SAFE_ON
! cutoff in terms of velocity difference sonsidered too small
DOUBLE PRECISION, PARAMETER :: cutoff = 1.0d-18
! structure function is velocity squared (so we need square cutoff)
DOUBLE PRECISION, PARAMETER :: cutoff2 = 1.0d-36
DOUBLE PRECISION, PARAMETER, PRIVATE :: pi = &
3.141592653589793238462643383279502884197d0
! This arrays are temporarily disabled since they are not really usede
! anyware. They can be reactivated once the models are validated.
! SGS tau_ij tensor
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgs11,sgs12,sgs13
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgs22,sgs23
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgs33
! General passive scalar transport symetric tensor
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive11
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive12
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive13
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive22
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive23
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgspassive33
! Structure functions
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: fstruc
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:) :: sgstemp
CONTAINS
!-----------------------------------------------------------
!
! subroutine: SpiralModel
! created:
! author: D.J. Hill
! C. Pantano
! description: actual sgs model
!
! Copyright (C) California Institute of Technology 2000-2005
!
! ---- requires data on the domain
! ---- [ixlo,ixhi] x [iylo,iyhi] x [izlo,izhi]
! ---- Calculates model terms on the domain
! ---- [ixlo+1,ixhi-1] x [iylo+1,iyhi-1] x [izlo+1,izhi-1]
! ----- This routine computes the subgrid model terms
! ----- including the subgrid scalar variance of the flow.
! ----- The energy spectrum in each cell E(k) is assumed to
! ----- result from two T.S.Lundgren style subgrid spiral vortices:
! ----- one aligned with the background strain, the other with the
! ----- background vorticity.
! ----- E(k) =K_0\epsilon^{3/2} ( lam E_strn(k) + (1-lam)E_vor(k) )
! ----- where E_{strn||vort)(k) k^{-5/3}exp(-2/3 \nu k^2 /a{s||v})
! ----- using the alingment model
! ----- lam P(strain) + (1-lam) P(vorticity)
! ----- where lam = major_eigen/( major_eigen + || w || )
! ----- The energy in a subgrid structure is given by
! ----- Int(E(k),k=kc..infty) = K_0\epsilon^{3/2} \Gamma[a,Y]
! ----- where a = -1/3. and Y = 2 k_c^2 nu /( 3 |a|)
! ----- and \Gamma[a,y] = Int( t^(a-1) exp(-t),t=y..inifnity)
! -----
! ----- the subgrid of the scalar spectrum Es(k) is assumed to be of
! ----- the Obukov-Corrsin form (inertial-convective subrange)
! ----- modified to take into account the vortex alignment.
! ----- Es_{sgs}(k) = \beta\epsilon^{-1/3}\epsilon^{c}(
! lam Es_strn(k) + (1-lam)Es_vor(k))
! ----- where \epsilon is the local energy dissipation
! ----- and \epsilon^{c} is the local scalar dissipation
! ----- and \beta is the Obukov-Corrsin prefactor
! ----- and Es_{strn||vort)(k) k^{-5/3}exp(-2/3 (2\nu+D) k^2 /a{s||v})
! ----- and D is the diffusivity (assumed to be zero)
! ----- On return: ux(nvars + 2 + is, :,:,:) holds the subgrid scalar
! ----- variance for the is passive scalar ( found in
! ----- vx(5+is,:,:,:) )
! ---------------------------------------------------------------------
SUBROUTINE SpiralModel(ux,vx)
use mesh
use array_bounds
use method_parms
USE Generic_Transport
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION, INTENT(INOUT):: ux(ncomps,ixlo:ixhi,iylo:iyhi,izlo:izhi)
DOUBLE PRECISION, INTENT(IN):: vx(nvars,ixlo:ixhi,iylo:iyhi,izlo:izhi)
INTEGER :: i,j,k
DOUBLE PRECISION, DIMENSION(2,3,3) :: etze
DOUBLE PRECISION, DIMENSION(3,3) :: vgrad
DOUBLE PRECISION :: sij11,sij12,sij13,sij22,sij23,sij33
DOUBLE PRECISION :: eval, nu,scalarnu,ke_strn,ke_vort
DOUBLE PRECISION :: e3x, e3y, e3z, lam, oabs, ox, oy, oz, k_c
DOUBLE PRECISION :: StrainWeight, VorticityWeight
DOUBLE PRECISION :: StructureAveStrain, StructureAveVorticity
DOUBLE PRECISION :: weighted_alignments, dx2, dy2, dz2
DOUBLE PRECISION :: strnweighted_ke, vortweighted_ke
DOUBLE PRECISION :: variancevort,variancestrn,strucscal
DOUBLE PRECISION :: ratiox,ratioy,ratioz,rootStrn,rootVort,fluxf
DOUBLE PRECISION :: del,localdel,s3prm,fact,facta,s3hat
DOUBLE PRECISION :: s3prmv, s3prms, rho, minsize, norm
DOUBLE PRECISION :: k0grp,c13,c23,sk2
DOUBLE PRECISION :: cutoffv, epssgs
! external modules
interface
Subroutine IntegrateStructure(ex,ey,ez,rx,ry,rz,StructureAve)
double precision ex, ey, ez
double precision rx, ry, rz
double precision StructureAve
end Subroutine IntegrateStructure
end interface
call cleslog_log_enter('SpiralModel')
! avoid calling this function before initializing the patch
if ( dx .eq. 0.0d0 .or. dy .eq. 0.0d0 .or. dz .eq. 0.0d0 ) goto 10
! --- a regularization parameter that represents the smallest
! --- possible strain rate : in nondim form:
! --- a/(kc kc nu) = max(abs(a)/(kc*kc*nu),minsize)
! --- :vaule 1/600 results in (del/pi)^(2/3) * 4.4d-264
! --- for the energy integrals.
minsize = 1.d0/600.d0
! length scales: strucr defined as in StructureFunction() !!
! --- the finite seperation in the structure function
c13 = (1.0d0/3.0d0)
c23 = (2.0d0/3.0d0)
localdel = (dx*dy*dz)**c13
! --- establish local aspect ratio since the structure function is
ratiox = dsqrt(dy*dz)/dx ! --- evaluated in the yz-plane
ratioy = dsqrt(dx*dz)/dy ! --- evaluated in the xz-plane
ratioz = dsqrt(dx*dy)/dz ! --- evaluated in the xy-plane
! ---- initialize the (delta_ij - e_i e_j) tensor
! ---- extended to hold two versions...
! ---- etze(1,i,j) is for e_i, e_j from resolved rate of strain
! ---- etze(2,i,j) is for e_i, e_j from resolved vorticity
etze(:,:,:) = 0.d0
! Initialize sgske locations. This is important because of the
! multiple cycle statements that can bypass setting the correct values.
if ( lessvm_version .eq. LESSVM_COMPRESSIBLE ) then
ux(nvars+3,:,:,:) = 0.0d0
endif
sgs11 = 0.0d0
sgs22 = 0.0d0
sgs33 = 0.0d0
sgs12 = 0.0d0
sgs13 = 0.0d0
sgs23 = 0.0d0
sgspassive11 = 0.0d0
sgspassive22 = 0.0d0
sgspassive33 = 0.0d0
sgspassive12 = 0.0d0
sgspassive13 = 0.0d0
sgspassive23 = 0.0d0
dx2 = 0.5d0/dx
dy2 = 0.5d0/dy
dz2 = 0.5d0/dz
! vorticity is velocity over lenghtscale (so we need gradient cutoff)
cutoffv = cutoff/localdel
DO k=izlo+1,izhi-1
DO j=iylo+1,iyhi-1
DO i=ixlo+1,ixhi-1
! -----------------------------------------------------------
! ----- COMPUTE RESOLVED SCALE QUANTITIES TO BE USED IN MODEL
! ---- cycle if the sturcture function is too small
if (fstruc(i,j,k,1) .LT. cutoff2) then
sgstemp(i,j,k) = 0.0d0
cycle
endif
! ---- get local velocity gradient
! ---- vgrad(i,j) = du_i/dx_j
vgrad(1,1) = (vx(2,i+1,j,k) - vx(2,i-1,j,k))*dx2
vgrad(1,2) = (vx(2,i,j+1,k) - vx(2,i,j-1,k))*dy2
vgrad(1,3) = (vx(2,i,j,k+1) - vx(2,i,j,k-1))*dz2
vgrad(2,1) = (vx(3,i+1,j,k) - vx(3,i-1,j,k))*dx2
vgrad(2,2) = (vx(3,i,j+1,k) - vx(3,i,j-1,k))*dy2
vgrad(2,3) = (vx(3,i,j,k+1) - vx(3,i,j,k-1))*dz2
vgrad(3,1) = (vx(4,i+1,j,k) - vx(4,i-1,j,k))*dx2
vgrad(3,2) = (vx(4,i,j+1,k) - vx(4,i,j-1,k))*dy2
vgrad(3,3) = (vx(4,i,j,k+1) - vx(4,i,j,k-1))*dz2
! ---- resolved rate of strain tensor
! ---- S_ij = (du_i/dx_j + du_j/dx_i)/2
sij11 = vgrad(1,1)
sij12 = 0.5d0*(vgrad(1,2)+vgrad(2,1))
sij13 = 0.5d0*(vgrad(1,3)+vgrad(3,1))
sij22 = vgrad(2,2)
sij23 = 0.5d0*(vgrad(2,3)+vgrad(3,2))
sij33 = vgrad(3,3)
if ( lessvm_align_mode .eq. LESSVM_DYNAMIC ) then
e3x = vx(nvars-2,i,j,k)
e3y = vx(nvars-1,i,j,k)
e3z = vx(nvars ,i,j,k)
! to be absolutely sure, normalize this vector
norm = dsqrt(e3x*e3x+e3y*e3y+e3z*e3z)
if ( norm .gt. 0.0d0 ) then
e3x = e3x/norm
e3y = e3y/norm
e3z = e3z/norm
else
sgstemp(i,j,k) = 0.0d0
cycle
endif
else
! ---- compute the principle eigen vector/value for strain S_ij
call spiral_eigen(sij11,sij12,sij13,sij22,sij23,sij33, &
e3x,e3y,e3z,eval)
! --- we couldn't find a local strain direction
if (eval.eq.0.0d0) then
sgstemp(i,j,k) = 0.0d0
cycle
endif
eval = ABS(eval)
endif
! --- make (delta_ij - e_i e_j) based on principle rate of strain
etze(1,1,1) = 1.d0 - e3x*e3x
etze(1,1,2) = - e3x*e3y
etze(1,1,3) = - e3x*e3z
etze(1,2,2) = 1.d0 - e3y*e3y
etze(1,2,3) = - e3y*e3z
etze(1,3,3) = 1.d0 - e3z*e3z
etze(1,2,1) = etze(1,1,2)
etze(1,3,1) = etze(1,1,3)
etze(1,3,2) = etze(1,2,3)
! integral from structure function equation
! the type of averaging (circular or spherical)
! is defined in the /equations/structurefunctions.f file
CALL IntegrateStructure(e3x,e3y,e3z,ratiox,ratioy,ratioz, &
StructureAveStrain)
IF( lessvm_align_mode .eq. LESSVM_ANSATZ2 ) THEN
! contribution from vortex alignment with vorticity
! (only for non-zero vorticity)
! ---- construct the local resolved vorticity
ox = vgrad(3,2) - vgrad(2,3)
oy = vgrad(1,3) - vgrad(3,1)
oz = vgrad(2,1) - vgrad(1,2)
oabs= sqrt(ox*ox + oy*oy + oz*oz)
lam = eval + oabs
! --- another test to make sure the flow wasn't laminar here
if ( lam .le. cutoffv ) then
sgstemp(i,j,k) = 0.0d0
cycle
endif
! --- constructing the alignment factors
StrainWeight = eval/lam
VorticityWeight = oabs/lam
! --- normalize the vorticity
if ( oabs .gt. cutoffv ) then
ox = ox/oabs
oy = oy/oabs
oz = oz/oabs
! --- make (delta_ij - e_i e_j) bases on resolved vorcity
etze(2,1,1) = 1.d0 - ox*ox
etze(2,1,2) = - ox*oy
etze(2,1,3) = - ox*oz
etze(2,2,2) = 1.d0 - oy*oy
etze(2,2,3) = - oy*oz
etze(2,3,3) = 1.d0 - oz*oz
etze(2,2,1) = etze(2,1,2)
etze(2,3,1) = etze(2,1,3)
etze(2,3,2) = etze(2,2,3)
! integral for aligment with vorticity
! the type of averaging (circular or spherical)
! is defined in the /equations/structurefunctions.f file
CALL IntegrateStructure(ox,oy,oz,ratiox,ratioy,ratioz, &
StructureAveVorticity)
else
etze(2,:,:) = 0.0d0
StrainWeight = 1.d0
VorticityWeight = 0.0d0
StructureAveVorticity = 0.d0
endif
ELSE
StrainWeight = 1.d0
VorticityWeight = 0.0d0
StructureAveVorticity = 0.d0
ENDIF
weighted_alignments = (StrainWeight*StructureAveStrain + &
VorticityWeight*StructureAveVorticity)
! Ko\epsilon^(2/3) computed by matching the resolved scale
strucscal = 1.0d0/((localdel**c23) * weighted_alignments)
k0grp = fstruc(i,j,k,1)*strucscal
! -------------------------------------------------------
! ----- Now we compute the contribution from a subgrid
! ----- Lundgrend type vortex, need kinematic VISCOSITY ->
! ----- divide by \rho for run with zero viscosity, use
! ----- nominal value for model nu = 2.0/1.d4/vx(1,i,j,k)
rho = vx(1,i,j,k)
if ( trsp_visc_mode .eq. CLES_TRANSP_VISC_VAR ) then
nu = mu(i,j,k)/rho
else
nu = mu(1,1,1)/rho
endif
! ---- in the scalar model we assume
! ---- k^(-5/3)exp(-2/3 (2nu+D)/a k^2)
! ---- call 2nu+d the scalarnu. Assume d=nu for now
scalarnu = 2.0d0*nu + nu
! --- establish the lenght scale for the SGS terms in physical
! --- space and wavenumber space
if ( lessvm_cutoff_safe .eq. LESSVM_CUTOFF_SAFE_ON ) then
del = max(sgstemp(i,j,k), localdel)
else
del = sgstemp(i,j,k)
endif
k_c = pi/del
sk2 = 2.0d0*k_c*k_c/3.0d0
facta =.5d0*(1.d0/k_c)**c23
! ------- Alignment with Resolved scale Strain:
! ---- Get the axial rate of strain along the subgrind vortex axis
! ---- computed as (e3)^t [S_ij] (e3)
! ---- this should just be the eigen value corresponding
! ---- to the eigen vector e3 of S_ij
s3prm = sij11*e3x*e3x + sij22*e3y*e3y + sij33*e3z*e3z &
& + 2.0d0*(sij12*e3x*e3y+sij13*e3x*e3z+sij23*e3y*e3z)
! ---- Integrate the LundgrenSpiral to get K=Int(E(k),kc..inf)
! ----- The energy in a subgrid structure is given by
! ----- Int(E(k),k=kc..infty) =
! K_0\epsilon^{3/2}*.5*(1/kc)^(2/3)Y^(1/3)\Gamma[a,Y]
! ----- where a = -1/3. and Y = 2 k_c^2 nu /( 3 |a|)
! ----- and \Gamma[a,y] = Int( t^(a-1) exp(-t),t=y..inifnity)
! make sure the argument to the integral isn't too big
! for our simulation, - as would result from a zero strain
s3prmv = sk2*nu*max(abs(s3prm)/(sk2*nu),minsize)
s3hat = sk2*nu/s3prmv
fact = facta* s3hat**c13
ke_strn = k0grp*fact*gaminc(s3hat)
strnweighted_ke = StrainWeight*ke_strn
! ---- For the passive scalar variance - do same integral
! ---- but with a different 'viscosity'
! ---- Integrate the LundgrenSpiral to get K=Int(Es(k),kc..inf)
s3prms = sk2*scalarnu*max(abs(s3prm)/(sk2*scalarnu),minsize)
s3hat = sk2*scalarnu/s3prms
fact = facta* (s3hat)**c13
variancestrn = fact*gaminc(s3hat)
variancestrn = StrainWeight*variancestrn
! ------- Alignment with Resolved scale Vorticity
! ---- Get the axial rate of strain along the subgrind vortex axis
! ---- computed as (omega)^t [S_ij] (omega)
! ---- notic, this doesn't have to be an eigen value of S_ij
IF( lessvm_align_mode .eq. LESSVM_ANSATZ2 ) THEN
s3prm = sij11*ox*ox + sij22*oy*oy+ sij33*oz*oz &
& + 2.0*(sij12*ox*oy+ sij13*ox*oz+ sij23*oy*oz )
! ---- need to desingularize
! ---- because we can have the vorticity
! ---- in the null space for the strain tensor
s3prmv = sk2*nu*max(abs(s3prm)/(sk2*nu),minsize)
! ---- Integrate the LundgrenSpiral to get K=Int(E(k),kc..inf)
! ----- The energy in a subgrid structure is given by
! ----- Int(E(k),k=kc..infty) =
! K_0\epsilon^{3/2}*.5*(1/kc)^(2/3)Y^(1/3)\Gamma[a,Y]
! ----- where a = -1/3. and Y = 2 k_c^2 nu /( 3 |a|)
! ----- and \Gamma[a,y] = Int( t^(a-1) exp(-t),t=y..inifnity)
s3hat =sk2*nu/s3prmv
fact = facta* (s3hat)**c13
ke_vort = k0grp*fact*gaminc(s3hat)
vortweighted_ke = VorticityWeight*ke_vort
! ---- For the passive scalar variance - do same integral
! ---- but with a different 'viscosity'
! ---- Integrate the LundgrenSpiral to get K=Int(Es(k),kc..inf)
s3prms = sk2*scalarnu*max(abs(s3prm)/(sk2*scalarnu),minsize)
s3hat = sk2*scalarnu/s3prms
fact = facta* (s3hat)**c13
variancevort = fact*gaminc(s3hat)
variancevort = VorticityWeight*variancevort
ELSE
ke_vort = 0.0d0
variancevort = 0.0d0
vortweighted_ke = 0.0d0
ENDIF
strnweighted_ke = strnweighted_ke* rho
vortweighted_ke = vortweighted_ke* rho
! calculate SGS stresses
sgs11(i,j,k) = strnweighted_ke*etze(1,1,1)
sgs22(i,j,k) = strnweighted_ke*etze(1,2,2)
sgs33(i,j,k) = strnweighted_ke*etze(1,3,3)
sgs12(i,j,k) = strnweighted_ke*etze(1,1,2)
sgs13(i,j,k) = strnweighted_ke*etze(1,1,3)
sgs23(i,j,k) = strnweighted_ke*etze(1,2,3)
if ( lessvm_align_mode .eq. LESSVM_ANSATZ2 ) then
sgs11(i,j,k) = sgs11(i,j,k)+vortweighted_ke*etze(2,1,1)
sgs22(i,j,k) = sgs22(i,j,k)+vortweighted_ke*etze(2,2,2)
sgs33(i,j,k) = sgs33(i,j,k)+vortweighted_ke*etze(2,3,3)
sgs12(i,j,k) = sgs12(i,j,k)+vortweighted_ke*etze(2,1,2)
sgs13(i,j,k) = sgs13(i,j,k)+vortweighted_ke*etze(2,1,3)
sgs23(i,j,k) = sgs23(i,j,k)+vortweighted_ke*etze(2,2,3)
endif
if ( lessvm_bscatter .eq. LESSVM_BACKSCATTER_OFF ) then
epssgs = -sij11*sgs11(i,j,k) - sij22*sgs22(i,j,k) &
- sij33*sgs33(i,j,k) - 2.0d0*(sij12*sgs12(i,j,k)+ &
sij13*sgs13(i,j,k)+sij23*sgs23(i,j,k))
if ( epssgs .lt. 0.0d0 ) then
sgs11(i,j,k) = 0.0d0
sgs22(i,j,k) = 0.0d0
sgs33(i,j,k) = 0.0d0
sgs12(i,j,k) = 0.0d0
sgs13(i,j,k) = 0.0d0
sgs23(i,j,k) = 0.0d0
endif
endif
!----- subgrid TKE for statistics and pressure correction
!----- 1/2 (Tau_ii)
if ( lessvm_version .eq. LESSVM_COMPRESSIBLE ) then
ux(nvars+3,i,j,k) = 0.5d0*(sgs11(i,j,k)+sgs22(i,j,k)&
+sgs33(i,j,k))
endif
!---- compute subgrid temperature fluxes
rootStrn = StrainWeight*dsqrt(ke_strn)
rootVort = VorticityWeight*dsqrt(ke_vort)
fluxf = 0.5*pi/k_c
! ---- create the _symetric_ SGS-passive transport tensor
! ---- this is used later for the Tempurature and for all
! ---- the passive scalars.
! ----: PassiveSubgridFlux_i = -sgspassive_ij (d Passive/dx_j)
! ---- for instance the tempurature flux is given by
! rho q_i = -rho * sgspassive_ij (d T/dx_j)
! -rho q_i gamma/(gamma-1)
rootStrn = rootStrn * fluxf
rootVort = rootVort * fluxf
sgspassive11(i,j,k) = rootStrn*etze(1,1,1)
sgspassive12(i,j,k) = rootStrn*etze(1,1,2)
sgspassive13(i,j,k) = rootStrn*etze(1,1,3)
sgspassive22(i,j,k) = rootStrn*etze(1,2,2)
sgspassive23(i,j,k) = rootStrn*etze(1,2,3)
sgspassive33(i,j,k) = rootStrn*etze(1,3,3)
if ( lessvm_align_mode .eq. LESSVM_ANSATZ2 ) then
sgspassive11(i,j,k) = sgspassive11(i,j,k)+rootVort*etze(2,1,1)
sgspassive12(i,j,k) = sgspassive12(i,j,k)+rootVort*etze(2,1,2)
sgspassive13(i,j,k) = sgspassive13(i,j,k)+rootVort*etze(2,1,3)
sgspassive22(i,j,k) = sgspassive22(i,j,k)+rootVort*etze(2,2,2)
sgspassive23(i,j,k) = sgspassive23(i,j,k)+rootVort*etze(2,2,3)
sgspassive33(i,j,k) = sgspassive33(i,j,k)+rootVort*etze(2,3,3)
endif
! Scalar factor
! beta\epsilon^(-1/3)\epsilon_o :
! matching the resolved scale
! ---- compute Int(E(k),k=kc..infinity)
! ----
! ---- where Es(k)
! = Betagrp*(\lamda Es(k,vortex aligned) +
! (1-\lamda) Es(k,strain aligned) )
! ---- put the scalar variance in the structure function
! ---- and multiply by 2.0d0 since the structure function
! ---- matching implicitly assumes Int(E) = 1/2 |s|^2
sgstemp(i,j,k) = 2.0d0*(variancevort+variancestrn)*strucscal
END DO
END DO
END DO
10 call cleslog_log_exit('SpiralModel')
RETURN
END SUBROUTINE SpiralModel
!---------------------------------------------------------------------------
!
! AUXILIARY FUNCTIONS: gaminc, eigsrt, eigen, evector
!
!---------------------------------------------------------------------------
!-----------------------------------------------------------
!
! subroutine: gaminc
! created:
! author: D.I. Pullin
! D.J. Hill
! description: computes the incomplete gamma function
! \Gamma[a,y] = Int( t^(a-1) exp(-t),t=y..inifnity)
! this integral assumes a = -1/3.
!
! Copyright (C) California Institute of Technology 2000-2005
!
!-----------------------------------------------------------
DOUBLE PRECISION FUNCTION gaminc(y)
IMPLICIT NONE
DOUBLE PRECISION y,ythd,ysq,ycub,denom,num
DOUBLE PRECISION g(42),p,pm1,pp1,pm2
INTEGER j
DATA g/0.,0.,1.053318515,&
& 0.7635093981, 0.576870, 0.4479712295,0.3547708363,&
& 0.285157594, 0.2318870423, 0.1903499802,0.1574731302,&
& 0.1311305945,0.1098077007,0.09239878998,0.07808029672,&
& 0.06622821442,0.05636279002,0.0481104998,0.04117731868,&
& 0.03532956022,0.03037990663,0.02617706753,0.02259802105,&
& 0.01954212176,0.0169265774,0.01468294222,0.01275437411,&
& 0.01109347179,0.009660556665,0.00842229874,0.007350610742,&
& 0.006421753093,0.005615605641,0.004915072194,0.004305591432,&
& 0.003774733502,0.003311865981,0.002907876269,0.002554940089,&
& 0.002246327831,0.001976242056,0.001739680777/
! if(0.5= tol) THEN
eta = min(max(r/(q*sqrt(q)),-1.0d0),1.0d0)
b = ACOS(eta)
eval1 = -2.0*SQRT(q)*COS(b/3.0) - a/3.0
eval2 = -2.0*SQRT(q)*COS((b + 2.0*pi)/3.0) - a/3.0
eval3 = -2.0*SQRT(q)*COS((b - 2.0*pi)/3.0) - a/3.0
ELSE
eval1 = a11
eval2 = a22
eval3 = a33
END IF
! Sort eigenvalues: eval1 < eval2 < eval3
IF (eval2 < eval1) THEN
b = eval1
eval1 = eval2
eval2 = b
END IF
IF (eval3 < eval2) THEN
b = eval2
eval2 = eval3
eval3 = b
END IF
IF (eval2 < eval1) THEN
b = eval1
eval1 = eval2
eval2 = b
END IF
! Compute eigenvectors
CALL spiral_evector(eval1,a11,a12,a13,a22,a23,a33,1,evec1)
CALL spiral_evector(eval2,a11,a12,a13,a22,a23,a33,2,evec2)
CALL spiral_evector(eval3,a11,a12,a13,a22,a23,a33,3,evec3)
! copy to array for better sorting
evals(1) = eval1
evals(2) = eval2
evals(3) = eval3
do i = 1, 3
evecs(i,1) = evec1(i)
evecs(i,2) = evec2(i)
evecs(i,3) = evec3(i)
enddo
! check for non-convergence
IF((evals(1) .EQ. 0.) .AND. (evals(2) .EQ. 0.) .AND. &
(evals(3) .EQ. 0.)) THEN
evecs = 0.0d0
evecs(1,1)=1.0d0
evecs(2,2)=1.0d0
evecs(3,3)=1.0d0
ELSE
CALL spiral_eigsrt(evals, evecs, 3, 3)
END IF
! ---- the principal eigenvector
e3x = evecs(1,1)
e3y = evecs(2,1)
e3z = evecs(3,1)
eval = evals(1)
RETURN
END SUBROUTINE spiral_eigen
!-----------------------------------------------------------
!
! subroutine: spiral_eigsrt
! created:
! author: D.I. Pullin
! D.J. Hill
! description: sort eigensystem of resolved
! rate of strain tensor
!
! Copyright (C) California Institute of Technology 2000-2005
!
!-----------------------------------------------------------
SUBROUTINE spiral_eigsrt(d,v,n,np)
IMPLICIT NONE
INTEGER n,np
DOUBLE PRECISION d(np),v(np,np)
INTEGER i,j,k
DOUBLE PRECISION p
DO i=1,n-1
k=i
p=d(i)
DO j=i+1,n
IF(d(j).GE.p)THEN
k=j
p=d(j)
ENDIF
END DO
IF(k.NE.i)THEN
d(k)=d(i)
d(i)=p
DO j=1,n
p=v(j,i)
v(j,i)=v(j,k)
v(j,k)=p
END DO
ENDIF
END DO
END SUBROUTINE spiral_eigsrt
!-----------------------------------------------------------
!
! subroutine: spiral_evector
! created:
! author: D.I. Pullin
! D.J. Hill
! description: compute normalized eigenvectors of resolved
! rate of strain tensor
!
! Copyright (C) California Institute of Technology 2000-2005
!
!-----------------------------------------------------------
SUBROUTINE spiral_evector(eval,a11,a12,a13,a22,a23,a33,iunit,evec)
! Assumes eigenvalues are DOUBLE PRECISION and distinct, unless eval == 0
IMPLICIT NONE
DOUBLE PRECISION, INTENT(IN) :: a11,a12,a13,a22,a23,a33
DOUBLE PRECISION, INTENT(IN) :: eval
DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: evec
INTEGER, INTENT(IN) :: iunit
DOUBLE PRECISION ::det1,det2,det3,adet1,adet2,adet3
DOUBLE PRECISION, PARAMETER :: tol = 1.0d-15
det1 = (a22 - eval)*(a33 - eval) - a23**2
det2 = (a11 - eval)*(a33 - eval) - a13**2
det3 = (a11 - eval)*(a22 - eval) - a12**2
adet1 = ABS(det1)
adet2 = ABS(det2)
adet3 = ABS(det3)
IF (adet1 == MAX(adet1,adet2,adet3).AND.adet1>=tol)THEN
evec(1) = 1.0
evec(2) = (-a12*(a33 - eval) + a13*a23)/det1
evec(3) = (-a13*(a22 - eval) + a12*a23)/det1
END IF
IF(adet2 == MAX(adet1,adet2,adet3).AND.adet2 >= tol)THEN
evec(1) = (-a12*(a33 - eval) + a23*a13)/det2
evec(2) = 1.0
evec(3) = (-a23*(a11 - eval) + a12*a13)/det2
END IF
IF(adet3 == MAX(adet1,adet2,adet3).AND.adet3 >= tol)THEN
evec(1) = (-a13*(a22 - eval) + a23*a12)/det3
evec(2) = (-a23*(a11 - eval) + a13*a12)/det3
evec(3) = 1.0
END IF
IF(adet1