subroutine tkestats(tsamplep,tstatsdumpp) ! change of variable names not yet translated across to here ! tg3315 30/11/17
use modfields, only : u0,v0,w0,thlm,uyt,vyt,wyt,thlyt,pres0,&
tvmx,tvmy,tvmz,strain2av,tsgsmx1,tsgsmx2,tsgsmy1,tsgsmy2,&
tsgsmz1,tsgsmz2,pres0
use modglobal, only : ib,ie,ih,jb,je,jgb,jge,dy,jh,ke,kb,kh,rk3step,cexpnr,tsample,tstatsdump,dzf,zh,dxf,dzf,numol,&
dzfi,dxfi,dyi,dy2i,dxfiq,dxhiq,dyiq,dzfi5,dzh,dzf,dzhi,dzhiq,dxf,dxhi
use modstat_nc, only : writestat_nc
use modsurfdata, only : thls
use modsubgriddata, only : ekm
implicit none
real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: tekm ! turbulent viscosity
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: emom
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: eomm
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: eopm
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: epom
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: emmo
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: eomp
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: epmo
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: emop
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: empo
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: tkesgs
! real, dimension(ib-ih:ie+ih,jb-jh:je+jh,kb-kh:ke+kh) :: nusgs
! real :: dummy
integer i,j,k,im,ip,jm,jp,km,kp
real tstatsdumppi,tsamplep,tstatsdumpp,strain2,tkesgs,nusgs,&
emom,eomm,eopm,epom,emmo,eomp,epmo,emop,empo,dummy
tekm(:,:,:) = ekm(:,:,:) - numol
tstatsdumppi = 1./tstatsdumpp
!---------------------------------------
! Viscous transport TKE
!---------------------------------------
!> Time averaged viscous transport in x,y and z to be used to calculate the total viscous transport for TKE
! Tvmx at u-locations (ib:ih+ih:jb:je,kb:ke)
! This is similar to routine diffu time u_i
do k = kb,ke
kp=k+1
km=k-1
do j = jb,je
jp=j+1
jm=j-1
do i = ib,ie
ip=i+1
im=i-1
dummy = u0(i,j,k)*( &
( numol * (u0(i+1,j,k)-u0(i,j,k))*dxfi(i) &
-numol * (u0(i,j,k)-u0(i-1,j,k))*dxfi(i-1) ) * 2. * dxhi(i) &
+ &
( numol * ( (u0(i,jp,k)-u0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i-1,jp,k))*dxhi(i)) &
-numol * ( (u0(i,j,k)-u0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i-1,j,k)) *dxhi(i)) &
) * dyi &
+ &
( numol * ( (u0(i,j,kp)-u0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i-1,j,kp))*dxhi(i)) &
-numol * ( (u0(i,j,k)-u0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i-1,j,k)) *dxhi(i)) &
) *dzfi(k) )
tvmx(i,j,k)=(tvmx(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average tvmx
! Tvmv at v-locations (ib:ih:jb:je+1,kb:ke)
! This is similar to routine diffv time v
dummy = v0(i,j,k) * ( &
( numol * ( (v0(i+1,j,k)-v0(i,j,k)) *dxhi(i+1) &
+(u0(i+1,j,k)-u0(i+1,jm,k))*dyi) &
-numol * ( (v0(i,j,k)-v0(i-1,j,k)) *dxhi(i) &
+(u0(i,j,k)-u0(i,jm,k)) *dyi) &
) * dxfi(i) & ! = d/dx( Km*(dv/dx + du/dy) )
+ &
(numol * (v0(i,jp,k)-v0(i,j,k)) &
-numol * (v0(i,j,k)-v0(i,jm,k)) ) * 2. * dy2i & ! = d/dy( 2*Km*(dv/dy) )
+ &
(numol * ( (v0(i,j,kp)-v0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i,jm,kp)) *dyi) &
-numol * ( (v0(i,j,k)-v0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i,jm,k)) *dyi) &
) * dzfi(k) ) ! = d/dz( Km*(dv/dz + dw/dy) )
tvmy(i,j,k)=(tvmy(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average tvmy
! Tvmz at w-locations (ib:ih:jb:je,kb:ke+kh)
! This is similar to routine diffw time w
dummy = w0(i,j,k) * ( &
( numol * ( (w0(i+1,j,k)-w0(i,j,k)) *dxhi(i+1) &
+(u0(i+1,j,k)-u0(i+1,j,km)) *dzhi(k) ) &
-numol * ( (w0(i,j,k)-w0(i-1,j,k)) *dxhi(i) &
+(u0(i,j,k)-u0(i,j,km)) *dzhi(k) ) &
)*dxfi(i) &
+ &
( numol * ( (w0(i,jp,k)-w0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i,jp,km)) *dzhi(k) ) &
-numol * ( (w0(i,j,k)-w0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i,j,km)) *dzhi(k) ) &
)*dyi &
+ &
( numol * (w0(i,j,kp)-w0(i,j,k)) *dzfi(k) &
-numol * (w0(i,j,k)-w0(i,j,km)) *dzfi(km) ) * 2. &
* dzhi(k) )
tvmz(i,j,k)=(tvmz(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average uwsgsav
! Compute stresses and fluxes at c.c. , also used in total viscous transport
strain2 = ( &
((u0(ip,j,k)-u0(i,j,k)) *dxfi(i) )**2 + &
((v0(i,jp,k)-v0(i,j,k)) *dyi )**2 + &
((w0(i,j,kp)-w0(i,j,k)) *dzfi(k) )**2 )
strain2 = strain2 + 0.125 * ( &
((w0(i,j,kp)-w0(im,j,kp)) *dxhi(i) + &
(u0(i,j,kp)-u0(i,j,k)) *dzhi(kp) )**2 + &
((w0(i,j,k)-w0(im,j,k)) *dxhi(i) + &
(u0(i,j,k)-u0(i,j,km)) *dzhi(k) )**2 + &
((w0(ip,j,k)-w0(i,j,k)) *dxhi(ip) + &
(u0(ip,j,k)-u0(ip,j,km)) *dzhi(k) )**2 + &
((w0(ip,j,kp)-w0(i,j,kp)) *dxhi(ip) + &
(u0(ip,j,kp)-u0(ip,j,k)) *dzhi(kp) )**2 )
strain2 = strain2 + 0.125 * ( &
((u0(i,jp,k)-u0(i,j,k)) *dyi + &
(v0(i,jp,k)-v0(im,jp,k)) *dxhi(i) )**2 + &
((u0(i,j,k)-u0(i,jm,k)) *dyi + &
(v0(i,j,k)-v0(im,j,k)) *dxhi(i) )**2 + &
((u0(ip,j,k)-u0(ip,jm,k)) *dyi + &
(v0(ip,j,k)-v0(i,j,k)) *dxhi(ip) )**2 + &
((u0(ip,jp,k)-u0(ip,j,k)) *dyi + &
(v0(ip,jp,k)-v0(i,jp,k)) *dxhi(ip) )**2 )
strain2 = strain2 + 0.125 * ( &
((v0(i,j,kp)-v0(i,j,k)) *dzhi(kp) + &
(w0(i,j,kp)-w0(i,jm,kp)) *dyi )**2 + &
((v0(i,j,k)-v0(i,j,km)) *dzhi(k)+ &
(w0(i,j,k)-w0(i,jm,k)) *dyi )**2 + &
((v0(i,jp,k)-v0(i,jp,km)) *dzhi(k)+ &
(w0(i,jp,k)-w0(i,j,k)) *dyi )**2 + &
((v0(i,jp,kp)-v0(i,jp,k)) *dzhi(kp) + &
(w0(i,jp,kp)-w0(i,j,kp)) *dyi )**2 )
strain2av(i,j,k)=(strain2av(i,j,k)*(tstatsdumpp-tsamplep) + strain2*tsamplep)*tstatsdumppi ! update average strain2av
!--------------------------------------------------
!> SGS TKE
!--------------------------------------------------
! x-direction
emom = ( dzf(km) * ( tekm(i,j,k)*dxf(i-1) + tekm(i-1,j,k)*dxf(i) ) + & ! dx is non-equidistant
dzf(k) * ( tekm(i,j,km)*dxf(i-1) + tekm(i-1,j,km)*dxf(i) ) )*dxhi(i) * dzhiq(k)
emop = ( dzf(kp) * ( tekm(i,j,k)*dxf(i-1) + tekm(i-1,j,k)*dxf(i) ) + & ! dx is non-equidistant
dzf(k) * ( tekm(i,j,kp)*dxf(i-1) + tekm(i-1,j,kp)*dxf(i) ) )*dxhi(i) * dzhiq(kp)
empo = 0.25 * ( ( tekm(i,j,k)+tekm(i,jp,k))*dxf(i-1) + (tekm(i-1,j,k)+tekm(i-1,jp,k))*dxf(i) ) * dxhi(i) ! dx is non-equidistant
emmo = 0.25 * ( ( tekm(i,j,k)+tekm(i,jm,k))*dxf(i-1) +(tekm(i-1,jm,k)+tekm(i-1,j,k))*dxf(i) ) * dxhi(i) ! dx is non-equidistant
! dummy = u0(i,j,k)*( &
dummy = ( &
( tekm(i,j,k) * (u0(i+1,j,k)-u0(i,j,k))*dxfi(i) &
-tekm(i-1,j,k)* (u0(i,j,k)-u0(i-1,j,k))*dxfi(i-1) ) * 2. * dxhi(i) &
+ &
( empo * ( (u0(i,jp,k)-u0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i-1,jp,k))*dxhi(i)) &
-emmo * ( (u0(i,j,k)-u0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i-1,j,k)) *dxhi(i)) &
) * dyi &
+ &
( emop * ( (u0(i,j,kp)-u0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i-1,j,kp))*dxhi(i)) &
-emom * ( (u0(i,j,k)-u0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i-1,j,k)) *dxhi(i)) &
) *dzfi(k) )
tsgsmx1(i,j,k)=(tsgsmx1(i,j,k)*(tstatsdumpp-tsamplep) + dummy*u0(i,j,k)*tsamplep)*tstatsdumppi ! update average tsgsmx1
tsgsmx2(i,j,k)=(tsgsmx2(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average tsgsmx2
! y-direction
eomm = ( dzf(km) * ( tekm(i,j,k) + tekm(i,jm,k) ) + &
dzf(k) * ( tekm(i,j,km) + tekm(i,jm,km) ) ) * dzhiq(k)
eomp = ( dzf(kp) * ( tekm(i,j,k) + tekm(i,jm,k) ) + &
dzf(k) * ( tekm(i,j,kp) + tekm(i,jm,kp) ) ) * dzhiq(kp)
emmo = 0.25 * ( ( tekm(i,j,k)+tekm(i,jm,k))*dxf(i-1) +(tekm(i-1,jm,k)+tekm(i-1,j,k))*dxf(i) ) * dxhi(i) ! dx is non-equidistant
epmo = 0.25 * ( ( tekm(i,j,k)+tekm(i,jm,k))*dxf(i+1) + (tekm(i+1,jm,k)+tekm(i+1,j,k))*dxf(i) ) * dxhi(i+1) ! dx is non-equidistant
! dummy = v0(i,j,k) * ( &
dummy = ( &
( epmo * ( (v0(i+1,j,k)-v0(i,j,k)) *dxhi(i+1) &
+(u0(i+1,j,k)-u0(i+1,jm,k))*dyi) &
-emmo * ( (v0(i,j,k)-v0(i-1,j,k)) *dxhi(i) &
+(u0(i,j,k)-u0(i,jm,k)) *dyi) &
) * dxfi(i) & ! = d/dx( Km*(dv/dx + du/dy) )
+ &
(tekm(i,j,k) * (v0(i,jp,k)-v0(i,j,k)) &
-tekm(i,jm,k)* (v0(i,j,k)-v0(i,jm,k)) ) * 2. * dy2i & ! = d/dy( 2*Km*(dv/dy) )
+ &
( eomp * ( (v0(i,j,kp)-v0(i,j,k)) *dzhi(kp) &
+(w0(i,j,kp)-w0(i,jm,kp)) *dyi) &
-eomm * ( (v0(i,j,k)-v0(i,j,km)) *dzhi(k) &
+(w0(i,j,k)-w0(i,jm,k)) *dyi) &
) * dzfi(k) ) ! = d/dz( Km*(dv/dz + dw/dy) )
tsgsmy1(i,j,k)=(tsgsmy1(i,j,k)*(tstatsdumpp-tsamplep) + dummy*v0(i,j,k)*tsamplep)*tstatsdumppi ! update average tsgsmy1 = <v*d/dxj(2*nu*S2j)>
tsgsmy2(i,j,k)=(tsgsmy2(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average tsgsmy2 = <d/dxj(2*nu*S2j)>
! z-direction
emom = ( dzf(km) * ( tekm(i,j,k)*dxf(i-1) + tekm(i-1,j,k)*dxf(i) )*dxhi(i) + &
dzf(k) * ( tekm(i,j,km)*dxf(i-1) + tekm(i-1,j,km)*dxf(i) )*dxhi(i) ) * dzhiq(k)
eomm = ( dzf(km) * ( tekm(i,j,k) + tekm(i,jm,k) ) + &
dzf(k) * ( tekm(i,j,km) + tekm(i,jm,km) ) ) * dzhiq(k)
eopm = ( dzf(km) * ( tekm(i,j,k) + tekm(i,jp,k) ) + &
dzf(k) * ( tekm(i,j,km) + tekm(i,jp,km) ) ) * dzhiq(k)
epom = ( dzf(km) * ( tekm(i,j,k)*dxf(i+1) + tekm(i+1,j,k)*dxf(i) )*dxhi(i+1) + &
dzf(k) * ( tekm(i,j,km)*dxf(i+1) + tekm(i+1,j,km)*dxf(i) )*dxhi(i+1) ) * dzhiq(k)
! dummy = w0(i,j,k) * ( &
dummy = ( &
( epom * ( (w0(i+1,j,k)-w0(i,j,k)) *dxhi(i+1) &
+(u0(i+1,j,k)-u0(i+1,j,km)) *dzhi(k) ) &
-emom * ( (w0(i,j,k)-w0(i-1,j,k)) *dxhi(i) &
+(u0(i,j,k)-u0(i,j,km)) *dzhi(k) ) &
)*dxfi(i) &
+ &
( eopm * ( (w0(i,jp,k)-w0(i,j,k)) *dyi &
+(v0(i,jp,k)-v0(i,jp,km)) *dzhi(k) ) &
-eomm * ( (w0(i,j,k)-w0(i,jm,k)) *dyi &
+(v0(i,j,k)-v0(i,j,km)) *dzhi(k) ) &
)*dyi &
+ &
( tekm(i,j,k) * (w0(i,j,kp)-w0(i,j,k)) *dzfi(k) &
-tekm(i,j,km)* (w0(i,j,k)-w0(i,j,km)) *dzfi(km) ) * 2. &
* dzhi(k))
tsgsmz1(i,j,k)=(tsgsmz1(i,j,k)*(tstatsdumpp-tsamplep) + dummy*w0(i,j,k)*tsamplep)*tstatsdumppi ! update average tsgsmz1 = <w*d/dxj(2*nu*S3j)>
tsgsmz2(i,j,k)=(tsgsmz2(i,j,k)*(tstatsdumpp-tsamplep) + dummy*tsamplep)*tstatsdumppi ! update average tsgsmz2 = <d/dxj(2*nu*S3j)>
end do
end do
end do
end subroutine tkestats