subroutine tkestatsdump
use modfields, only : u0,v0,w0,thl0,uav,vav,wav,uuav,vvav,wwav,uvav,uwav,vwav,thlav,thlthlav,pres0,thluav,thlvav,thlwav,&
upupav,vpvpav,wpwpav,thlpthlpav,upvpav,upwpav,vpwpav,thlpupav,thlpvpav,thlpwpav,presav,&
strain2av,disssgsav,t_vav,tvmx,tvmy,tvmz,tsgsmx1,tsgsmx2,tsgsmy1,tsgsmy2,tsgsmz1,t_sgsav,nusgsav,&
tpm,t_pav,ttmx,ttmy,ttmz,t_tav,p_bav,d_sgsav,p_tav,tkeadv,tsgsmz1,tsgsmz2,t_t,t_v,t_p,t_sgs,d_sgs,&
p_b,p_t,adv,IIc,IIcs
use modglobal, only : ib,ie,ih,jb,je,jgb,jge,dy,jh,ke,kb,kh,rk3step,timee,cexpnr,tsample,tstatsdump,jtot,imax,dzf,&
dzf,dzfi,dzhi,dxf,dxfi,dyi,dxhi,dy2i,grav,numol
use modmpi, only : myid,cmyid,my_real,mpi_sum,avey_ibm,mpierr,comm3d,excjs,avexy_ibm
use modsurfdata, only : thls
use modsubgrid, only : ekh
implicit none
real, dimension(ib:ie,jb:je,kb:ke) :: disssgsfl ! average subgrid visc. * average rate of strain squared : 2*<nu_t>*<Sij>*<Sij>
real, dimension(ib:ie,jb:je,kb:ke) :: dissresav ! average resolved dissipation: 2*nu*<Sij'*Sij'> = 2*nu*( <Sij*Sij> - <Sij>*<Sij> )
real, dimension(ib:ie,jb:je,kb:ke) :: tke ! tke = 0.5*<ui'ui'>
real, dimension(ib:ie,jb:je,kb:ke) :: mke ! = <ui>d/dxj(<ui><uj>) + <ui>d/dxj(<ui'uj'>) = <ui>d/dxj(<ui*uj>)
real, dimension(ib:ie+1,jb :je, kb:ke) :: dummyx
real, dimension(ib:ie, jb-1:je+1,kb:ke) :: dummyy
real, dimension(ib:ie, jb :je, kb:ke+1) :: dummyz
integer i,j,k,ip,im,jp,jm,kp,km
real strainav2
real dummy
! Tvav = (Tvm - <ui>*d/dxj(<Sij>) ) + 2*nu*<Sij'Sij'>
! Tvm = Tvmx + Tvmy + Tvmz -> therefore: subtraction, then interpolation,
! then addition of 2*nu*<Sij'Sij'>
do k=kb,ke
km = k-1
kp = k+1
do j=jb,je
jp = j+1
jm = j-1
do i=ib,ie
im = i-1
ip = i+1
! t_vav(i,j,k) = 0.5*( (tvmx(i,j,k) - ( &
dummyx(i,j,k) = ( &
( numol * (uav(i+1,j,k)-uav(i,j,k))*dxfi(i) &
-numol * (uav(i,j,k)-uav(i-1,j,k))*dxfi(i-1) ) * 2. * dxhi(i) &
+ &
( numol * ( (uav(i,jp,k)-uav(i,j,k)) *dyi &
+(vav(i,jp,k)-vav(i-1,jp,k))*dxhi(i)) &
- numol * ( (uav(i,j,k)-uav(i,jm,k)) *dyi &
+(vav(i,j,k)-vav(i-1,j,k)) *dxhi(i)) &
) * dyi &
+ &
( numol * ( (uav(i,j,kp)-uav(i,j,k)) *dzhi(kp) &
+(wav(i,j,kp)-wav(i-1,j,kp))*dxhi(i)) &
- numol * ( (uav(i,j,k)-uav(i,j,km)) *dzhi(k) &
+(wav(i,j,k)-wav(i-1,j,k)) *dxhi(i)) &
) *dzfi(k) )
! y-direction
dummyy(i,j,k) = ( &
( numol * ( (vav(i+1,j,k)-vav(i,j,k)) *dxhi(i+1) &
+(uav(i+1,j,k)-uav(i+1,jm,k))*dyi) &
-numol * ( (vav(i,j,k)-vav(i-1,j,k)) *dxhi(i) &
+(uav(i,j,k)-uav(i,jm,k)) *dyi) &
) * dxfi(i) & ! = d/dx( Km*(dv/dx + du/dy) )
+ &
(numol * (vav(i,jp,k)-vav(i,j,k)) &
-numol * (vav(i,j,k)-vav(i,jm,k)) ) * 2. * dy2i & ! =d/dy( 2*Km*(dv/dy) )
+ &
( numol * ( (vav(i,j,kp)-vav(i,j,k)) *dzhi(kp) &
+(wav(i,j,kp)-wav(i,jm,kp)) *dyi) &
-numol * ( (vav(i,j,k)-vav(i,j,km)) *dzhi(k) &
+(wav(i,j,k)-wav(i,jm,k)) *dyi) &
) * dzfi(k) ) ! = d/dz( Km*(dv/dz + dw/dy) )
! z-direction
dummyz(i,j,k) = ( &
( numol * ( (wav(i+1,j,k)-wav(i,j,k)) *dxhi(i+1) &
+(uav(i+1,j,k)-uav(i+1,j,km)) *dzhi(k) ) &
-numol * ( (wav(i,j,k)-wav(i-1,j,k)) *dxhi(i) &
+(uav(i,j,k)-uav(i,j,km)) *dzhi(k) ) &
)*dxfi(i) &
+ &
( numol * ( (wav(i,jp,k)-wav(i,j,k)) *dyi &
+(vav(i,jp,k)-vav(i,jp,km)) *dzhi(k) ) &
-numol * ( (wav(i,j,k)-wav(i,jm,k)) *dyi &
+(vav(i,j,k)-vav(i,j,km)) *dzhi(k) ) &
)*dyi &
+ &
( numol * (wav(i,j,kp)-wav(i,j,k)) *dzfi(k) &
-numol * (wav(i,j,k)-wav(i,j,km)) *dzfi(km) ) * 2. &
* dzhi(k))
strainav2 = ( &
((uav(ip,j,k)-uav(i,j,k)) *dxfi(i) )**2 + &
((vav(i,jp,k)-vav(i,j,k)) *dyi )**2 + &
((wav(i,j,kp)-wav(i,j,k)) *dzfi(k) )**2 )
strainav2 = strainav2 + 0.125 * ( &
((wav(i,j,kp)-wav(im,j,kp)) *dxhi(i) + &
(uav(i,j,kp)-uav(i,j,k)) *dzhi(kp) )**2 + &
((wav(i,j,k)-wav(im,j,k)) *dxhi(i) + &
(uav(i,j,k)-uav(i,j,km)) *dzhi(k) )**2 + &
((wav(ip,j,k)-wav(i,j,k)) *dxhi(ip) + &
(uav(ip,j,k)-uav(ip,j,km)) *dzhi(k) )**2 + &
((wav(ip,j,kp)-wav(i,j,kp)) *dxhi(ip) + &
(uav(ip,j,kp)-uav(ip,j,k)) *dzhi(kp) )**2 )
strainav2 = strainav2 + 0.125 * ( &
((uav(i,jp,k)-uav(i,j,k)) *dyi + &
(vav(i,jp,k)-vav(im,jp,k)) *dxhi(i) )**2 + &
((uav(i,j,k)-uav(i,jm,k)) *dyi + &
(vav(i,j,k)-vav(im,j,k)) *dxhi(i) )**2 + &
((uav(ip,j,k)-uav(ip,jm,k)) *dyi + &
(vav(ip,j,k)-vav(i,j,k)) *dxhi(ip) )**2 + &
((uav(ip,jp,k)-uav(ip,j,k)) *dyi + &
(vav(ip,jp,k)-vav(i,jp,k)) *dxhi(ip) )**2 )
strainav2 = strainav2 + 0.125 * ( &
((vav(i,j,kp)-vav(i,j,k)) *dzhi(kp) + &
(wav(i,j,kp)-wav(i,jm,kp)) *dyi )**2 + &
((vav(i,j,k)-vav(i,j,km)) *dzhi(k)+ &
(wav(i,j,k)-wav(i,jm,k)) *dyi )**2 + &
((vav(i,jp,k)-vav(i,jp,km)) *dzhi(k)+ &
(wav(i,jp,k)-wav(i,j,k)) *dyi )**2 + &
((vav(i,jp,kp)-vav(i,jp,k)) *dzhi(kp) + &
(wav(i,jp,kp)-wav(i,j,kp)) *dyi )**2 )
dissresav(i,j,k) = 2.*numol *(strain2av(i,j,k) - strainav2) !resolved dissipation
end do
end do
end do
! BC's
tvmx (ie+1,:,:) = tvmx (ie,:,:)
tsgsmx1(ie+1,:,:) = tsgsmx1(ie,:,:)
tsgsmx2(ie+1,:,:) = tsgsmx2(ie,:,:)
dummyx (ie+1,:,:) = dummyx (ie,:,:)
ttmx (ie+1,:,:) = ttmx (ie,:,:)
call excjs( tvmy , ib,ie,jb,je,kb,ke,0,1) ! jb-1 is not used
call excjs( tsgsmy1, ib,ie,jb,je,kb,ke,0,1) ! jb-1 is not used
call excjs( tsgsmy2, ib,ie,jb,je,kb,ke,0,1) ! jb-1 is not used
call excjs( dummyy, ib,ie,jb,je,kb,ke,0,1) ! jb-1 is not used
call excjs( ttmy , ib,ie,jb,je,kb,ke,0,1) ! jb-1 is not used
tvmz (:,:,ke+1) = tvmz (:,:,ke)
tsgsmz1(:,:,ke+1) = tsgsmz1(:,:,ke)
tsgsmz2(:,:,ke+1) = tsgsmz2(:,:,ke)
dummyz (:,:,ke+1) = dummyz(:,:,ke)
ttmz (:,:,ke+1) = ttmz (:,:,ke)
do k=kb,ke
km = k-1
kp = k+1
do j=jb,je
jp = j+1
jm = j-1
do i=ib,ie
im = i-1
ip = i+1
! Total viscous dissipation
t_vav(i,j,k) = 0.5*( (tvmx(i, j,k) - dummyx(i,j,k) *uav(i, j,k)) + &
(tvmx(ip,j,k) - dummyx(ip,j,k)*uav(ip,j,k))) &
+ 0.5*( (tvmy(i,j, k) - dummyy(i,j,k) *vav(i,j, k)) + &
(tvmy(i,jp,k) - dummyy(i,jp,k)*vav(i,jp,k))) &
+ 0.5*( (tvmz(i,j,k ) - dummyz(i,j,k) *wav(i,j,k )) + &
(tvmz(i,j,kp) - dummyz(i,j,kp)*wav(i,j,kp))) &
+ dissresav(i,j,k) ! d/dxj(2*nu*<ui'Sij'>) = <u_i*d/dxj(2*nu*Sij')> +2*nu*<Sij'Sij'>
! Now the same for subgrid stress
! <d/dxj(2*u_i'*nu_t*Sij)'> = <u_i'*d/dxj(2*nu_t*Sij)'> + <(2*nu_t*Sij)'*Sij'>
! = <u_i*d/dxj(2*nu_t*Sij)> -
! <u_i>*<d/dxj(2*nu_t*Sij)> + <2*nu_t*Sij*Sij> -
! <(2*nu_t*Sij)>*<Sij>
! = <u_i*d/dxj(2*nu_t*Sij)> -
! <u_i>*<d/dxj(2*nu_t*Sij)> + <2*nu_t*Sij*Sij> -
! 2*<nu_t>*<Sij>*<Sij> - 2*<nu_t'*Sij'>*<Sij>
!---------------------------------------
!Total subgrid TKE
!---------------------------------------
! Mean SGS dissipation
disssgsfl(i,j,k) = 2.*nusgsav(i,j,k)*strainav2 ! = 2*<nu_sgs>*<sij>*<sij>
! TKE
tke(i,j,k) = 0.5*(0.5*(upupav(ip,j,k)+upupav(i,j,k)) + &
0.5*(vpvpav(i,jp,k)+vpvpav(i,j,k)) + &
0.5*(wpwpav(i,j,kp)+wpwpav(i,j,k)))
! total SGS
t_sgsav(i,j,k) = 0.5*( (tsgsmx1(i,j,k) - uav(i,j,k) *tsgsmx2(i,j,k)) + &
(tsgsmx1(ip,j,k) - uav(ip,j,k)*tsgsmx2(ip,j,k))) &
+ & ! = <2*nu_t*SijSij> - <2*nu_t*Sij>*<Sij>
0.5*( (tsgsmy1(i,j,k) - vav(i,j,k) *tsgsmy2(i,j,k)) + &
(tsgsmy1(i,jp,k) - vav(i,jp,k)*tsgsmy2(i,jp,k))) &
+ & ! = <2*nu_t*SijSij> - <2*nu_t*Sij>*<Sij>
0.5*( (tsgsmz1(i,j,k) - vav(i,j,k) *tsgsmz2(i,j,k)) + &
(tsgsmz1(i,j,kp) - vav(i,j,kp)*tsgsmz2(i,j,kp))) &
+ disssgsav(i,j,k) - disssgsfl(i,j,k)
! -2*<nu_t'Sij'>*<Sij> should still be added!
! SGS dissipation
d_sgsav(i,j,k)= - disssgsav(i,j,k) + disssgsfl(i,j,k)
! +2*<nu_t'Sij'>*<Sij> should still be added! (is compensated with above)
!---------------------------------------
!Total pressure TKE
!---------------------------------------
! Pressure correlation term
! - <uj'*dp'/dxj> = - <uj*dp/dxj> + <uj>*d<p>/dxj
t_pav(i,j,k) = tpm(i,j,k) + &
0.5*(uav(i,j,k)*(presav(i,j,k)-presav(i-1,j,k))*dxhi(i) + &
uav(i+1,j,k)*(presav(i+1,j,k)-presav(i,j,k))*dxhi(i+1)) &
+ &
0.5*(vav(i,j,k)*(presav(i,j,k)-presav(i,j-1,k))*dyi + &
vav(i,j+1,k)*(presav(i,j+1,k)-presav(i,j,k))*dyi) &
+ &
0.5*(wav(i,j,k)*(presav(i,j,k)-presav(i,j,k-1))*dzhi(k) + &
wav(i,j,k+1)*(presav(i,j,k+1)-presav(i,j,k))*dzhi(k+1))
! - d/dxj(<0.5*ui'ui'uj'>) = -<uj'd/dxj(<0.5*ui'ui'>) + <ui'uj'><Sij>
! = -<uj*d/dxj(0.5*ui'ui')> + <uj>*d/dxj(<0.5*ui'ui'> +
! <ui'uj'><Sij>)
! ttav(i,j,k) = ttm(i,j,k) -
!---------------------------------------
!Total advection TKE
!---------------------------------------
! <advection term N.S. times ui> = MKE + A - Pshear - Tt
! Tt = -<ui'd/dxj(ui'uj')> = -<d/dxj(0.5*ui'ui'uj')> = A + MKE - Pshear - Total
!Pshear =Ptav = -<ui'uj'>d/dxj(<Sij>) = -<ui'uj'>d<ui>/dxj
! mechanical or shear production
p_tav(i,j,k) = - ( &
0.5 *(upupav(i,j,k)+upupav(ip,j,k))* (uav(ip,j,k)-uav(i,j,k))*dxfi(i) + & ! <u'u'>*d<u>/dx
0.25*(upvpav(i,j,k) *(uav(i, j, k)-uav(i, jm,k) )*dyi + &
upvpav(i,jp,k) *(uav(i, jp,k)-uav(i, j, k) )*dyi + &
upvpav(ip,j,k) *(uav(ip,j, k)-uav(ip,jm,k) )*dyi + &
upvpav(ip,jp,k)*(uav(ip,jp,k)-uav(ip,j, k) )*dyi) + & ! <u'v'>*d<u>/dy
0.25*(upwpav(i, j,k ) *(uav(i, j,k )-uav(i,j,km))*dzhi(k) + &
upwpav(i, j,kp) *(uav(i, j,kp)-uav(i, j,k))*dzhi(kp) + &
upwpav(ip,j,k ) *(uav(ip,j,k )-uav(ip,j,km))*dzhi(k) + &
upwpav(ip,j,kp) *(uav(ip,j,kp)-uav(ip,j,k))*dzhi(kp)) + & ! <u'w'>*d<u>/dz
0.25*(upvpav(i, j, k) *(vav(i, j, k)-vav(im,j,k))*dxhi(i) + &
upvpav(ip,j, k) *(vav(ip,j, k)-vav(i, j,k))*dxhi(ip) + &
upvpav(i, jp,k) *(vav(i,jp,k)-vav(im,jp,k))*dxhi(i) + &
upvpav(ip,jp,k) *(vav(ip,jp,k)-vav(i,jp,k))*dxhi(ip)) + & ! <u'v'>*d<v>/dx
0.5 *(vpvpav(i,j,k)+vpvpav(i,jp,k))*(vav(i,jp,k)-vav(i,j,k))*dyi + & ! <v'v'>*d<v>/dy
0.5 *(vpvpav(i,j,k)+vpvpav(i,jp,k))*(vav(i,jp,k)-vav(i,j,k))*dyi + & ! <v'v'>*d<v>/dy
0.25*(vpwpav(i,j ,k ) *(vav(i,j ,k )-vav(i,j,km))*dzhi(k) + &
vpwpav(i,j ,kp) *(vav(i,j ,kp)-vav(i,j ,k))*dzhi(kp) + &
vpwpav(i,jp,k ) *(vav(i,jp,k)-vav(i,jp,km))*dzhi(k) + &
vpwpav(i,jp,kp) *(vav(i,jp,kp)-vav(i,jp,k))*dzhi(kp)) + & ! <v'w'>*d<v>/dz
0.25*(upwpav(i, j, k) *(wav(i, j,k )-wav(im,j,k))*dxhi(i) + &
upwpav(ip,j, k) *(wav(ip,j,k )-wav(i, j,k))*dxhi(ip) + &
upwpav(i, j,kp) *(wav(i,j,kp)-wav(im,j,kp))*dxhi(i) + &
upwpav(ip,j,kp) *(wav(ip,j,kp)-wav(i,j,kp))*dxhi(ip)) + & ! <u'w'>*d<w>/dx
0.25*(vpwpav(i,j,k) *(wav(i,j, k )-wav(i,jm,k ) )*dyi + &
vpwpav(i,jp,k) *(wav(i,jp,k )-wav(i,j, k ) )*dyi + &
vpwpav(ip,j,k) *(wav(i,j, kp)-wav(i,jm,kp) )*dyi + &
vpwpav(ip,jp,k)*(wav(i,jp,kp)-wav(i,j, kp) )*dyi) + & ! <v'w'>*d<w>/dy
0.5 *(wpwpav(i,j,k)+wpwpav(i,j,kp))*(wav(i,j,kp)-wav(i,j,k))*dzfi(k) ) ! <w'w'>*d<w>/dz
! Mean kinetic energy term (expected to be small).
mke(i,j,k) = 0.5*(uav(ip,j,k)+uav(i,j,k))*(uuav(ip,j,k)-uuav(i,j,k))*dxfi(i) + & !<u>*d<uu>/dx
0.5*(uav(i, j,k)*(uvav(i ,jp,k)-uvav(i ,j,k))*dyi + & ! <u>*d<uv>/dy
uav(ip,j,k)*(uvav(ip,jp,k)-uvav(ip,j,k))*dyi) + &
0.5*(uav(i, j,k)*(uwav(i ,j,kp)-uwav(i ,j,k))*dzfi(k) + & ! <u>*d<uw>/dz
uav(ip,j,k)*(uwav(ip,j,kp)-uwav(ip,j,k))*dzfi(k)) + &
0.5*(vav(i,j, k)*(uvav(ip,j ,k)-uvav(i,j ,k))*dxfi(i) + & ! <v>*d<uv>/dx
vav(i,jp,k)*(uvav(ip,jp,k)-uvav(i,jp,k))*dxfi(i)) + &
0.5*(vav(i,jp,k)+vav(i,j,k))*(vvav(i,jp,k)-vvav(i,j,k))*dyi + & ! <v>*d<vv>/dy
0.5*(vav(i,j ,k)*(vwav(i,j ,kp)-vwav(i,j ,k))*dzfi(k) + & ! <v>*d<vw>/dz
vav(i,jp,k)*(vwav(i,jp,kp)-vwav(i,jp,k))*dzfi(k)) + &
0.5*(wav(i,j,k )*(uwav(ip,j,k )-uwav(i,j,k ))*dxfi(i) + & ! <w>*d<uw>/dx
wav(i,j,kp)*(uwav(ip,j,kp)-uwav(i,j,kp))*dxfi(i)) + &
0.5*(wav(i,j,k )*(vwav(i,jp,k )-vwav(i,j,k ))*dyi + & ! <w>*d<vw>/dy
wav(i,j,kp)*(vwav(i,jp,kp)-vwav(i,j,kp))*dyi) + &
0.5*(wav(i,j,kp)+wav(i,j,k))*(wwav(i,j,kp)-wwav(i,j,k))*dzfi(k) ! <w>*d<ww>/dz
! Advection of TKE
tkeadv(i,j,k) = 0.5*(uav(i, j,k)*(tke(i, j,k)-tke(im,j,k))*dxhi(i) + & ! <u>*de/dx
uav(ip,j,k)*(tke(ip,j,k)-tke(i ,j,k))*dxhi(ip)) + & !
0.5*(vav(i, j,k)*(tke(i,j ,k)-tke(i,jm,k))*dyi + & ! <v>*de/dy
vav(i,jp,k)*(tke(i,jp,k)-tke(i,j ,k))*dyi) + &
0.5*(wav(i,j,k )*(tke(i,j,k )-tke(i,j,km))*dzhi(k) + & ! <w>*de/dz
wav(i,j,kp)*(tke(i,j,kp)-tke(i,j,k ))*dzhi(kp))
! <advection term N.S. times ui> = MKE + A - Pshear - Tt
! Tt = -<ui'd/dxj(ui'uj')> = -<d/dxj(0.5*ui'ui'uj')> = A + MKE -
! Pshear - Total
! = tkeadv + mke -
! p_tav - ttm
! t_tav(i,j,k) = tkeadv(i,j,k) + mke(i,j,k) - p_tav(i,j,k) - ttm(i,j,k)
t_tav(i,j,k) = tkeadv(i,j,k) + mke(i,j,k) - p_tav(i,j,k) &
- 0.5*(ttmx(i,j,k) + ttmx(ip,j,k)) &
- 0.5*(ttmy(i,j,k) + ttmy(i,jp,k)) &
- 0.5*(ttmz(i,j,k) + ttmz(i,j,kp))
p_bav(i,j,k) = (grav/thls)*0.5*(thlpwpav(i,j,k)+thlpwpav(i,j,kp)) !use of thls here...????
end do
end do
end do
! need updating tg3315
call avexy_ibm(p_b(kb:ke+kh),p_bav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(t_p(kb:ke+kh),t_pav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(adv(kb:ke+kh),tkeadv(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(t_t(kb:ke+kh),t_tav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(t_sgs(kb:ke+kh),t_sgsav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(p_t(kb:ke+kh),p_tav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(d_sgs(kb:ke+kh),d_sgsav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
call avexy_ibm(t_v(kb:ke+kh),t_vav(:,:,kb:ke+kh),ib,ie,jb,je,kb,ke,ih,jh,kh,IIc,IIcs,.true.)
end subroutine tkestatsdump