! ==========================================================
SUBROUTINE stdcdat()
use method_parms
! ==========================================================
implicit none
include 'cles.i'
integer :: lmechout, maxsize
parameter( lmechout = 11, maxsize=256 )
double precision :: Wk(maxsize)
double precision :: RU, PA
character(len=16) :: cwork(maxsize)
character(len=16) :: stmp
integer :: k, i, nmax, ip
if ( ncomps .gt. maxsize ) then
print *, 'Error: we did not expect such large vector of state'
stop
endif
! scalars
do i=1, nscal
write(stmp,398) i
cwork(i) = stmp
enddo
ip = nscal+1
! dcflag
cwork(ip)= 'DCFlag'
ip = ip + 1
! if using LES, have access to all other fields
if ( useLES .eq. CLES_TRUE ) then
! sgs kinetic energy
cwork(ip)= 'SGSKinEnergy'
ip = ip + 1
! sgs component 11
cwork(ip)= 'sgs11'
ip = ip + 1
! sgs component 22
cwork(ip)= 'sgs22'
ip = ip + 1
! sgs component 33
cwork(ip)= 'sgs33'
ip = ip + 1
! sgs component 12
cwork(ip)= 'sgs12'
ip = ip + 1
! sgs component 13
cwork(ip)= 'sgs13'
ip = ip + 1
! sgs component 23
cwork(ip)= 'sgs23'
ip = ip + 1
! velocity structure function
cwork(ip)= 'SGSDissipation'
ip = ip + 1
! scalar variance and structure functions
do i=1,nscal
write(stmp,399) i
cwork(ip) = stmp
ip = ip + 1
enddo
do i=1, nscal+1
write(stmp,404) i
cwork(ip) = stmp
ip = ip + 1
enddo
endif
nmax = ip-1
do i=1, nmax
Wk(i) = 1.d0
enddo
RU = 1.d0
PA = 1.d0
open(unit=lmechout, status='unknown', form='formatted', &
& file='chem.dat')
write (lmechout,400) RU
write (lmechout,401) PA
write (lmechout,402) (k,cwork(k),k=1,nmax)
write (lmechout,403) (k,Wk(k),k=1,nmax)
close (lmechout)
398 format('Scalar',i2.2)
399 format('ScalarVar',i2.2)
404 format('Structure',i2.2)
400 format('RU ',e16.8)
401 format('PA ',e16.8)
402 format('Sp(',i2.2,') ',a16)
403 format('W(',i2.2,') ',e16.8)
RETURN
END SUBROUTINE stdcdat
! ==========================================================
SUBROUTINE out1eu(q,mx_,lb,ub,qo,mxo,lbo,ubo,&
& lbr,ubr,shaper,meqn,nc,time)
! ==========================================================
! ---- share variables
USE method_parms
use cles_interfaces
! ---- share procedures
USE Generic_EvalGamma
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION :: vx(nvars)
INTEGER meqn, mx_, mxo, nc
DOUBLE PRECISION q(meqn,mx_), qo(mxo), time
INTEGER lb(1), ub(1), lbo(1), ubo(1), lbr(1), ubr(1), shaper(1),&
& stride, imin(1), imax(1), i, getindx
stride = (ub(1) - lb(1))/(mx_-1)
imin(1) = MAX(lb(1), lbo(1), lbr(1))
imax(1) = MIN(ub(1), ubo(1), ubr(1))
IF (MOD(imin(1)-lb(1),stride) .NE. 0) THEN
imin(1) = imin(1) + stride - MOD(imin(1)-lb(1),stride)
ENDIF
imin(1) = getindx(imin(1), lb(1), stride)
IF (MOD(imax(1)-lb(1),stride) .NE. 0) THEN
imax(1) = imax(1) - MOD(imax(1)-lb(1),stride)
ENDIF
imax(1) = getindx(imax(1), lb(1), stride)
DO i = imin(1), imax(1)
if ( nc .ge. 4 .and. nc .le. 6 ) then
vx(1) = q(1,i)
vx(2:nvars) = q(2:nvars,i)/vx(1)
call cles_eqstate(q(1,i), ncomps, vx, nvars, 1, CLES_FALSE)
endif
IF (nc.EQ.1) qo(i) = q(1,i)
IF (nc.EQ.2) qo(i) = q(2,i)/q(1,i)
IF (nc.EQ.3) qo(i) = q(5,i)
IF (nc.EQ.4) qo(i) = q(nvars+1,i)
IF (nc.EQ.5) qo(i) = vx(5)
IF (nc.EQ.6) CALL GetGamma(q(:,i), vx, qo(i))
IF (nc.ge.(nvars-nscal+2) .and. nc.le.(nvars+1)) qo(i) = q(nc-1,i)/q(1,i)
IF (nc.EQ.(nvars+2)) qo(i) = q(nvars+2,i)
if (nc.gt.(nvars+2) ) qo(i) = q(nc,i)
END DO
! set extendend output to false so we force a call to set array bounds
useExOutput = CLES_FALSE
RETURN
END SUBROUTINE out1eu
! ==========================================================
SUBROUTINE out2eu(q,mx_,my_,lb,ub,qo,mxo,myo,lbo,ubo,&
& lbr,ubr,shaper,meqn,nc,time)
! ==========================================================
! ---- share variables
USE method_parms
use cles_interfaces
! ---- share procedures
USE Generic_EvalGamma
IMPLICIT NONE
include 'cles.i'
DOUBLE PRECISION :: vx(nvars)
INTEGER meqn, mx_, my_, mxo, myo, nc
DOUBLE PRECISION q(meqn,mx_,my_), qo(mxo,myo), time
INTEGER lb(2), ub(2), lbo(2), ubo(2), lbr(2), ubr(2), shaper(2),&
& stride, imin(2), imax(2), i, j, getindx, d
stride = (ub(1) - lb(1))/(mx_-1)
DO d = 1, 2
imin(d) = MAX(lb(d), lbr(d))
imax(d) = MIN(ub(d), ubr(d))
IF (MOD(imin(d)-lb(d),stride) .NE. 0) THEN
imin(d) = imin(d) + stride - MOD(imin(d)-lb(d),stride)
ENDIF
imin(d) = getindx(imin(d), lb(d), stride)
IF (MOD(imax(d)-lb(d),stride) .NE. 0) THEN
imax(d) = imax(d) - MOD(imax(d)-lb(d),stride)
ENDIF
imax(d) = getindx(imax(d), lb(d), stride)
END DO
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
if ( nc .ge. 5 .and. nc .le. 7 ) then
vx(1) = q(1,i,j)
vx(2:nvars) = q(2:nvars,i,j)/vx(1)
call cles_eqstate(q(1,i,j), ncomps, vx, nvars, 1, CLES_FALSE)
endif
IF (nc.EQ.1) qo(i,j) = q(1,i,j)
IF (nc.EQ.2) qo(i,j) = q(2,i,j)/q(1,i,j)
IF (nc.EQ.3) qo(i,j) = q(3,i,j)/q(1,i,j)
IF (nc.EQ.4) qo(i,j) = q(5,i,j)
IF (nc.EQ.5) qo(i,j) = q(nvars+1,i,j)
IF (nc.EQ.6) qo(i,j) = vx(5)
IF (nc.EQ.7) CALL GetGamma(q(:,i,j), vx, qo(i,j))
IF (nc.ge.(nvars-nscal+3) .and. nc.le.(nvars+2) ) &
qo(i,j) = q(nc-2,i,j)/q(1,i,j)
IF (nc.EQ.(nvars+3)) qo(i,j) = q(nvars+2,i,j)
if (nc.gt.(nvars+3) ) qo(i,j) = q(nc-1,i,j)
END DO
END DO
! set extendend output to false so we force a call to set array bounds
useExOutput = CLES_FALSE
RETURN
END SUBROUTINE out2eu
! ==========================================================
SUBROUTINE out3eu(q,mx_,my_,mz_,lb,ub,qo,mxo,myo,mzo,lbo,ubo,&
& lbr,ubr,shaper,meqn,nc,time)
! ==========================================================
!
! nc value | field
! 1 r
! 2 u
! 3 v
! 4 w
! 5 E
! 6 T
! 7 p
! 8 gamma
! 9 Y1
! 10 Y2
! ... ...
! 8+nscal Yn
! 9+nscal dcflag
! 10+nscal sgske
! 11+nscal sgs11
! 12+nscal sgs22
! 13+nscal sgs33
! 14+nscal sgs12
! 15+nscal sgs13
! 16+nscal sgs23
! 17+nscal sgseps
! 18+nscal sgsvar1
! 19+nscal sgsvar2
! ... ...
! 17+2*nscal sgsvarn
! 18+2*nscal fstruct
! 19+2*nscal sstruct1
! 20+2*nscal sstruct2
! ... ...
! 18+3*nscal sstructn
!
! ---- share variables
USE method_parms
USE array_bounds
use cles_interfaces
! ---- share procedures
USE Generic_EvalGamma
USE Generic_Transport
IMPLICIT NONE
include 'cles.i'
INTEGER meqn, mx_, my_, mz_, mxo, myo, mzo, nc
DOUBLE PRECISION time, q(meqn,mx_,my_,mz_), qo(mxo,myo,mzo)
INTEGER lb(3), ub(3), lbo(3), ubo(3), lbr(3), ubr(3), shaper(3),&
& stride, imin(3), imax(3), i, j, k, getindx, d
INTEGER :: nn, iret
EXTERNAL cles_hook_exist, cles_hook4
INTEGER :: has_hooks, cles_hook_exist
INTEGER :: ipar(4), npar, nq, nvx, tmp_alloc, ierr
DOUBLE PRECISION, ALLOCATABLE :: vx(:,:,:,:)
DOUBLE PRECISION, ALLOCATABLE :: tmp(:,:,:)
call cleslog_log_enter('out3eu')
stride = (ub(1) - lb(1))/(mx_-1)
DO d = 1, 3
imin(d) = MAX(lb(d), lbr(d))
imax(d) = MIN(ub(d), ubr(d))
IF (MOD(imin(d)-lb(d),stride) .NE. 0) THEN
imin(d) = imin(d) + stride - MOD(imin(d)-lb(d),stride)
ENDIF
imin(d) = getindx(imin(d), lb(d), stride)
IF (MOD(imax(d)-lb(d),stride) .NE. 0) THEN
imax(d) = imax(d) - MOD(imax(d)-lb(d),stride)
ENDIF
imax(d) = getindx(imax(d), lb(d), stride)
END DO
! density or total energy
if ( nc .eq. 1 .or. nc .eq. 5) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nc,i,j,k)
enddo
enddo
enddo
goto 10
endif
! temperature
if ( nc .eq. 6 ) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nvars+1,i,j,k)
enddo
enddo
enddo
goto 10
endif
! dcflag
IF (nc.EQ.(9+nscal)) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nvars+2,i,j,k)
enddo
enddo
enddo
goto 10
endif
! velocity components
if ( nc .ge. 2 .and. nc .le. 4 ) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nc,i,j,k)/q(1,i,j,k)
enddo
enddo
enddo
goto 10
endif
! scalars
IF (nc.ge.9 .and. nc.le.(8+nscal)) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nc-3,i,j,k)/q(1,i,j,k)
enddo
enddo
enddo
goto 10
endif
ALLOCATE( vx(nvars,mx_,my_,mz_), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'out3eu/vx', iret)
endif
! determine primitive variables
vx(1,:,:,:) = q(1,:,:,:)
do i=2, nvars
vx(i,:,:,:) = q(i,:,:,:)/q(1,:,:,:)
enddo
! if LES is active (we need to compute all this before
! computing p since it depends on the availability of sgske).
if ( useLES .eq. CLES_TRUE .and. useExOutput .eq. CLES_TRUE ) then
ALLOCATE( tmp(mx_,my_,mz_), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'out3eu/tmp', iret)
endif
call AllocateTransport()
call GetTransport(q, vx)
call AllocateLES()
call InitializeLES(q, vx)
! sgs stuff
if ( nc .eq. 11+nscal ) then
call LES_Spiral_Sgs(tmp, 1, 1)
else if ( nc .eq. 12+nscal ) then
call LES_Spiral_Sgs(tmp, 2, 2)
else if ( nc .eq. 13+nscal ) then
call LES_Spiral_Sgs(tmp, 3, 3)
else if ( nc .eq. 14+nscal ) then
call LES_Spiral_Sgs(tmp, 1, 2)
else if ( nc .eq. 15+nscal ) then
call LES_Spiral_Sgs(tmp, 1, 3)
else if ( nc .eq. 16+nscal ) then
call LES_Spiral_Sgs(tmp, 2, 3)
else if ( nc .eq. 17+nscal ) then
call LES_Spiral_Dissipation(vx, tmp)
else if ( nc .ge. 18+nscal .and. nc .le. 17+2*nscal ) then
call LES_Spiral_Variance(tmp, nc-17-nscal)
else if ( nc .ge. 18+2*nscal .and. nc .le. 18+3*nscal ) then
call LES_Spiral_Structure(tmp, nc-17-2*nscal)
else
goto 40
endif
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = tmp(i,j,k)
enddo
enddo
enddo
goto 30
endif
40 nn = mx_*my_*mz_
call cles_eqstate(q, ncomps, vx, nvars, nn, useLES)
IF (nc.EQ.7) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = vx(5,i,j,k)
END DO
END DO
END DO
else if ( nc .eq. 8 ) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
CALL GetGamma(q(:,i,j,k), vx(:,i,j,k), qo(i,j,k))
END DO
END DO
END DO
else if ( nc .eq. 10+nscal .and. useLES .eq. CLES_TRUE ) then
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = q(nvars+3,i,j,k)/q(1,i,j,k)
END DO
END DO
END DO
else if ( useExOutput .eq. CLES_TRUE ) then
! User provided hook
has_hooks = cles_hook_exist(CLES_HOOK_OUTPUT)
if ( has_hooks .eq. CLES_TRUE ) then
ipar(1) = nc
ipar(2) = mx_
ipar(3) = my_
ipar(4) = mz_
npar = 4
nq = ncomps*nn
nvx = nvars*nn
if ( .NOT. ALLOCATED(tmp) ) then
ALLOCATE( tmp(mx_,my_,mz_), stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_ALLOCATE, &
'out3eu/tmp/hook', iret)
endif
tmp_alloc = CLES_TRUE
else
tmp_alloc = CLES_FALSE
endif
call cles_hook4(CLES_HOOK_OUTPUT, ipar, npar, q, nq, vx, nvx, &
tmp, nn, ierr)
DO k = imin(3), imax(3)
DO j = imin(2), imax(2)
DO i = imin(1), imax(1)
qo(i,j,k) = tmp(i,j,k)
END DO
END DO
END DO
if ( tmp_alloc .eq. CLES_TRUE ) then
DEALLOCATE( tmp, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'out3eu/tmp/hook', iret)
endif
endif
endif
endif
30 if ( useLES .eq. CLES_TRUE .and. useExOutput .eq. CLES_TRUE ) then
DEALLOCATE( tmp, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'out3eu/tmp', iret)
endif
call DeallocateLES()
call DeallocateTransport()
endif
DEALLOCATE( vx, stat=iret )
if ( iret .ne. 0 ) then
call cleslog_error(CLESLOG_ERROR_DEALLOCATE, &
'out3eu/vx', iret)
endif
! set extended output to false so we force a call to set array bounds
10 useExOutput = CLES_FALSE
call cleslog_log_exit('out3eu')
RETURN
END SUBROUTINE out3eu