1664 lines
66 KiB
Fortran
1664 lines
66 KiB
Fortran
MODULE post
|
|
USE Compact
|
|
USE m_parameters
|
|
|
|
IMPLICIT NONE
|
|
PRIVATE
|
|
INTEGER :: countnum
|
|
INTEGER :: num_, dummyu_ ! hybrid
|
|
REAL :: tnow
|
|
REAL :: hxp,hyp,hzp
|
|
REAL, PARAMETER :: pi=3.14159265358979323846
|
|
REAL, PARAMETER :: me=1.00e-20
|
|
REAL, DIMENSION(:), ALLOCATABLE :: favg_ndata
|
|
!REAL, DIMENSION(:,:), ALLOCATABLE :: SMF,SMC,FMS,KMS
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: SMF,SMC
|
|
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: CM_b,CM_u,RMS,RMS_b,RMS_u,COV
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: RMS_g,RMS_b_g,RMS_u_g,COV_g
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: RMS_b_2,RMS_u_2 !uPrime
|
|
|
|
!REAL, DIMENSION(:,:), ALLOCATABLE :: IRE1,IRE2,IRE3,IRE4,IRE5,IRE6,IRE7,IRE8,IRE9 !edge-cold-bc-3
|
|
REAL, DIMENSION(:,:), ALLOCATABLE :: IRE1,IRE3,IRE4,IRE5,IRE6
|
|
|
|
!REAL, DIMENSION(:,:,:), ALLOCATABLE :: cgm,c,Wc,y,DivN,Div_V !,pres
|
|
!REAL, DIMENSION(:,:,:), ALLOCATABLE :: c_dot,c_dot_g,c_g,FSD_dot
|
|
!REAL, DIMENSION(:,:,:), ALLOCATABLE :: vn,GN_vn,GN_sd,G2N_c
|
|
!REAL, DIMENSION(:,:,:), ALLOCATABLE :: G2_Y
|
|
REAL, DIMENSION(:,:,:), ALLOCATABLE :: c,Wc,y
|
|
REAL, DIMENSION(:,:,:), ALLOCATABLE :: c_dot,c_dot_g,c_g,FSD_dot
|
|
|
|
REAL, DIMENSION(:,:,:), ALLOCATABLE :: u,v,w
|
|
|
|
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: old_scalar, new_scalar
|
|
|
|
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: m_v,NV,G_C,sd,u_dot,G_V,G_FSD,ub_dot,uu_dot
|
|
!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: G_vn,G_sd,G_Y !G_Y is added for edge-cold-bc-3
|
|
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: m_v_new
|
|
|
|
REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: CM_b_g,CM_u_g,u_g,u_dot_g,ub_dot_g,uu_dot_g
|
|
!REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: d2gc
|
|
|
|
PUBLIC :: main
|
|
CONTAINS
|
|
|
|
SUBROUTINE main
|
|
INTEGER :: fread,i
|
|
|
|
CALL READ_INTRO
|
|
|
|
CALL ALLOCATE_ARRAYS
|
|
|
|
CALL PRINT_BANNER
|
|
|
|
countnum=0
|
|
|
|
|
|
firstloop: DO fread=startnum,endnum,skipnum
|
|
|
|
IF ( to_omit(fread) ) THEN
|
|
|
|
WRITE(*,'(a47,i7,a4,i5,a3,i5)') &
|
|
' Current fullsavenum = ', fread, ' || ', (fread-startnum+1), ' / ', (endnum-startnum+1)
|
|
|
|
WRITE(*,'(a12,i6,a20,i6)') &
|
|
' Skip. ', omit_t(i,1), ' <= fullsavenum <= ', omit_t(i,2)
|
|
|
|
ELSE
|
|
countnum=countnum+1
|
|
|
|
CALL READ_FILE(fread)
|
|
CALL CAL_Yrs
|
|
! CALL CAL_CGM_N
|
|
! CALL CAL_Grad_Div
|
|
! CALL CAL_Sds
|
|
! CALL CAL_vnsd
|
|
CALL CAL_SUM ! Sum for each fort.xxxx
|
|
CALL SAVE_SUM ! Total sum
|
|
ENDIF
|
|
|
|
ENDDO firstloop
|
|
|
|
write(*,*) '1st loop finished'
|
|
|
|
CALL AVERAGING
|
|
|
|
|
|
|
|
|
|
secondloop: DO fread=startnum,endnum,skipnum
|
|
|
|
IF ( to_omit(fread) ) THEN
|
|
|
|
WRITE(*,'(a40,i7,a4,i4,a3,i4)') &
|
|
' Current fullsavenum = ', fread, ' || ', (fread-startnum+1), ' / ', (endnum-startnum+1)
|
|
|
|
WRITE(*,'(a12,i6,a20,i6)') &
|
|
' Skip. ', omit_t(i,1), ' <= fullsavenum <= ', omit_t(i,2)
|
|
|
|
ELSE
|
|
|
|
CALL READ_FILE(fread)
|
|
CALL CAL_Yrs
|
|
!CALL CAL_CGM_N !dhkim
|
|
CALL CAL_FLUCTUATION
|
|
CALL SAVE_SUM_FLUCTUATION
|
|
|
|
ENDIF
|
|
|
|
ENDDO secondloop
|
|
|
|
CALL FLUCTUATION_AVG
|
|
|
|
! CALL FINAL_AVG ! edge-cold-bc-3
|
|
CALL SAVE_AVG_RESULTS
|
|
|
|
CALL DEALLOCATES_CLOSE
|
|
|
|
WRITE(*,*) ' Avergaing RAW data is FINISHED'
|
|
WRITE(*,*) 'qEdge_X.dat is generated'
|
|
END SUBROUTINE main
|
|
|
|
!========================================================================================
|
|
! End of main routine
|
|
!========================================================================================
|
|
|
|
SUBROUTINE PRINT_BANNER
|
|
WRITE(*,*) ' This program, x-edge-cold-bc-5-hybrid, is written by D. Kim, 2018'
|
|
WRITE(*,*) ' It is to study the statistics of the flame parameters at the leading edge'
|
|
WRITE(*,*) ' in turbulent premixed flames.'
|
|
WRITE(*,'(a40,i5,a11,i5,a1)') ' Postprocess will be done from "FORT.',startnum,'" to "FORT.',endnum,'"'
|
|
END SUBROUTINE PRINT_BANNER
|
|
|
|
SUBROUTINE READ_FILE(num)
|
|
INTEGER, INTENT(IN) :: num
|
|
REAL, DIMENSION(2) :: tmpr
|
|
INTEGER :: nx, ny, nz
|
|
REAL :: tmp1,tmp2
|
|
REAL :: dt,dummyu
|
|
INTEGER :: ncyc
|
|
|
|
OPEN(num,FORM='unformatted',STATUS='unknown')
|
|
|
|
READ (num) tnow,nx,ny,nz,tmp1,tmp2
|
|
|
|
IF ((nx .ne. nxp) .or. (ny .ne. nyp) .or. (nz .ne. nzp)) THEN
|
|
WRITE(0,*) "Array dimension mismatch", nx, ny, nz, " != ", nxp, nyp, nzp
|
|
STOP -1
|
|
ENDIF
|
|
|
|
READ (num) ncyc,dt,dummyu
|
|
READ (num) tmpr(1:2)
|
|
READ (num) tmpr(1:2)
|
|
READ (num) tmpr(1:2)
|
|
|
|
WRITE(*,'(a40,f8.3,a2,i7,a2,i5,a4,i5,a3,i5)') ' Current time / NCYC / fullsavenum = ',tnow,&
|
|
' /',ncyc,' /',num,' || ',(num-startnum+1),' / ',(endnum-startnum+1)
|
|
WRITE(*,*) ' Reading current data file and processing'
|
|
|
|
m_v(1,:,:,:)=1. ! density is fixed at 1.
|
|
|
|
num_=num
|
|
|
|
IF(num.le.shiftnum) THEN
|
|
WRITE(*,*) ' with an old fort data from Nueman-0X'
|
|
READ (num) u,v,w,old_scalar
|
|
m_v(6,:,:,:) = old_scalar(2,:,:,:)
|
|
ELSE
|
|
WRITE(*,*) ' with a new fort data from Comb-Cluster'
|
|
READ (num) u,v,w,new_scalar
|
|
u = u + dummyu
|
|
m_v_new(:,:,:,2) = new_scalar(:,:,:,2)
|
|
ENDIF
|
|
|
|
m_v(2,:,:,:) = u
|
|
m_v(3,:,:,:) = v
|
|
m_v(4,:,:,:) = w
|
|
|
|
CLOSE (num)
|
|
|
|
END SUBROUTINE READ_FILE
|
|
|
|
|
|
SUBROUTINE CAL_Yrs
|
|
INTEGER :: i,j,k
|
|
REAL :: wrate,yi,rpr
|
|
rpr=1./prp/lep ! 1/Sc
|
|
|
|
refwr=pre*1.*exp(-ac/(1.+bc*c_ref)) ! wrate at c_ref
|
|
|
|
DO k=1,nzp
|
|
DO j=1,nyp
|
|
DO i=1,nxp
|
|
|
|
IF(num_.le.shiftnum) THEN
|
|
yi=m_v(6,i,j,k)
|
|
ELSE
|
|
yi=m_v_new(i,j,k,2)
|
|
ENDIF
|
|
|
|
c(i,j,k)=1.-yi
|
|
y(i,j,k)=yi
|
|
if(c(i,j,k).lt.0.) c(i,j,k)=0. ! 141011
|
|
if(c(i,j,k).gt.1.) c(i,j,k)=1. ! 141011
|
|
|
|
! wrate=pre*yi*exp(-ac/(1.+bcc*(1.-yi))) !wrate
|
|
! IF ((1.-yi).le.ccut) THEN
|
|
! wrate=pre*yi*exp(-cs/(1.+bcc*(1.-yi))) !wrate
|
|
! ENDIF
|
|
|
|
wrate=pre*yi*exp(-ac/(1.+bc*(1.-yi))) !wrate
|
|
|
|
! Cold boundary treatment, Kwon ------------------------------------------------
|
|
! min_wr=0.
|
|
|
|
IF (c(i,j,k).le.c_ref) THEN
|
|
wrate=min_wr
|
|
IF (c(i,j,k).gt.c_cut) wrate=((refwr-min_wr)*exp(prof_wr*(c(i,j,k)-c_ref))+ &
|
|
min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref)))
|
|
ENDIF
|
|
! ------------------------------------------------------------------------------
|
|
Wc(i,j,k)=wrate
|
|
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
|
|
rod=vis0p*rpr ! vis0p dynamic viscosity, rod= density*mass diffusivity = rho*D
|
|
|
|
END SUBROUTINE CAL_Yrs
|
|
|
|
! SUBROUTINE CAL_CGM_N
|
|
! REAL, DIMENSION(2,nxp) :: ux !,dux
|
|
! REAL, DIMENSION(2,nyp) :: uy !,duy
|
|
! REAL, DIMENSION(2,nzp) :: uz !,duz
|
|
! REAL, DIMENSION(3,nxp) :: dux
|
|
! REAL, DIMENSION(3,nyp) :: duy
|
|
! REAL, DIMENSION(3,nzp) :: duz
|
|
! INTEGER :: i,j,k
|
|
!
|
|
! ux=0.;dux=0.;uy=0.;duy=0.;uz=0.;duz=0.;G_C=0.; cgm=0.; NV=0.; G_Y=0.; d2gc=0.
|
|
!
|
|
!
|
|
!!$omp parallel do private(ux,dux,uy,duy)
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! ux(1,i)=c(i,j,k)
|
|
! ux(2,i)=1.-c(i,j,k) !edge-cold-bc-3
|
|
! ENDDO
|
|
! !CALL dfnonp(nxp,hxp,ux,dux(1,:),1,1)
|
|
! CALL dfnonp(nxp,hxp,ux(1:2,:),dux(1:2,:),2,1) !edge-cold-bc-3
|
|
! DO i=1,nxp
|
|
! IF(c(i,j,k).le.0.) dux(1,i)=0.
|
|
! G_C(1,i,j,k)=dux(1,i) ! dc/dx
|
|
! G_Y(1,i,j,k)=dux(2,i) ! d(1-c/)dx
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
! DO i=1,nxp
|
|
! DO j=1,nyp
|
|
! uy(1,j)=c(i,j,k)
|
|
! uy(2,j)=1.-c(i,j,k) !edge-cold-bc-3
|
|
! ENDDO
|
|
! !CALL dfp(nyp,hyp,uy,duy(1,:),1,2)
|
|
! CALL dfp(nyp,hyp,uy(1:2,:),duy(1:2,:),2,2) !edge-cold-bc-3
|
|
! DO j=1,nyp
|
|
! IF(c(i,j,k).le.0.) duy(1,j)=0.
|
|
! G_C(2,i,j,k)=duy(1,j) ! dc/dy
|
|
! G_Y(2,i,j,k)=duy(2,j) ! d(1-c/)dy
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
! IF (twod.eq.0) THEN
|
|
!!$omp parallel do private(uz,duz)
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! DO k=1,nzp
|
|
! uz(1,k)=c(i,j,k)
|
|
! uz(2,k)=1.-c(i,j,k) !edge-cold-bc-3
|
|
! ENDDO
|
|
! !CALL dfp(nzp,hzp,uz,duz(1,:),1,3)
|
|
! CALL dfp(nzp,hzp,uz(1:2,:),duz(1:2,:),2,3) !edge-cold-bc-3
|
|
! DO k=1,nzp
|
|
! IF(c(i,j,k).le.0.) duz(1,k)=0.
|
|
! G_C(3,i,j,k)=duz(1,k) ! dc/dz
|
|
! G_Y(3,i,j,k)=duz(2,k) ! d(1-c)/dz
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDIF
|
|
!
|
|
!!$omp parallel do
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! cgm(i,j,k)=SQRT(G_C(1,i,j,k)*G_C(1,i,j,k)+ &
|
|
! G_C(2,i,j,k)*G_C(2,i,j,k)+ &
|
|
! G_C(3,i,j,k)*G_C(3,i,j,k) ) ! |Grad(c)|
|
|
!! IF (c(i,j,k).gt.min_c.and.c(i,j,k).le.1.) THEN
|
|
! NV(1,i,j,k)=-G_C(1,i,j,k)/cgm(i,j,k) ! Nx
|
|
! NV(2,i,j,k)=-G_C(2,i,j,k)/cgm(i,j,k) ! Ny
|
|
! NV(3,i,j,k)=-G_C(3,i,j,k)/cgm(i,j,k) ! Nz
|
|
! IF(cgm(i,j,k).eq.0.) NV(1:3,i,j,k)=0.
|
|
!! ELSE
|
|
!! NV(1:3,i,j,k)=0.
|
|
!! ENDIF
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
!!---edge-cold-bc-3---------------------------------------------------
|
|
!
|
|
!!!$omp parallel do private(ux,dux,uy,duy)
|
|
!! DO k=1,nzp
|
|
!! DO j=1,nyp
|
|
!! DO i=1,nxp
|
|
!! ux(1,i)=G_C(1,i,j,k)
|
|
!! ENDDO
|
|
!! CALL dfnonp(nxp,hxp,ux,dux(1,:),1,1)
|
|
!! DO i=1,nxp
|
|
!! IF(c(i,j,k).le.0.) dux(1,i)=0.
|
|
!! d2gc(2,i,j,k)=dux(1,i) ! d2c/dx2
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+dux(1,i) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!!
|
|
!! DO i=1,nxp
|
|
!! DO j=1,nyp
|
|
!! uy(1,j)=G_C(2,i,j,k)
|
|
!! ENDDO
|
|
!! CALL dfp(nyp,hyp,uy,duy(1,:),1,2)
|
|
!! DO j=1,nyp
|
|
!! IF(c(i,j,k).le.0.) duy(1,j)=0.
|
|
!! d2gc(3,i,j,k)=duy(1,j) ! d2c/dc2
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+duy(1,j) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!!
|
|
!! IF (twod.eq.0) THEN
|
|
!!!$omp parallel do private(uz,duz)
|
|
!! DO j=1,nyp
|
|
!! DO i=1,nxp
|
|
!! DO k=1,nzp
|
|
!! uz(1,k)=G_C(3,i,j,k)
|
|
!! ENDDO
|
|
!! CALL dfp(nzp,hzp,uz,duz(1,:),1,3)
|
|
!! DO k=1,nzp
|
|
!! IF(c(i,j,k).le.0.) duz(1,k)=0.
|
|
!! d2gc(4,i,j,k)=duz(1,k)
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+duz(1,k) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDIF
|
|
!
|
|
!!!$omp parallel do private(ux,dux,uy,duy)
|
|
!! DO k=1,nzp
|
|
!! DO j=1,nyp
|
|
!! DO i=1,nxp
|
|
!! ux(1,i)=c(i,j,k)
|
|
!! ENDDO
|
|
!! CALL d2fnonp(nxp,hxp,ux,dux(1,:),1,1)
|
|
!! DO i=1,nxp
|
|
!! IF(c(i,j,k).le.min_c.or.c(i,j,k).gt.1.) dux(1,i)=0.
|
|
!! d2gc(2,i,j,k)=dux(1,i) ! d2c/dx2
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+dux(1,i) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!!
|
|
!! DO i=1,nxp
|
|
!! DO j=1,nyp
|
|
!! uy(1,j)=c(i,j,k)
|
|
!! ENDDO
|
|
!! CALL d2fp(nyp,hyp,uy,duy(1,:),1,2)
|
|
!! DO j=1,nyp
|
|
!! IF(c(i,j,k).le.min_c.or.c(i,j,k).gt.1.) duy(1,j)=0.
|
|
!! d2gc(3,i,j,k)=duy(1,j) ! d2c/dc2
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+duy(1,j) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!!
|
|
!! IF (twod.eq.0) THEN
|
|
!!!$omp parallel do private(uz,duz)
|
|
!! DO j=1,nyp
|
|
!! DO i=1,nxp
|
|
!! DO k=1,nzp
|
|
!! uz(1,k)=c(i,j,k)
|
|
!! ENDDO
|
|
!! CALL d2fp(nzp,hzp,uz,duz(1,:),1,3)
|
|
!! DO k=1,nzp
|
|
!! IF(c(i,j,k).le.min_c.or.c(i,j,k).gt.1.) duz(1,k)=0.
|
|
!! d2gc(4,i,j,k)=duz(1,k)
|
|
!! d2gc(1,i,j,k)=d2gc(1,i,j,k)+duz(1,k) ! Local Laplacian(c)
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDDO
|
|
!! ENDIF
|
|
!!===edge-cold-bc-3===================================================
|
|
!
|
|
! END SUBROUTINE CAL_CGM_N
|
|
!
|
|
! SUBROUTINE CAL_Sds
|
|
! INTEGER :: i,j,k
|
|
! REAL, DIMENSION(4,nxp) :: ux,dux,d2ux
|
|
! REAL, DIMENSION(3,nyp) :: uy,duy,d2uy
|
|
! REAL, DIMENSION(3,nzp) :: uz,duz,d2uz
|
|
!
|
|
!! sd = sdr + sdd
|
|
! sd=0.
|
|
! d2gc=0.
|
|
!
|
|
! G2_Y=0.
|
|
!!========Sdr==================================================
|
|
!!$omp parallel do
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
!
|
|
! sd(2,i,j,k)=Wc(i,j,k)/cgm(i,j,k)/m_v(1,i,j,k) ! sdr
|
|
! IF (c(i,j,k).le.0.) sd(2,i,j,k)=0.
|
|
!
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!!=============================================================
|
|
!
|
|
!!========Sdn==================================================
|
|
!!$omp parallel do private(ux,dux,d2ux,uy,duy,d2uy)
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! ux(1,i)=rod ! rho*D
|
|
! ux(2,i)=c(i,j,k) ! c
|
|
! ux(3,i)=rod*cgm(i,j,k) ! rho*D*fsd'
|
|
! ux(4,i)=G_Y(1,i,j,k) ! d(1-c)/dx
|
|
! ENDDO
|
|
!
|
|
! CALL dfnonp(nxp,hxp,ux(1:4,:),dux(1:4,:),4,1)
|
|
! CALL d2fnonp(nxp,hxp,ux(2,:),d2ux(1,:),1,1)
|
|
!
|
|
! DO i=1,nxp
|
|
! IF(c(i,j,k).le.0.) THEN
|
|
! dux(1:4,i)=0.
|
|
! d2ux(1,i)=0.
|
|
! ENDIF
|
|
! sd(3,i,j,k)=sd(3,i,j,k)+(( ux(1,i)*d2ux(1,i)+dux(1,i)*dux(2,i) )) ! sdd
|
|
! sd(4,i,j,k)=sd(4,i,j,k) - NV(1,i,j,k)*dux(3,i) ! sdn
|
|
! d2gc(2,i,j,k)=d2ux(1,i) ! d2c/dx2
|
|
! d2gc(1,i,j,k)=d2gc(1,i,j,k)+d2ux(1,i) ! Local Laplacian(c)
|
|
! G2_Y(i,j,k)=dux(4,i) ! d2(1-c)/dx2
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
! DO i=1,nxp
|
|
! DO j=1,nyp
|
|
! uy(1,j)=rod ! rho*D
|
|
! uy(2,j)=c(i,j,k) ! c
|
|
! uy(3,j)=rod*cgm(i,j,k) ! rho*D*fsd'
|
|
! ENDDO
|
|
!
|
|
! CALL dfp(nyp,hyp,uy(1:3,:),duy(1:3,:),3,2)
|
|
! CALL d2fp(nyp,hyp,uy(2,:),d2uy(1,:),1,2)
|
|
!
|
|
! DO j=1,nyp
|
|
! IF(c(i,j,k).le.0.) THEN
|
|
! duy(1:3,j)=0.
|
|
! d2uy(1,j)=0.
|
|
! ENDIF
|
|
! sd(3,i,j,k)=sd(3,i,j,k)+(( uy(1,j)*d2uy(1,j)+duy(1,j)*duy(2,j) )) ! sdn
|
|
! sd(4,i,j,k)=sd(4,i,j,k) - NV(2,i,j,k)*duy(3,j) ! sdn
|
|
! d2gc(3,i,j,k)=d2uy(1,j) ! d2c/dc2
|
|
! d2gc(1,i,j,k)=d2gc(1,i,j,k)+d2uy(1,j) ! Local Laplacian(c)
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
! IF (twod.eq.0) THEN
|
|
!!$omp parallel do private(uz,duz,d2uz)
|
|
! DO i=1,nxp
|
|
! DO j=1,nyp
|
|
! DO k=1,nzp
|
|
! uz(1,k)=rod ! rho*D
|
|
! uz(2,k)=c(i,j,k) ! c
|
|
! uz(3,k)=rod*cgm(i,j,k) ! rho*D*fsd'
|
|
! ENDDO
|
|
!
|
|
! CALL dfp(nzp,hzp,uz(1:3,:),duz(1:3,:),3,3)
|
|
! CALL d2fp(nzp,hzp,uz(2,:),d2uz(1,:),1,3)
|
|
!
|
|
! DO k=1,nzp
|
|
! IF(c(i,j,k).le.0.) THEN
|
|
! duz(1:3,k)=0.
|
|
! d2uz(1,k)=0.
|
|
! ENDIF
|
|
! sd(3,i,j,k)=sd(3,i,j,k)+(( uz(1,k)*d2uz(1,k)+duz(1,k)*duz(2,k) ))
|
|
! sd(4,i,j,k)=sd(4,i,j,k) - NV(3,i,j,k)*duz(3,k) ! sdn
|
|
! d2gc(4,i,j,k)=d2uz(1,k)
|
|
! d2gc(1,i,j,k)=d2gc(1,i,j,k)+d2uz(1,k) ! Local Laplacian(c)
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDIF
|
|
!
|
|
!!$omp parallel do
|
|
! DO i=1,nxp
|
|
! DO j=1,nyp
|
|
! DO k=1,nzp
|
|
! sd(3:4,i,j,k)=sd(3:4,i,j,k)/m_v(1,i,j,k)/cgm(i,j,k)
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!!=============================================================
|
|
!
|
|
!!========Sd===================================================
|
|
!!$omp parallel do
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! sd(1,i,j,k)=sd(2,i,j,k)+sd(3,i,j,k)
|
|
! sd(4,i,j,k)=sd(4,i,j,k)+sd(2,i,j,k) ! sdn
|
|
! sd(5,i,j,k)=-rod/m_v(1,i,j,k)*DivN(i,j,k) ! sdt
|
|
! sd(6,i,j,k)=sd(4,i,j,k)+sd(5,i,j,k)-sd(2,i,j,k) ! sdd2
|
|
! sd(7,i,j,k)=sd(4,i,j,k)+sd(5,i,j,k) ! sd2
|
|
! IF(c(i,j,k).le.0.) sd(1:7,i,j,k)=0.
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
!! write(*,*) 'sdd, sdd2, sd1, sd2'
|
|
!! do i=1,nxp
|
|
!! write(*,'(4e20.10)') sd(3,i,128,128),sd(6,i,128,128),sd(1,i,128,128),sd(7,i,128,128)
|
|
!! enddo
|
|
!
|
|
!
|
|
! END SUBROUTINE CAL_Sds
|
|
!
|
|
! SUBROUTINE CAL_vnsd
|
|
! REAL, DIMENSION(2,nxp) :: ux,dux
|
|
! REAL, DIMENSION(2,nyp) :: uy,duy
|
|
! REAL, DIMENSION(2,nzp) :: uz,duz
|
|
! INTEGER :: i,j,k
|
|
! REAL :: ui,vi,wi
|
|
!
|
|
!! Compute Vn
|
|
! vn=0.
|
|
!!$omp parallel do private(ui,vi,wi)
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! ui=m_v(2,i,j,k)/m_v(1,i,j,k) ! u
|
|
! vi=m_v(3,i,j,k)/m_v(1,i,j,k) ! v
|
|
! wi=m_v(4,i,j,k)/m_v(1,i,j,k) ! w
|
|
! vn(i,j,k)=ui*NV(1,i,j,k)+vi*NV(2,i,j,k)+wi*NV(3,i,j,k)
|
|
! IF (c(i,j,k).le.0.) vn(i,j,k)=0.
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
!
|
|
!! Compute Grad. of Vn and sd
|
|
!
|
|
!!$omp parallel do private(ux,dux,uy,duy)
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! ux=0.
|
|
! DO i=1,nxp
|
|
! ux(1,i)=sd(1,i,j,k) ! Sd
|
|
! ux(2,i)=vn(i,j,k) ! vn
|
|
! IF (c(i,j,k).le.0.) ux(1:2,i)=0.
|
|
! ENDDO
|
|
!
|
|
! CALL dfnonp(nxp,hxp,ux(1:2,:),dux(1:2,:),2,1)
|
|
!
|
|
! DO i=1,nxp
|
|
! IF (c(i,j,k).le.0.) dux(1:2,i)=0.
|
|
! G_sd(1,i,j,k)=dux(1,i) ! d(Sd)/dx
|
|
! G_vn(1,i,j,k)=dux(2,i) ! d(vn)/dx
|
|
! ENDDO
|
|
! ENDDO ! y-loop
|
|
!
|
|
! DO i=1,nxp
|
|
! uy=0.
|
|
! DO j=1,nyp
|
|
! uy(1,j)=sd(1,i,j,k) ! Sd
|
|
! uy(2,j)=vn(i,j,k) ! vn
|
|
! IF (c(i,j,k).le.0.) uy(1:2,j)=0.
|
|
! ENDDO
|
|
!
|
|
! CALL dfp(nyp,hyp,uy(1:2,:),duy(1:2,:),2,2) ! y-dir. is periodic in FPs
|
|
!
|
|
! DO j=1,nyp
|
|
! IF (c(i,j,k).le.0.) duy(1:2,j)=0.
|
|
! G_sd(2,i,j,k)=duy(1,j) ! d(Sd)/dy
|
|
! G_vn(2,i,j,k)=duy(2,j) ! d(vn)/dy
|
|
! ENDDO
|
|
! ENDDO ! x-loop
|
|
! ENDDO ! z-loop
|
|
!
|
|
! IF (twod.eq.0) THEN
|
|
!!$omp parallel do private(uz,duz)
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! uz=0.
|
|
! DO k=1,nzp
|
|
! uz(1,k)=sd(1,i,j,k) ! Sd
|
|
! uz(2,k)=vn(i,j,k) ! vn
|
|
! IF (c(i,j,k).le.0.) uz(1:2,k)=0.
|
|
! ENDDO
|
|
! CALL dfp(nzp,hzp,uz(1:2,:),duz(1:2,:),2,3) ! z-dir. is periodic in FPs
|
|
!
|
|
! DO k=1,nzp
|
|
! IF (c(i,j,k).le.0.) duz(1:2,k)=0.
|
|
! G_sd(3,i,j,k)=duz(1,k) ! d(Sd)/dy
|
|
! G_vn(3,i,j,k)=duz(2,k) ! d(vn)/dy
|
|
! ENDDO
|
|
! ENDDO ! x-lopp
|
|
! ENDDO ! y-loop
|
|
! ENDIF
|
|
!
|
|
!!$omp parallel do
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! GN_sd(i,j,k)=NV(1,i,j,k)*G_sd(1,i,j,k)+NV(2,i,j,k)*G_sd(2,i,j,k)+ &
|
|
! NV(3,i,j,k)*G_sd(3,i,j,k) ! d(sd)/dn
|
|
! GN_vn(i,j,k)=NV(1,i,j,k)*G_vn(1,i,j,k)+NV(2,i,j,k)*G_vn(2,i,j,k)+ &
|
|
! NV(3,i,j,k)*G_vn(3,i,j,k) ! d(vn)/dn
|
|
! G2N_c(i,j,k) = -(NV(1,i,j,k)*G_FSD(1,i,j,k)+NV(2,i,j,k)*G_FSD(2,i,j,k) &
|
|
! +NV(3,i,j,k)*G_FSD(3,i,j,k)) ! d2c/dn2 = -d(fsd')/dn
|
|
! ENDDO
|
|
! ENDDO
|
|
! ENDDO
|
|
! END SUBROUTINE CAL_vnsd
|
|
!
|
|
! SUBROUTINE CAL_Grad_Div
|
|
! REAL, DIMENSION(3,nxp) :: ux,dux
|
|
! REAL, DIMENSION(3,nyp) :: uy,duy
|
|
! REAL, DIMENSION(3,nzp) :: uz,duz
|
|
! INTEGER :: i,j,k
|
|
!
|
|
! Div_V=0. ; DivN=0. ; G_V=0. ; G_FSD=0.
|
|
!
|
|
!!$omp parallel do private(ux,dux,uy,duy)
|
|
! DO k=1,nzp
|
|
! DO j=1,nyp
|
|
! ux=0.
|
|
! DO i=1,nxp
|
|
! ux(1,i)=m_v(2,i,j,k)/m_v(1,i,j,k) ! u
|
|
! ux(2,i)=NV(1,i,j,k) ! Nx
|
|
! ux(3,i)=cgm(i,j,k) ! FSD`
|
|
! IF (c(i,j,k).le.0.) ux(2:3,i)=0.
|
|
! ENDDO
|
|
!
|
|
! CALL dfnonp(nxp,hxp,ux(1:3,:),dux(1:3,:),3,1)
|
|
!
|
|
! DO i=1,nxp
|
|
! IF (c(i,j,k).le.0.) dux(2:3,i)=0.
|
|
! G_V(1,i,j,k)=dux(1,i) ! du/dx
|
|
! Div_V(i,j,k)=dux(1,i) ! Div(V)
|
|
! DivN(i,j,k)=dux(2,i) ! d(Nx)/dx
|
|
! G_FSD(1,i,j,k)=dux(3,i) ! d(FSD`)/dx
|
|
! ENDDO
|
|
! ENDDO ! y-loop
|
|
!
|
|
! DO i=1,nxp
|
|
! uy=0.
|
|
! DO j=1,nyp
|
|
! uy(1,j)=m_v(3,i,j,k)/m_v(1,i,j,k) ! v
|
|
! uy(2,j)=NV(2,i,j,k) ! Ny
|
|
! uy(3,j)=cgm(i,j,k) ! FSD`
|
|
! IF (c(i,j,k).le.0.) uy(2:3,j)=0.
|
|
! ENDDO
|
|
!
|
|
! CALL dfp(nyp,hyp,uy(1:3,:),duy(1:3,:),3,2) ! y-dir. is periodic in FPs
|
|
!
|
|
! DO j=1,nyp
|
|
! IF (c(i,j,k).le.0.) duy(2:3,j)=0.
|
|
! G_V(2,i,j,k)=duy(1,j) ! dv/dy
|
|
! Div_V(i,j,k)=Div_V(i,j,k)+duy(1,j) ! Div(rho*V)
|
|
! DivN(i,j,k)=DivN(i,j,k)+duy(2,j) ! d(Ny)/dy
|
|
! G_FSD(2,i,j,k)=duy(3,j) ! d(FSD`)/dy
|
|
! ENDDO
|
|
! ENDDO ! x-loop
|
|
! ENDDO ! z-loop
|
|
!
|
|
! IF (twod.eq.0) THEN
|
|
!!$omp parallel do private(uz,duz)
|
|
! DO j=1,nyp
|
|
! DO i=1,nxp
|
|
! uz=0.
|
|
! DO k=1,nzp
|
|
! uz(1,k)=m_v(4,i,j,k)/m_v(1,i,j,k) ! w
|
|
! uz(2,k)=NV(3,i,j,k) ! Nz
|
|
! uz(3,k)=cgm(i,j,k) ! FSD`
|
|
! IF (c(i,j,k).le.0.) uz(2:3,k)=0.
|
|
! ENDDO
|
|
! CALL dfp(nzp,hzp,uz(1:3,:),duz(1:3,:),3,3) ! z-dir. is periodic in FPs
|
|
!
|
|
! DO k=1,nzp
|
|
! IF (c(i,j,k).le.0.) duz(2:3,k)=0.
|
|
! G_V(3,i,j,k)=duz(1,k) ! dv/dy
|
|
! Div_V(i,j,k)=Div_V(i,j,k)+duz(1,k) ! Div(rho*V)
|
|
! DivN(i,j,k)=DivN(i,j,k)+duz(2,k) ! d(Nz)/dz
|
|
! G_FSD(3,i,j,k)=duz(3,k) ! d(FSD`)/dz
|
|
! ENDDO
|
|
! ENDDO ! x-lopp
|
|
! ENDDO ! y-loop
|
|
! ENDIF
|
|
!
|
|
! END SUBROUTINE CAL_Grad_Div
|
|
|
|
|
|
SUBROUTINE CAL_SUM
|
|
INTEGER :: i,j,k,jj
|
|
REAL :: ui,vi,wi
|
|
|
|
SMF=0.; SMC=0.
|
|
CM_b=0.; CM_u=0. !; CM_b_g=0.; CM_u_g=0.
|
|
|
|
DO i=1,nxp
|
|
DO j=syp,eyp
|
|
jj=j-syp+1
|
|
DO k=1,nzp
|
|
ui=m_v(2,i,j,k)/m_v(1,i,j,k) ! u
|
|
vi=m_v(3,i,j,k)/m_v(1,i,j,k) ! v
|
|
wi=m_v(4,i,j,k)/m_v(1,i,j,k) ! w
|
|
|
|
!WRITE(210,'(4e20.10)') REAL(i),REAL(j),REAL(k),ui
|
|
|
|
! For local averages
|
|
u_g(1,i,jj,k)=u_g(1,i,jj,k)+m_v(2,i,j,k)/m_v(1,i,j,k) ! Sum [u(i,j,k)]
|
|
u_g(2,i,jj,k)=u_g(2,i,jj,k)+m_v(3,i,j,k)/m_v(1,i,j,k) ! Sum [v(i,j,k)]
|
|
u_g(3,i,jj,k)=u_g(3,i,jj,k)+m_v(4,i,j,k)/m_v(1,i,j,k) ! Sum [w(i,j,k)]
|
|
! c_g(i,jj,k)=c_g(i,jj,k)+c(i,j,k)
|
|
|
|
! IF (c(i,j,k).gt.min_c.and.c(i,j,k).le.max_c) THEN ! 131103
|
|
CM_b_g(1,i,jj,k)=CM_b_g(1,i,jj,k)+c(i,j,k)*m_v(2,i,j,k)/m_v(1,i,j,k) ! sum [c*u]
|
|
CM_b_g(2,i,jj,k)=CM_b_g(2,i,jj,k)+c(i,j,k)*m_v(3,i,j,k)/m_v(1,i,j,k) ! sum [c*v]
|
|
CM_b_g(3,i,jj,k)=CM_b_g(3,i,jj,k)+c(i,j,k)*m_v(4,i,j,k)/m_v(1,i,j,k) ! sum [c*w]
|
|
CM_b_g(4,i,jj,k)=CM_b_g(4,i,jj,k)+c(i,j,k)*m_v(1,i,j,k) ! Sum [c*rho]
|
|
|
|
CM_u_g(1,i,jj,k)=CM_u_g(1,i,jj,k)+(1.-c(i,j,k))*m_v(2,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*u]
|
|
CM_u_g(2,i,jj,k)=CM_u_g(2,i,jj,k)+(1.-c(i,j,k))*m_v(3,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*v]
|
|
CM_u_g(3,i,jj,k)=CM_u_g(3,i,jj,k)+(1.-c(i,j,k))*m_v(4,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*w]
|
|
CM_u_g(4,i,jj,k)=CM_u_g(4,i,jj,k)+(1.-c(i,j,k))*m_v(1,i,j,k) ! Sum [(1-c)*rho]
|
|
! ENDIF
|
|
|
|
! For Simple averages, < >
|
|
SMF(1,i) =SMF(1,i)+ui ! Sum [u]
|
|
SMF(2,i) =SMF(2,i)+vi ! Sum [v]
|
|
SMF(3,i) =SMF(3,i)+wi ! Sum [w]
|
|
! SMF(4,i) =SMF(4,i)+m_v(1,i,j,k) ! Sum [rho]
|
|
! SMF(5,i) =SMF(5,i)+0. ! Sum [T]
|
|
! SMF(6,i) =SMF(6,i)+rod/m_v(1,i,j,k) ! Sum [D]
|
|
! SMF(7,i) =SMF(7,i)+G_C(1,i,j,k) ! Sum [dc/dx]
|
|
! SMF(8,i) =SMF(8,i)+G_C(2,i,j,k) ! Sum [dc/dy]
|
|
! SMF(9,i) =SMF(9,i)+G_C(3,i,j,k) ! Sum [dc/dz]
|
|
! SMF(10,i)=SMF(10,i)+G_V(1,i,j,k) ! Sum [du/dx]
|
|
! SMF(11,i)=SMF(11,i)+G_V(2,i,j,k) ! Sum [dv/dy]
|
|
! SMF(12,i)=SMF(12,i)+G_V(3,i,j,k) ! Sum [dw/dz]
|
|
! SMF(13,i)=SMF(13,i)+Div_V(i,j,k) ! Sum [Div(V)]
|
|
! SMF(14,i)=SMF(14,i)+G2N_c(i,j,k) ! Sum [ d2c/dn2 ]
|
|
! SMF(15,i) =SMF(15,i)+G_Y(1,i,j,k) ! Sum [d(1-c)/dx]
|
|
! SMF(16,i) =SMF(16,i)+G_Y(2,i,j,k) ! Sum [d(1-c)/dy]
|
|
! SMF(17,i) =SMF(17,i)+G_Y(3,i,j,k) ! Sum [d(1-c)/dz]
|
|
! SMF(18,i) =SMF(18,i)+d2gc(1,i,j,k) ! Sum [Lap.(c)]
|
|
! SMF(19,i) =SMF(19,i)+d2gc(2,i,j,k) ! Sum [d2c/dx2]
|
|
! SMF(20,i) =SMF(20,i)+d2gc(3,i,j,k) ! Sum [d2c/dy2]
|
|
! SMF(21,i) =SMF(21,i)+d2gc(4,i,j,k) ! Sum [d2c/dz2]
|
|
! SMF(22,i) =SMF(22,i)+NV(1,i,j,k)*G2N_c(i,j,k) ! Sum [ nx * d2c/dn2 ]
|
|
! SMF(23,i) =SMF(23,i)+G2_Y(i,j,k) ! Sum [d2(1-c)/dx2 ]
|
|
|
|
SMC(1,i)=SMC(1,i)+c(i,j,k) ! Sum [c]
|
|
!SMC(2,i)=SMC(2,i)+cgm(i,j,k) ! Sum [FSD`] 131031
|
|
SMC(3,i)=SMC(3,i)+y(i,j,k) ! Sum [y]
|
|
SMC(4,i)=SMC(4,i)+Wc(i,j,k)/m_v(1,i,j,k) ! Sum [Wc/rho]
|
|
|
|
|
|
! For surface averages, < >f
|
|
! IF (c(i,j,k).gt.min_c.and.c(i,j,k).le.max_c) THEN ! 131103
|
|
! For conditional average
|
|
! For burned quantities
|
|
CM_b(1,i)=CM_b(1,i)+c(i,j,k)*m_v(2,i,j,k)/m_v(1,i,j,k) ! sum [c*u]
|
|
CM_b(2,i)=CM_b(2,i)+c(i,j,k)*m_v(3,i,j,k)/m_v(1,i,j,k) ! sum [c*v]
|
|
CM_b(3,i)=CM_b(3,i)+c(i,j,k)*m_v(4,i,j,k)/m_v(1,i,j,k) ! sum [c*w]
|
|
CM_b(4,i)=CM_b(4,i)+c(i,j,k)*m_v(1,i,j,k) ! Sum [c*rho]
|
|
|
|
! For unburned quantities
|
|
CM_u(1,i)=CM_u(1,i)+(1.-c(i,j,k))*m_v(2,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*u]
|
|
CM_u(2,i)=CM_u(2,i)+(1.-c(i,j,k))*m_v(3,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*v]
|
|
CM_u(3,i)=CM_u(3,i)+(1.-c(i,j,k))*m_v(4,i,j,k)/m_v(1,i,j,k) ! sum [(1-c)*w]
|
|
CM_u(4,i)=CM_u(4,i)+(1.-c(i,j,k))*m_v(1,i,j,k) ! Sum [(1-c)*rho]
|
|
! ENDIF ! 131031
|
|
|
|
! For surface averages, < >f , < >k
|
|
! IF (c(i,j,k).gt.min_c) THEN ! 131103
|
|
! favg_ndata(i)=favg_ndata(i)+1.
|
|
! FMS(1,i) =FMS(1,i)+cgm(i,j,k) ! Sum [ FSD` ]
|
|
! FMS(2,i) =FMS(2,i)+sd(1,i,j,k)*cgm(i,j,k) ! Sum [Sd*FSD`]
|
|
! FMS(3,i) =FMS(3,i)+sd(2,i,j,k)*cgm(i,j,k) ! Sum [Sdr*FSD`]
|
|
! FMS(4,i) =FMS(4,i)+sd(3,i,j,k)*cgm(i,j,k) ! Sum [Sdd*FSD`]
|
|
! FMS(5,i) =FMS(5,i)+sd(4,i,j,k)*cgm(i,j,k) ! Sum [Sdn*FSD`]
|
|
! FMS(6,i) =FMS(6,i)+sd(5,i,j,k)*cgm(i,j,k) ! Sum [Sdt*FSD`]
|
|
! FMS(7,i) =FMS(7,i)+sd(6,i,j,k)*cgm(i,j,k) ! Sum [Sdd2*FSD`]
|
|
! FMS(8,i) =FMS(8,i)+sd(7,i,j,k)*cgm(i,j,k) ! Sum [Sd2*FSD`]
|
|
! FMS(9,i) =FMS(9,i)+DivN(i,j,k)*cgm(i,j,k) ! Sum [DivN*FSD`]
|
|
! FMS(10,i)=FMS(10,i)+abs(DivN(i,j,k))*cgm(i,j,k) ! Sum [|DivN|*FSD`]
|
|
! FMS(11,i)=FMS(11,i)+G_C(1,i,j,k)/c(i,j,k)*cgm(i,j,k) ! Sum [(dc/dx)/c*FSD`]
|
|
! FMS(12,i)=FMS(12,i)+cgm(i,j,k)**2./c(i,j,k) ! Sum [-(dc/dn)/c*FSD`]
|
|
! FMS(13,i)=FMS(13,i)+G_FSD(1,i,j,k) ! Sum [(1/FSD`)*(d(FSD`)/dx)*FSD`]
|
|
! FMS(14,i)=FMS(14,i)+G_FSD(2,i,j,k) ! Sum [(1/FSD`)*(d(FSD`)/dy)*FSD`]
|
|
! FMS(15,i)=FMS(15,i)+G_FSD(3,i,j,k) ! Sum [(1/FSD`)*(d(FSD`)/dz)*FSD`]
|
|
! FMS(16,i)=FMS(16,i)+NV(1,i,j,k)*cgm(i,j,k) ! Sum [Nx * FSD`]
|
|
! FMS(17,i)=FMS(17,i)+NV(2,i,j,k)*cgm(i,j,k) ! Sum [Ny * FSD`]
|
|
! FMS(18,i)=FMS(18,i)+NV(3,i,j,k)*cgm(i,j,k) ! Sum [Nz * FSD`]
|
|
! FMS(19,i)=FMS(19,i)+ui*cgm(i,j,k) ! Sum [u * FSD`]
|
|
! FMS(20,i)=FMS(20,i)+vi*cgm(i,j,k) ! Sum [v * FSD`]
|
|
! FMS(21,i)=FMS(21,i)+wi*cgm(i,j,k) ! Sum [w * FSD`]
|
|
! FMS(22,i)=FMS(22,i)+vn(i,j,k)*cgm(i,j,k) ! Sum [ vn * FSD`]
|
|
! FMS(23,i)=FMS(23,i)+(NV(1,i,j,k)**2.+NV(2,i,j,k)**2.+NV(3,i,j,k)**2.)*cgm(i,j,k) ! Sum [(N dot N)*FSD`]
|
|
! FMS(24,i)=FMS(24,i)+GN_vn(i,j,k)*cgm(i,j,k) ! Sum [d(vn)/dn * FSD`]
|
|
! FMS(25,i)=FMS(25,i)+GN_sd(i,j,k)*cgm(i,j,k) ! Sum [d(Sd)/dn * FSD`]
|
|
! FMS(26,i)=FMS(26,i)+(NV(1,i,j,k)*G_FSD(1,i,j,k)+&
|
|
! NV(2,i,j,k)*G_FSD(2,i,j,k)+&
|
|
! NV(3,i,j,k)*G_FSD(3,i,j,k)) ! Sum [-(N dot (grad.(FSD`))/FSD`) * FSD`]
|
|
!
|
|
! KMS(1,i) =KMS(1,i) + vn(i,j,k)*G2N_c(i,j,k) ! Sum [ Vn * d2c/dn2 ]
|
|
! KMS(2,i) =KMS(2,i) + sd(1,i,j,k)*G2N_c(i,j,k) ! Sum [ Sd * d2c/dn2 ]
|
|
! ENDIF
|
|
|
|
ENDDO ! z-loop
|
|
ENDDO ! y-loop
|
|
ENDDO ! x-loop
|
|
|
|
END SUBROUTINE CAL_SUM
|
|
|
|
SUBROUTINE SAVE_SUM
|
|
INTEGER :: i
|
|
!#############################################################
|
|
!####### Abbreviation ##################
|
|
!
|
|
! Lap. A - Laplacian of A
|
|
! Grad. A - Gradient of A, cf.) variable name -> G_X means Grad. X
|
|
!
|
|
!#############################################################
|
|
|
|
! IRE1 : 1 / 2 / 3 / 4 / 5 / 6 / 7 /
|
|
! <C> / <FSD> / <Yr> / <u> / <v> / <w> / <rho> /
|
|
|
|
! IRE2 : 1 / 2 / 3 /
|
|
! <rho*Sd>f / <rho*Sdr>f / <rho*Sdn>f /
|
|
|
|
! IRE3 : 1 / 2 / 3 / 4 / 5 / 6 /
|
|
! RMS(u') / RMS(v') / RMS(w') / RMS(u')_g / RMS(v')_g / RMS(w')_g /
|
|
|
|
! IRE4 : 1 / 2 / 3 / 4 /
|
|
! <u>b / <v>b / <w>b / <rho>b /
|
|
|
|
! IRE5 : 1 / 2 / 3 / 4 /
|
|
! <u>u / <v>u / <w>u / <rho>u /
|
|
|
|
DO i=1,nxp
|
|
IRE1(1,i)=IRE1(1,i)+SMC(1,i) ! Sum [c]
|
|
IRE1(2,i)=IRE1(2,i)+SMC(2,i) ! Sum [FSD`]
|
|
IRE1(3,i)=IRE1(3,i)+SMC(3,i) ! Sum [yr]
|
|
IRE1(4,i)=IRE1(4,i)+SMC(4,i) ! Sum [Wc/rho]
|
|
|
|
IRE1(5,i) =IRE1(5,i)+SMF(1,i) ! Sum [u]
|
|
IRE1(6,i) =IRE1(6,i)+SMF(2,i) ! Sum [v]
|
|
IRE1(7,i) =IRE1(7,i)+SMF(3,i) ! Sum [w]
|
|
! IRE1(8,i) =IRE1(8,i)+SMF(4,i) ! Sum [rho]
|
|
! IRE1(9,i) =IRE1(9,i)+SMF(5,i) ! Sum [T] - dummy
|
|
! IRE1(10,i)=IRE1(10,i)+SMF(6,i) ! Sum [D]
|
|
! IRE1(11,i)=IRE1(11,i)+SMF(7,i) ! Sum [dc/dx]
|
|
! IRE1(12,i)=IRE1(12,i)+SMF(8,i) ! Sum [dc/dy]
|
|
! IRE1(13,i)=IRE1(13,i)+SMF(9,i) ! Sum [dc/dz]
|
|
! IRE1(14,i)=IRE1(14,i)+SMF(10,i) ! Sum [du/dx]
|
|
! IRE1(15,i)=IRE1(15,i)+SMF(11,i) ! Sum [dv/dy]
|
|
! IRE1(16,i)=IRE1(16,i)+SMF(12,i) ! Sum [dw/dz]
|
|
! IRE1(17,i)=IRE1(17,i)+SMF(13,i) ! Sum [Div(V)]
|
|
! IRE1(19,i)=IRE1(19,i)+SMF(14,i) ! Sum [ d2c/dn2 ]
|
|
! IRE1(20,i)=IRE1(20,i)+SMF(15,i) ! Sum [d(1-c)/dx]
|
|
! IRE1(21,i)=IRE1(21,i)+SMF(16,i) ! Sum [d(1-c)/dy]
|
|
! IRE1(22,i)=IRE1(22,i)+SMF(17,i) ! Sum [d(1-c)/dz]
|
|
! IRE1(23,i)=IRE1(23,i)+SMF(18,i) ! Sum [Lap.(c)]
|
|
! IRE1(24,i)=IRE1(24,i)+SMF(19,i) ! Sum [d2c/dx2]
|
|
! IRE1(25,i)=IRE1(25,i)+SMF(20,i) ! Sum [d2c/dy2]
|
|
! IRE1(26,i)=IRE1(26,i)+SMF(21,i) ! Sum [d2c/dz2]
|
|
!
|
|
! IRE1(27,i)=IRE1(27,i)+SMF(22,i) ! Sum [ Nx*(d2c/dn2) ]
|
|
! IRE1(28,i)=IRE1(28,i)+SMF(23,i) ! Sum [d2(1-c)/dx2]
|
|
!
|
|
! IRE2(1,i) =IRE2(1,i) +FMS(1,i) ! Sum [ FSD` ]
|
|
! IRE2(2,i) =IRE2(2,i) +FMS(2,i) ! Sum [Sd*FSD`]
|
|
! IRE2(3,i) =IRE2(3,i) +FMS(3,i) ! Sum [Sdr*FSD`]
|
|
! IRE2(4,i) =IRE2(4,i) +FMS(4,i) ! Sum [Sdd*FSD`]
|
|
! IRE2(5,i) =IRE2(5,i) +FMS(5,i) ! Sum [Sdn*FSD`]
|
|
! IRE2(6,i) =IRE2(6,i) +FMS(6,i) ! Sum [Sdt*FSD`]
|
|
! IRE2(7,i) =IRE2(7,i) +FMS(7,i) ! Sum [Sdd2*FSD`]
|
|
! IRE2(8,i) =IRE2(8,i) +FMS(8,i) ! Sum [Sd2*FSD`]
|
|
! IRE2(9,i) =IRE2(9,i) +FMS(9,i) ! Sum [DivN*FSD`]
|
|
! IRE2(10,i)=IRE2(10,i)+FMS(10,i) ! Sum [|DivN|*FSD`]
|
|
! IRE2(11,i)=IRE2(11,i)+FMS(11,i) ! Sum [(dc/dx)/c*FSD`]
|
|
! IRE2(12,i)=IRE2(12,i)+FMS(12,i) ! Sum [-(dc/dn)/c*FSD`]
|
|
! IRE2(13,i)=IRE2(13,i)+FMS(13,i) ! Sum [(1/FSD`)*(d(FSD`)/dx)*FSD`]
|
|
! IRE2(14,i)=IRE2(14,i)+FMS(14,i) ! Sum [(1/FSD`)*(d(FSD`)/dy)*FSD`]
|
|
! IRE2(15,i)=IRE2(15,i)+FMS(15,i) ! Sum [(1/FSD`)*(d(FSD`)/dz)*FSD`]
|
|
! IRE2(16,i)=IRE2(16,i)+FMS(16,i) ! Sum [Nx*FSD`]
|
|
! IRE2(17,i)=IRE2(17,i)+FMS(17,i) ! Sum [Ny*FSD`]
|
|
! IRE2(18,i)=IRE2(18,i)+FMS(18,i) ! Sum [Nz*FSD`]
|
|
! IRE2(19,i)=IRE2(19,i)+FMS(19,i) ! Sum [u*FSD`]
|
|
! IRE2(20,i)=IRE2(20,i)+FMS(20,i) ! Sum [v*FSD`]
|
|
! IRE2(21,i)=IRE2(21,i)+FMS(21,i) ! Sum [w*FSD`]
|
|
! IRE2(22,i)=IRE2(22,i)+FMS(22,i) ! Sum [ (V dot N) *FSD`]
|
|
! IRE2(23,i)=IRE2(23,i)+FMS(23,i) ! Sum [ (N dot N) *FSD`]
|
|
! IRE2(24,i)=IRE2(24,i)+FMS(24,i) ! Sum [ d(vn)/dn * FSD`]
|
|
! IRE2(25,i)=IRE2(25,i)+FMS(25,i) ! Sum [ d(sd)/dn * FSD`]
|
|
!
|
|
! IRE2(34,i)=IRE2(34,i)+FMS(26,i) ! Sum [ (N dot (grad.(FSD`))/FSD`) * FSD`]
|
|
|
|
IRE3(1,i)=IRE3(1,i)+CM_b(1,i) ! Sum [c*u]
|
|
IRE3(2,i)=IRE3(2,i)+CM_b(2,i) ! Sum [c*v]
|
|
IRE3(3,i)=IRE3(3,i)+CM_b(3,i) ! Sum [c*w]
|
|
IRE3(4,i)=IRE3(4,i)+CM_b(4,i) ! Sum [c*rho]
|
|
|
|
IRE4(1,i)=IRE4(1,i)+CM_u(1,i) ! Sum [(1-c)*u]
|
|
IRE4(2,i)=IRE4(2,i)+CM_u(2,i) ! Sum [(1-c)*v]
|
|
IRE4(3,i)=IRE4(3,i)+CM_u(3,i) ! Sum [(1-c)*w]
|
|
IRE4(4,i)=IRE4(4,i)+CM_u(4,i) ! Sum [(1-c)*rho]
|
|
|
|
! IRE7(1,i)=IRE7(1,i)+KMS(1,i) ! Sum [ vn * d2c/dn2 ]
|
|
! IRE7(2,i)=IRE7(2,i)+KMS(2,i) ! Sum [ sd * d2c/dn2 ]
|
|
|
|
ENDDO
|
|
|
|
END SUBROUTINE SAVE_SUM
|
|
|
|
SUBROUTINE AVERAGING
|
|
REAL, DIMENSION(1,nxp) :: ux,dux
|
|
INTEGER :: i,j,k,jj
|
|
REAL :: ndata,nfile
|
|
|
|
WRITE(*,*) 'countnum_AVG',countnum
|
|
nfile=REAL(countnum)
|
|
ndata=nfile*REAL((eyp-syp+1)*nzp)
|
|
write(*,'(a30,4i5)')'AVERAGING,syp,eyp,nzp,nfile',syp,eyp,nzp,nint(nfile)
|
|
|
|
! Simple average, < >
|
|
IRE1(1:7,:)=IRE1(1:7,:)/ndata
|
|
!IRE1(1:17,:)=IRE1(1:17,:)/ndata
|
|
!IRE1(19:28,:)=IRE1(19:28,:)/ndata
|
|
|
|
! Surface average, < >k (IRE7(1~6))
|
|
!IRE7(1:2,:)=IRE7(1:2,:)/ndata
|
|
! DO i=1,nxp
|
|
! IRE7(3,i) = IRE7(1,i) + IRE7(2,i) ! <vn+sd>k * <d2c/dn2>
|
|
! IRE7(4:5,i)= IRE7(1:2,i)/IRE1(19,i)
|
|
! IRE7(6,i) = IRE7(4,i) + IRE7(5,i)
|
|
!
|
|
! if(IRE1(19,i).eq.0.) IRE7(4:6,i)=0.
|
|
! ENDDO
|
|
|
|
|
|
!! Surface average, < >f ================================
|
|
! DO i=1,nxp
|
|
! IRE2(1,i)=IRE2(1,i)/favg_ndata(i) ! < fsd2 >
|
|
! if(favg_ndata(i).eq.0.) IRE2(1,i)=0.
|
|
! IF(IRE2(1,i).ne.0.) THEN
|
|
! IRE2(2:25,i)=IRE2(2:25,i)/favg_ndata(i)/IRE2(1,i)
|
|
! IRE2(34,i)=-1*IRE2(34,i)/favg_ndata(i)/IRE2(1,i)
|
|
! IRE2(26,i)=IRE2(24,i)+IRE2(25,i) ! <d(Vn+Sd)/dn>f
|
|
!
|
|
! IRE2(27,i)=IRE2(19,i)*IRE2(16,i)+IRE2(20,i)*IRE2(17,i)+IRE2(21,i)*IRE2(18,i) ! <V>f dot <N>f
|
|
! IRE2(28,i)=IRE2(22,i)-IRE2(27,i) ! <V" dot N">f
|
|
! IRE2(29,i)=IRE2(16,i)**2.+IRE2(17,i)**2.+IRE2(18,i)**2. ! <N>f dot <N>f
|
|
! IRE2(30,i)=IRE2(23,i)-IRE2(29,i) ! <N" dot N">f
|
|
! IRE2(31,i)=SQRT(IRE2(29,i)) ! |<N>f| = abs(<N>f)
|
|
! IRE2(32,i)=IRE2(2,i)+IRE2(22,i) ! <vn+Sd>f
|
|
! ELSE
|
|
! IRE2(2:32,i)=0. ! To avoid infinity
|
|
! ENDIF
|
|
! ENDDO
|
|
|
|
!! d<Nx>f/dx
|
|
! DO i=1,nxp
|
|
! ux(1,i)=IRE2(16,i) ! <Nx>f
|
|
! ENDDO
|
|
! CALL dfnonp(nxp,hxp,ux,dux,1,1)
|
|
! DO i=1,nxp
|
|
! IRE2(33,i)=dux(1,i) ! d<Nx>f/dx
|
|
! ENDDO
|
|
|
|
! Conditional average, < >b, < >u
|
|
IRE3(1:4,:)=IRE3(1:4,:)/ndata ! <cA>
|
|
IRE4(1:4,:)=IRE4(1:4,:)/ndata ! <(1-c)A>
|
|
|
|
DO i=1,nxp
|
|
IRE3(1:4,i)=IRE3(1:4,i)/IRE1(1,i) ! < >b
|
|
IRE4(1:4,i)=IRE4(1:4,i)/(1.-IRE1(1,i)) ! < >u
|
|
IF(IRE1(1,i).eq.0.) IRE3(1:4,i)=0.
|
|
IF(IRE1(1,i).eq.1.) IRE4(1:4,i)=0.
|
|
ENDDO
|
|
|
|
!$omp parallel do
|
|
DO i=1,nxp
|
|
DO jj=1,(eyp-syp+1)
|
|
DO k=1,nzp
|
|
!IRE1(18,i)=IRE1(18,i)+c_g(i,jj,k)
|
|
IRE3(5:8,i)=IRE3(5:8,i)+CM_b_g(1:4,i,jj,k)
|
|
IRE4(5:8,i)=IRE4(5:8,i)+CM_u_g(1:4,i,jj,k)
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
|
|
DO i=1,nxp
|
|
!IRE1(18,i)=IRE1(18,i)/ndata
|
|
IRE3(5:8,i)=IRE3(5:8,i)/ndata/IRE1(1,i) ! < >b_g
|
|
IRE4(5:8,i)=IRE4(5:8,i)/ndata/(1.-IRE1(1,i)) ! < >u_g
|
|
IF(IRE1(1,i).eq.0.) IRE3(5:8,i)=0.
|
|
IF(IRE1(1,i).eq.1.) IRE4(5:8,i)=0.
|
|
ENDDO
|
|
|
|
!$omp parallel do
|
|
DO i=1,nxp
|
|
DO jj=1,(eyp-syp+1)
|
|
DO k=1,nzp
|
|
CM_b_g(1:4,i,jj,k)=CM_b_g(1:4,i,jj,k)/nfile/IRE1(1,i) ! < >b_g
|
|
CM_u_g(1:4,i,jj,k)=CM_u_g(1:4,i,jj,k)/nfile/(1.-IRE1(1,i)) ! < >u_g
|
|
IF(IRE1(1,i).eq.0.) CM_b_g(1:4,i,jj,k)=0.
|
|
IF(IRE1(1,i).eq.1.) CM_u_g(1:4,i,jj,k)=0.
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
|
|
! Local averages
|
|
u_g=u_g/nfile
|
|
c_g=c_g/nfile
|
|
|
|
END SUBROUTINE AVERAGING
|
|
|
|
|
|
SUBROUTINE CAL_FLUCTUATION
|
|
INTEGER :: i,j,k,jj
|
|
|
|
DO i=1,nxp
|
|
DO j=syp,eyp
|
|
jj=j-syp+1
|
|
DO k=1,nzp
|
|
! Total mean based fluctuations
|
|
u_dot(1:3,i,jj,k)=m_v(2:4,i,j,k)/m_v(1,i,j,k)-IRE1(5:7,i) ! u',v',w'
|
|
|
|
ub_dot(1:3,i,jj,k)=m_v(2:4,i,j,k)/m_v(1,i,j,k)-IRE3(1:3,i) ! (u_b)',(v_b)',(w_b)'
|
|
uu_dot(1:3,i,jj,k)=m_v(2:4,i,j,k)/m_v(1,i,j,k)-IRE4(1:3,i) ! (u_u)',(v_u)',(w_u)'
|
|
!c_dot(i,jj,k)=c(i,j,k)-IRE1(1,i) ! c'
|
|
!FSD_dot(i,jj,k)=cgm(i,j,k)-IRE1(2,i) ! (FSD')' !dhkim
|
|
|
|
! Local mean based fluctuations
|
|
u_dot_g(1,i,jj,k)=m_v(2,i,j,k)/m_v(1,i,j,k)-u_g(1,i,jj,k) ! u' at each grid
|
|
u_dot_g(2,i,jj,k)=m_v(3,i,j,k)/m_v(1,i,j,k)-u_g(2,i,jj,k) ! u' at each grid
|
|
u_dot_g(3,i,jj,k)=m_v(4,i,j,k)/m_v(1,i,j,k)-u_g(3,i,jj,k) ! u' at each grid
|
|
|
|
ub_dot_g(1:3,i,jj,k)=m_v(2:4,i,j,k)/m_v(1,i,j,k)-CM_b_g(1:3,i,jj,k) ! (u_b)'_g,(v_b)'_g,(w_b)'_g
|
|
uu_dot_g(1:3,i,jj,k)=m_v(2:4,i,j,k)/m_v(1,i,j,k)-CM_u_g(1:3,i,jj,k) ! (u_u)'_g,(v_u)'_g,(w_u)'_g
|
|
!c_dot_g(i,jj,k)=c(i,j,k)-c_g(i,jj,k) ! c'_g
|
|
|
|
ENDDO
|
|
ENDDO
|
|
ENDDO
|
|
|
|
END SUBROUTINE CAL_FLUCTUATION
|
|
|
|
SUBROUTINE SAVE_SUM_FLUCTUATION
|
|
INTEGER :: i,j,k,jj
|
|
RMS=0.; RMS_g=0.; RMS_b=0.; RMS_u=0.; COV=0.
|
|
RMS_b_g=0.; RMS_u_g=0. ; COV_g=0.
|
|
RMS_b_2=0.; RMS_u_2=0.
|
|
|
|
!$omp parallel do
|
|
DO i=1,nxp
|
|
DO j=syp,eyp
|
|
jj=j-syp+1
|
|
DO k=1,nzp
|
|
RMS(1,i)=RMS(1,i)+u_dot(1,i,jj,k)**2. ! Sum [u'^2]
|
|
RMS(2,i)=RMS(2,i)+u_dot(2,i,jj,k)**2. ! Sum [v'^2]
|
|
RMS(3,i)=RMS(3,i)+u_dot(3,i,jj,k)**2. ! Sum [w'^2]
|
|
RMS(4,i)=RMS(4,i)+0.5*(u_dot(1,i,jj,k)**2.+u_dot(2,i,jj,k)**2.+&
|
|
u_dot(3,i,jj,k)**2.) ! Sum [tke]
|
|
!--------------------uPrime--------------------------------------------------------------------------------
|
|
! RMS(5,i)=RMS(5,i)+SQRT((u_dot(1,i,jj,k)**2.+u_dot(2,i,jj,k)**2.+&
|
|
! u_dot(3,i,jj,k)**2.)/3) ! Sum [SQRT(1/3(u'^2 + v'^2 + w'^2))]
|
|
! RMS(6,i)=RMS(6,i)+SQRT(u_dot(1,i,jj,k)**2) ! Sum [SQRT(u'^2)]
|
|
! RMS(7,i)=RMS(7,i)+SQRT(u_dot(2,i,jj,k)**2) ! Sum [SQRT(v'^2)]
|
|
! RMS(8,i)=RMS(8,i)+SQRT(u_dot(3,i,jj,k)**2) ! Sum [SQRT(w'^2)]
|
|
RMS(5,i)=RMS(5,i)+((u_dot(1,i,jj,k)**2.+u_dot(2,i,jj,k)**2.+&
|
|
u_dot(3,i,jj,k)**2.)/3) ! Sum [1/3(*u'^2 + v'^2 + w'^2)]
|
|
RMS(6,i)=RMS(6,i)+(u_dot(1,i,jj,k)**2) ! Sum [(u'^2)]
|
|
RMS(7,i)=RMS(7,i)+(u_dot(2,i,jj,k)**2) ! Sum [(v'^2)]
|
|
RMS(8,i)=RMS(8,i)+(u_dot(3,i,jj,k)**2) ! Sum [(w'^2)]
|
|
!====================uPrime================================================================================
|
|
|
|
! IF (c(i,j,k).gt.min_c) THEN
|
|
RMS_b(1:3,i)=RMS_b(1:3,i)+c(i,j,k)*ub_dot(1:3,i,jj,k)**2. ! Sum [c*[(u_b)'^2]]
|
|
RMS_b(4,i)=RMS_b(4,i)+c(i,j,k)*0.5*(ub_dot(1,i,jj,k)**2.+&
|
|
ub_dot(2,i,jj,k)**2.+ub_dot(3,i,jj,k)**2.) ! Sum [c*tke_b]
|
|
|
|
RMS_u(1:3,i)=RMS_u(1:3,i)+(1.-c(i,j,k))*uu_dot(1:3,i,jj,k)**2. ! Sum [(1-c)*[(u_u)'^2]]
|
|
RMS_u(4,i)=RMS_u(4,i)+(1.-c(i,j,k))*0.5*(uu_dot(1,i,jj,k)**2.+&
|
|
uu_dot(2,i,jj,k)**2.+uu_dot(3,i,jj,k)**2.) ! Sum [tke_u]
|
|
!--------------------uPrime--------------------------------------------------------------------------------
|
|
! RMS_b_2(1:3,i)=RMS_b_2(1:3,i)+c(i,j,k)*u_dot(1:3,i,jj,k)**2. ! Sum [c*[u'^2]]
|
|
! RMS_b_2(4,i)=RMS_b_2(4,i)+c(i,j,k)*(u_dot(1,i,jj,k)**2.+&
|
|
! u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3. ! Sum [c*1/3*[u'^2 + v'^2 + w'^2]]
|
|
!
|
|
! RMS_u_2(1:3,i)=RMS_u_2(1:3,i)+(1.-c(i,j,k))*u_dot(1:3,i,jj,k)**2. ! Sum [(1-c)*[u'^2]]
|
|
! RMS_u_2(4,i)=RMS_u_2(4,i)+(1.-c(i,j,k))*(u_dot(1,i,jj,k)**2.+&
|
|
! u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3. ! Sum [(1-c)*1/3*[u'^2 + v'^2 + w'^2]]
|
|
|
|
|
|
! RMS_b_2(1:3,i)=RMS_b_2(1:3,i)+c(i,j,k)*SQRT(u_dot(1:3,i,jj,k)**2.) ! Sum[c*SQRT[u'],c*[v'],c*[w']]
|
|
! RMS_b_2(4,i)=RMS_b_2(4,i)+c(i,j,k)*SQRT((u_dot(1,i,jj,k)**2.+&
|
|
! u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3.) ! Sum [c*rmsU`]
|
|
!
|
|
! RMS_u_2(1:3,i)=RMS_u_2(1:3,i)+(1.-c(i,j,k))*SQRT(u_dot(1:3,i,jj,k)**2.) ! Sum [(1-c)*SQRT[u'],[v'],[x']]
|
|
! RMS_u_2(4,i)=RMS_u_2(4,i)+(1.-c(i,j,k))*SQRT((u_dot(1,i,jj,k)**2.+&
|
|
! u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3.) ! Sum [(1-c)*rmsU`]
|
|
|
|
RMS_b_2(1:3,i)=RMS_b_2(1:3,i)+c(i,j,k)*(u_dot(1:3,i,jj,k)**2.)
|
|
! Sum[c*[u'^2],c*[v'^2],c*[w'^2]]
|
|
RMS_b_2(4,i)=RMS_b_2(4,i)+c(i,j,k)*((u_dot(1,i,jj,k)**2.+&
|
|
u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3.)
|
|
! Sum [c*(1/3*[u'^2 + v'^2 + w'^2])]
|
|
|
|
RMS_u_2(1:3,i)=RMS_u_2(1:3,i)+(1.-c(i,j,k))*(u_dot(1:3,i,jj,k)**2.)
|
|
! Sum [(1-c)*[u'],(1-c)*[v'],(1-c)*[x']]
|
|
RMS_u_2(4,i)=RMS_u_2(4,i)+(1.-c(i,j,k))*((u_dot(1,i,jj,k)**2.+&
|
|
u_dot(2,i,jj,k)**2.+u_dot(3,i,jj,k)**2.)/3.)
|
|
! Sum [(1-c)*1/3*[u'^2 + v'^2 + w'^2]]
|
|
!====================uPrime================================================================================
|
|
|
|
! COV(1,i)=COV(1,i)+u_dot(1,i,jj,k)*c_dot(i,jj,k) ! Sum [u'c']
|
|
! COV(2,i)=COV(2,i)+u_dot(2,i,jj,k)*c_dot(i,jj,k) ! Sum [v'c']
|
|
! COV(3,i)=COV(3,i)+u_dot(3,i,jj,k)*c_dot(i,jj,k) ! Sum [w'c']
|
|
! COV(4,i)=COV(4,i)+u_dot(1,i,jj,k)*FSD_dot(i,jj,k) ! Sum [u'(FSD')'] !dhkim
|
|
|
|
RMS_g(1,i)=RMS_g(1,i)+u_dot_g(1,i,jj,k)**2. ! Sum [u'^2_g]
|
|
RMS_g(2,i)=RMS_g(2,i)+u_dot_g(2,i,jj,k)**2. ! Sum [v'^2_g]
|
|
RMS_g(3,i)=RMS_g(3,i)+u_dot_g(3,i,jj,k)**2. ! Sum [w'^2_g]
|
|
RMS_g(4,i)=RMS_g(4,i)+0.5*(u_dot_g(1,i,jj,k)**2.+ &
|
|
u_dot_g(2,i,jj,k)**2.+u_dot_g(3,i,jj,k)**2.) ! Sum [tke_g]
|
|
|
|
RMS_b_g(1:3,i)=RMS_b_g(1:3,i)+c(i,j,k)*ub_dot_g(1:3,i,jj,k)**2. ! Sum [c*[(u_b)'^2]_g]
|
|
RMS_b_g(4,i)=RMS_b_g(4,i)+c(i,j,k)*0.5*(ub_dot_g(1,i,jj,k)**2.+&
|
|
ub_dot_g(2,i,jj,k)**2.+ub_dot_g(3,i,jj,k)**2.) ! Sum [tke_b_g]
|
|
|
|
RMS_u_g(1:3,i)=RMS_u_g(1:3,i)+(1.-c(i,j,k))*uu_dot_g(1:3,i,jj,k)**2. ! Sum [(1-c)*[(u_u)'^2]_g]
|
|
RMS_u_g(4,i)=RMS_u_g(4,i)+(1.-c(i,j,k))*0.5*(uu_dot_g(1,i,jj,k)**2.+&
|
|
uu_dot_g(2,i,jj,k)**2.+uu_dot_g(3,i,jj,k)**2.) ! Sum [tke_u_g]
|
|
|
|
! COV_g(1,i)=COV_g(1,i)+u_dot_g(1,i,jj,k)*c_dot_g(i,jj,k) ! Sum [u'c']_g
|
|
! COV_g(2,i)=COV_g(2,i)+u_dot_g(2,i,jj,k)*c_dot_g(i,jj,k) ! Sum [v'c']_g
|
|
! COV_g(3,i)=COV_g(3,i)+u_dot_g(3,i,jj,k)*c_dot_g(i,jj,k) ! Sum [w'c']_g
|
|
! ENDIF
|
|
ENDDO ! z -loop
|
|
ENDDO ! y -loop
|
|
ENDDO ! x -loop
|
|
|
|
DO i=1,nxp
|
|
IRE5(1,i)=IRE5(1,i)+RMS(1,i) ! Sum [u'^2]
|
|
IRE5(2,i)=IRE5(2,i)+RMS(2,i) ! Sum [v'^2]
|
|
IRE5(3,i)=IRE5(3,i)+RMS(3,i) ! Sum [w'^2]
|
|
IRE5(4,i)=IRE5(4,i)+RMS(4,i) ! Sum [tke]
|
|
|
|
IRE5(5,i) =IRE5(5,i)+RMS_b(1,i) ! Sum [c*(u_b)^2]
|
|
IRE5(6,i) =IRE5(6,i)+RMS_b(2,i) ! Sum [c*(v_b)^2]
|
|
IRE5(7,i) =IRE5(7,i)+RMS_b(3,i) ! Sum [c*(w_b)^2]
|
|
IRE5(8,i) =IRE5(8,i)+RMS_b(4,i) ! Sum [tke_b]
|
|
|
|
IRE5(9,i) =IRE5(9,i) +RMS_u(1,i) ! Sum [(1-c)*(u_u)^2]
|
|
IRE5(10,i)=IRE5(10,i)+RMS_u(2,i) ! Sum [(1-c)*(v_u)^2]
|
|
IRE5(11,i)=IRE5(11,i)+RMS_u(3,i) ! Sum [(1-c)*(w_u)^2]
|
|
IRE5(12,i)=IRE5(12,i)+RMS_u(4,i) ! Sum [tke_u]
|
|
|
|
! IRE5(13,i)=IRE5(13,i)+COV(1,i) ! Sum [u'c']
|
|
! IRE5(14,i)=IRE5(14,i)+COV(2,i) ! Sum [v'c']
|
|
! IRE5(15,i)=IRE5(15,i)+COV(3,i) ! Sum [w'c']
|
|
! IRE5(16,i)=IRE5(16,i)+COV(4,i) ! Sum [u'(FSD')'] !dhkim
|
|
!
|
|
IRE5(14,i)=IRE5(14,i)+RMS(6,i) ! Sum [u'^2]
|
|
IRE5(15,i)=IRE5(15,i)+RMS(7,i) ! Sum [v'^2]
|
|
IRE5(16,i)=IRE5(16,i)+RMS(8,i) ! Sum [w'^2]
|
|
|
|
IRE5(17,i)=IRE5(17,i)+RMS(5,i) ! Sum [1/3(*u'^2 + v'^2 + w'^2)]
|
|
|
|
IRE5(18,i) =IRE5(18,i)+RMS_b_2(1,i) ! Sum [c*(u'^2)]
|
|
IRE5(19,i) =IRE5(19,i)+RMS_b_2(2,i) ! Sum [c*(v'^2)]
|
|
IRE5(20,i) =IRE5(20,i)+RMS_b_2(3,i) ! Sum [c*(w'^2)]
|
|
IRE5(21,i) =IRE5(21,i)+RMS_b_2(4,i) ! Sum [c*(1/3*[u'^2 + v'^2 + w'^2])]
|
|
|
|
IRE5(22,i) =IRE5(22,i)+RMS_u_2(1,i) ! Sum [(1-c)*(u'^2)]
|
|
IRE5(23,i) =IRE5(23,i)+RMS_u_2(2,i) ! Sum [(1-c)*(v'^2)]
|
|
IRE5(24,i) =IRE5(24,i)+RMS_u_2(3,i) ! Sum [(1-c)*(w'^2)]
|
|
IRE5(25,i) =IRE5(25,i)+RMS_u_2(4,i) ! Sum [(1-c)*(1/3*[u'^2 + v'^2 + w'^2])]
|
|
|
|
|
|
IRE6(1,i)=IRE6(1,i)+RMS_g(1,i) ! Sum [u'^2_g]
|
|
IRE6(2,i)=IRE6(2,i)+RMS_g(2,i) ! Sum [v'^2_g]
|
|
IRE6(3,i)=IRE6(3,i)+RMS_g(3,i) ! Sum [w'^2_g]
|
|
IRE6(4,i)=IRE6(4,i)+RMS_g(4,i) ! Sum [tke_g]
|
|
|
|
IRE6(5,i) =IRE6(5,i)+RMS_b_g(1,i) ! Sum [c*(u_b)^2_g]
|
|
IRE6(6,i) =IRE6(6,i)+RMS_b_g(2,i) ! Sum [c*(v_b)^2_g]
|
|
IRE6(7,i) =IRE6(7,i)+RMS_b_g(3,i) ! Sum [c*(w_b)^2_G]
|
|
IRE6(8,i) =IRE6(8,i)+RMS_b_g(4,i) ! Sum [tke_b_g]
|
|
|
|
IRE6(9,i) =IRE6(9,i) +RMS_u_g(1,i) ! Sum [(1-c)*(u_u)^2_g]
|
|
IRE6(10,i)=IRE6(10,i)+RMS_u_g(2,i) ! Sum [(1-c)*(v_u)^2_g]
|
|
IRE6(11,i)=IRE6(11,i)+RMS_u_g(3,i) ! Sum [(1-c)*(w_u)^2_g]
|
|
IRE6(12,i)=IRE6(12,i)+RMS_u_g(4,i) ! Sum [tke_u_g]
|
|
|
|
! IRE6(13,i)=IRE6(13,i)+COV_g(1,i) ! Sum [u'c']_g
|
|
! IRE6(14,i)=IRE6(14,i)+COV_g(2,i) ! Sum [v'c']_g
|
|
! IRE6(15,i)=IRE6(15,i)+COV_g(3,i) ! Sum [w'c']_g
|
|
ENDDO
|
|
|
|
END SUBROUTINE SAVE_SUM_FLUCTUATION
|
|
|
|
SUBROUTINE FLUCTUATION_AVG
|
|
REAL, DIMENSION(1,nxp) :: ux,dux
|
|
INTEGER :: i
|
|
REAL :: ndata
|
|
WRITE(*,*) 'countnum_FLUC_AVG',countnum
|
|
|
|
ndata=REAL(countnum*(eyp-syp+1)*nzp)
|
|
|
|
! For rms quantities
|
|
|
|
IRE5(1:3,:)=SQRT(IRE5(1:3,:)/ndata)
|
|
IRE5(4,:)=IRE5(4,:)/ndata
|
|
!--------------------uPrime--------------------------------------------------------------------------------
|
|
IRE5(14:16,:)=SQRT(IRE5(14:16,:)/ndata) !RMS(ux`),RMS(uy`),RMS(uz`)
|
|
IRE5(17,:)=SQRT(IRE5(17,:)/ndata) ! [RMS(U`)]
|
|
!IRE5(17,:)=IRE5(17,:)/ndata ! [<U`>]
|
|
!====================uPrime================================================================================
|
|
|
|
DO i=1,nxp
|
|
IRE5(5:7,i)=SQRT(IRE5(5:7,i)/ndata/IRE1(1,i))
|
|
IRE5(8,i)=IRE5(8,i)/ndata/IRE1(1,i)
|
|
IRE5(9:11,i)=SQRT(IRE5(9:11,i)/ndata/(1.-IRE1(1,i)))
|
|
IRE5(12,i)=IRE5(12,i)/ndata/(1.-IRE1(1,i))
|
|
IF(IRE1(1,i).le.0.) IRE5(5:8,i)=0.
|
|
IF(IRE1(1,i).eq.1.) IRE5(9:12,i)=0.
|
|
|
|
!--------------------uPrime--------------------------------------------------------------------------------
|
|
IRE5(18:20,i)=SQRT(IRE5(18:20,i)/ndata/IRE1(1,i)) !RMS(ux`)_b,RMS(uy`)_b,RMS(uz`)_b
|
|
IRE5(21,i)=SQRT(IRE5(21,i)/ndata/IRE1(1,i)) ! [RMS(U`)_b]
|
|
IF(IRE1(1,i).le.0.) IRE5(18:21,i)=0.
|
|
|
|
IRE5(22:24,i)=SQRT(IRE5(22:24,i)/ndata/(1.-IRE1(1,i))) !RMS(ux`)_u,RMS(uy`)_u,RMS(uz`)_u
|
|
IRE5(25,i)=SQRT(IRE5(25,i)/ndata/(1.-IRE1(1,i))) ! [RMS(U`)_u]
|
|
IF(IRE1(1,i).eq.1.) IRE5(22:25,i)=0.
|
|
|
|
! IRE5(18:20,i)=IRE5(18:20,i)/ndata/IRE1(1,i)
|
|
! IRE5(21,i)=IRE5(21,i)/ndata/IRE1(1,i) ! [<U`>b]
|
|
! IF(IRE1(1,i).le.0.) IRE5(18:21,i)=0.
|
|
!
|
|
! IRE5(22:24,i)=IRE5(22:24,i)/ndata/(1.-IRE1(1,i))
|
|
! IRE5(25,i)=IRE5(25,i)/ndata/(1.-IRE1(1,i)) ! [<U`>u]
|
|
! IF(IRE1(1,i).eq.1.) IRE5(22:25,i)=0.
|
|
!====================uPrime================================================================================
|
|
ENDDO
|
|
|
|
|
|
|
|
! IRE5(17,i)=IRE5(17,i)+RMS(5,i) ! Sum [turbulent intensity]
|
|
!
|
|
! IRE5(18,i) =IRE5(18,i)+RMS_b_2(1,i) ! Sum [c*u`^2]
|
|
! IRE5(19,i) =IRE5(19,i)+RMS_b_2(2,i) ! Sum [c*v`^2]
|
|
! IRE5(20,i) =IRE5(20,i)+RMS_b_2(3,i) ! Sum [c*w`^2]
|
|
! IRE5(21,i) =IRE5(21,i)+RMS_b_2(4,i) ! Sum [turbulent intensity_b]
|
|
!
|
|
! IRE5(22,i) =IRE5(22,i)+RMS_u_2(1,i) ! Sum [(1-c)*u`^2]
|
|
! IRE5(23,i) =IRE5(23,i)+RMS_u_2(2,i) ! Sum [(1-c)*v`^2]
|
|
! IRE5(24,i) =IRE5(24,i)+RMS_u_2(3,i) ! Sum [(1-c)*w`^2]
|
|
! IRE5(25,i) =IRE5(25,i)+RMS_u_2(4,i) ! Sum [turbulent intensity_u]
|
|
|
|
|
|
|
|
|
|
IRE6(1:3,:)=SQRT(IRE6(1:3,:)/ndata)
|
|
IRE6(4,:)=IRE6(4,:)/ndata
|
|
|
|
DO i=1,nxp
|
|
IRE6(5:7,i)=SQRT(IRE6(5:7,i)/ndata/IRE1(1,i))
|
|
IRE6(8,i)=IRE6(8,i)/ndata/IRE1(1,i)
|
|
IRE6(9:11,i)=SQRT(IRE6(9:11,i)/ndata/(1.-IRE1(1,i)))
|
|
IRE6(12,i)=IRE6(12,i)/ndata/(1.-IRE1(1,i))
|
|
IF(IRE1(1,i).le.0.) IRE6(5:8,i)=0.
|
|
IF(IRE1(1,i).eq.1.) IRE6(9:12,i)=0.
|
|
ENDDO
|
|
|
|
IRE6(13:15,:)=IRE6(13:15,:)/ndata
|
|
|
|
!! For ST
|
|
! DO i=1,nxp
|
|
! ux(1,i)=IRE1(2,i) ! <FSD>
|
|
! ENDDO
|
|
! CALL dfnonp(nxp,hxp,ux(1,:),dux(1,:),1,1)
|
|
!
|
|
! DO i=1,nxp
|
|
! IRE8(1,i)=IRE1(11,i)/IRE1(1,i) ! (1/<c>)*(d<c>/dx) = 1/Lw1
|
|
! IF(IRE1(1,i).le.0.) IRE8(1,i)=0.
|
|
! IRE8(2,i)=dux(1,i)/IRE1(2,i) ! (1/<FSD>)*(d<FSD>/dx)= 1/Lw2
|
|
! IF(IRE1(2,i).eq.0.) IRE8(2,i)=0.
|
|
! IRE8(3,i)= -IRE5(13,i)/IRE1(11,i) ! Dtu
|
|
! IRE8(4,i)= -IRE6(13,i)/IRE1(11,i) ! Dtu_g
|
|
! IF(IRE1(11,i).eq.0.) IRE8(3:4,i)=0.
|
|
! IRE8(5,i) = SQRT(IRE1(10,i)/(IRE1(10,i)+IRE8(3,i)))*SL_u/IRE1(10,i) !SL_u/Dmu * SQRT[ Dmu/(Dmu+Dtu) ] = 1/Lw_3
|
|
! if(IRE1(10,i)/(IRE1(10,i)+IRE8(3,i)).lt.0.) IRE8(5,i)=0.
|
|
!
|
|
! IRE8(6,i) = (IRE1(10,i)+IRE8(3,i))*IRE8(1,i) ! ST1
|
|
! IRE8(7,i) = (IRE1(10,i)+IRE8(3,i))*IRE8(2,i) ! ST2
|
|
! IRE8(8,i) = (IRE1(10,i)+IRE8(3,i))*(IRE2(11,i)-IRE2(9,i)) ! ST3
|
|
! IRE8(9,i) = (IRE1(10,i)+IRE8(3,i))*(IRE2(12,i)-IRE2(9,i)) ! ST4
|
|
! IRE8(10,i) = SL_u * SQRT(1.+IRE8(3,i)/IRE1(10,i)) ! ST5
|
|
! if((1.+IRE8(3,i)/IRE1(10,i)).lt.0.) IRE8(10,i)=0.
|
|
!
|
|
! IRE8(11,i) = 1./IRE2(11,i) ! Lm*_x
|
|
! if(IRE2(11,i).eq.0.) IRE8(11,i)=0.
|
|
!
|
|
! IRE8(12,i)= 1./IRE2(12,i) ! Lm*_n
|
|
! if(IRE2(12,i).eq.0.) IRE8(12,i)=0.
|
|
!
|
|
! IRE8(13,i)= 1./IRE8(1,i) ! Lw
|
|
! if(IRE8(1,i).eq.0.) IRE8(13,i)=0.
|
|
!
|
|
! IRE8(14,i) = 1./IRE8(5,i) ! Lw_3
|
|
! if(IRE8(5,i).le.0.) IRE8(14,i)=0.
|
|
!
|
|
! IRE8(15,i)=dux(1,i) ! d<FSD>/dx !dhkim
|
|
! IRE8(16,i)= -IRE5(16,i)/IRE8(15,i) ! Dts !dhkim
|
|
! if(IRE8(15,i).eq.0.) IRE8(16,i)=0.
|
|
! IRE8(17,i)= IRE1(20,i)/(1.-IRE1(1,i)) ! 1/(1-c)*d(1-c)/dx !dhkim
|
|
! if(IRE1(1,i).eq.1.) IRE8(17,i)=0.
|
|
!
|
|
!
|
|
! IRE8(18,i) = SL_u/IRE1(10,i) * SQRT( IRE1(10,i)/(IRE1(10,i)+IRE8(3,i)) ) ! 1/L_LE_3
|
|
! if(IRE1(10,i).eq.0.) IRE8(18,i)=0.
|
|
! if((IRE1(10,i)/(IRE1(10,i)+IRE8(3,i))).lt.0.) IRE8(18,i)=0.
|
|
! IRE8(19,i) = (IRE2(34,i)-IRE2(9,i))/IRE2(31,i) ! 1/L_LE_4
|
|
! if(IRE2(31,i).eq.0.) IRE8(19,i)=0.
|
|
! IRE8(20,i) = (IRE8(3,i)+IRE1(10,i))*IRE8(18,i) ! ST_6
|
|
! IRE8(21,i) = (IRE8(3,i)+IRE1(10,i))*IRE8(19,i) ! ST_7
|
|
! IRE8(22,i) = (IRE8(3,i)+IRE1(10,i))*((IRE2(11,i)-IRE2(9,i))/IRE2(31,i)) ! ST_8 (ST_3/|<nx>f|)
|
|
! if(IRE2(31,i).eq.0.) IRE8(22,i)=0.
|
|
!
|
|
! IRE8(23,i) = SL_u/IRE1(10,i) - IRE2(9,i) ! 1/L_LE_5=1/Lm-<DivN>f
|
|
! IRE8(24,i) = (IRE8(3,i)+IRE1(10,i))*IRE8(23,i) ! ST_9
|
|
! IRE8(25,i) = IRE1(24,i)/IRE1(11,i) ! 1/L_LE_6=1/(d<c>/dx)*(d2<c>/dx2)
|
|
! if(IRE1(11,i).eq.0.) IRE8(25,i)=0.
|
|
! IRE8(26,i) = (IRE8(3,i)+IRE1(10,i))*IRE8(25,i) ! ST_10
|
|
! ENDDO
|
|
|
|
END SUBROUTINE FLUCTUATION_AVG
|
|
|
|
! SUBROUTINE FINAL_AVG
|
|
! REAL, DIMENSION(4,nxp) :: ux,dux
|
|
! INTEGER :: i
|
|
! DO i=1,nxp
|
|
! IRE9(1,i) = IRE1(5,i)*IRE1(11,i) +IRE1(6,i)*IRE1(12,i) +IRE1(7,i)*IRE1(13,i)
|
|
! ! <v> dot Grad.<c>
|
|
! IRE9(2,i) = IRE8(16,i)*IRE1(23,i) ! Dts*Lap.<c>
|
|
! IRE9(3,i) = IRE2(2,i)*IRE1(2,i) ! <Sd>f*<FSD>
|
|
! IRE9(4,i) = IRE8(16,i)*IRE1(2,i)*IRE2(33,i) ! Dts*<FSD>*Div.(<n>f)
|
|
! IRE9(5,i) = IRE2(28,i)*IRE1(2,i) ! <V`` dot N``>f*<FSD>
|
|
! IRE9(6,i) = IRE9(1,i)-IRE9(2,i)-IRE9(3,i)-IRE9(4,i)-IRE9(5,i) ! cEqn.Balance
|
|
!
|
|
! IRE9(12,i) = IRE1(20,i)/(1-IRE1(1,i)) ! 1/(1-<c>)*(d(1-<c>)/dx)
|
|
! if(IRE1(1,i).eq.1.) IRE9(12,i)=0. ! = 1/L_TE
|
|
!
|
|
! ENDDO
|
|
!
|
|
! DO i=1,nxp
|
|
! ux(1,i)=IRE8(3,i) ! <Dt>
|
|
! ux(2,i)=IRE8(3,i)+IRE1(10,i) ! <Dt>+<Dm>
|
|
! ux(3,i)=1/IRE8(1,i) ! L_LE
|
|
! if(IRE8(1,i).eq.0.) ux(3,i)=0. ! = 1/L_TE
|
|
! ux(4,i)=1/IRE9(12,i) ! L_TE
|
|
! if(IRE9(12,i).eq.0.) ux(4,i)=0. ! = 1/L_TE
|
|
! ENDDO
|
|
! CALL dfnonp(nxp,hxp,ux(1:4,:),dux(1:4,:),4,1)
|
|
!
|
|
! DO i=1,nxp
|
|
! IRE9(7,i) = IRE1(27,i)/IRE1(19,i) ! <Nx>K
|
|
! if(IRE1(19,i).eq.0.) IRE9(7,i)=0.
|
|
! IRE9(8,i) = IRE2(16,i)/IRE9(7,i) ! <Nx>f/<Nx>K
|
|
! if(IRE9(7,i).eq.0.) IRE9(8,i)=0.
|
|
! IRE9(9,i) = dux(1,i) ! d<Dt>/dx
|
|
! IRE9(10,i) = dux(2,i) ! d<Dm+Dt>/dx
|
|
! IRE9(11,i) = IRE1(28,i)/IRE1(20,i) ! 1/(d<1-c>/dx)*(d2<1-c>/dx2)
|
|
! ! IRE9(12,i) is located above Do loop.
|
|
! !IRE9(12,i) = IRE1(20,i)/(1-IRE1(1,i)) ! 1/(1-<c>)*(d(1-<c>)/dx)
|
|
! ! = 1/L_TE
|
|
! IRE9(13,i) = dux(3,i) ! d(L_LE)/dx
|
|
! if(ux(3,i).eq.0.) dux(3,i)=0. ! = 1/L_TE
|
|
! IRE9(14,i) = dux(4,i) ! d(L_TE)/dx
|
|
! if(ux(4,i).eq.0.) dux(4,i)=0. ! = 1/L_TE
|
|
! ENDDO
|
|
! END SUBROUTINE FINAL_AVG
|
|
|
|
|
|
SUBROUTINE SAVE_AVG_RESULTS
|
|
INTEGER :: i
|
|
OPEN (200,FILE="qEdge_X.dat")
|
|
! IRE1
|
|
WRITE(200,*) 'VARIABLES = "X","<C>","<FSD>","<Yr>","<Wc/rho>","<u>","<v>","<w>"' ! 8
|
|
!WRITE(200,*) 'VARIABLES = "X","<C>","<FSD>","<Yr>","<Wc/rho>","<u>","<v>","<w>","<rho>"' ! 9
|
|
! WRITE(200,*) '"<T>","<D>","<dc/dx>","<dc/dy>","<dc/dz>","<du/dx>","<dv/dy>","<dw/dz>"' ! 8 -> 17
|
|
! WRITE(200,*) '"<Div(rho*V)>","<c>_g","<d2c/dn2>"' ! 3 -> 20
|
|
! WRITE(200,*) '"d(1-<c>)/dx","d(1-<c>)/dy","d(1-<c>)/dz"' ! 3 -> 23
|
|
! WRITE(200,*) '"<Lap.(c)>","<d2c/dx2>","<d2c/dy2>","<d2c/dz2>"' ! 4 -> 27
|
|
! WRITE(200,*) '"<Nx*(d2c/dn2)>","<d2(1-c)/dx2>"' ! 2 -> 29
|
|
|
|
! IRE2
|
|
! WRITE(200,*) '"<FSD2>","<Sd>f","<Sdr>f","<Sdd>f","<Sdn>f","<Sdt>f","<Sdd2>f","<Sd2>f"' ! 8 -> 8
|
|
! WRITE(200,*) '"<DivN>f","<|DivN|>f","<(dc/dx)/c>f"' ! 3 -> 11
|
|
! WRITE(200,*) '"<-(dc/dn)/c>f","<(1/FSD`)*d(FSD`)/dx>f","<(1/FSD`)*d(FSD`)/dy>f"' ! 3 -> 14
|
|
! WRITE(200,*) '"<(1/FSD`)*d(FSD`)/dz>f","<Nx>f","<Ny>f","<Nz>f","<u>f","<v>f","<w>f"' ! 7 -> 21
|
|
! WRITE(200,*) '"<vn>f","<N dot N>f","<d(vn)/dn>f","<d(Sd)/dn>f","<d(Vn+Sd)/dn>f"' ! 5 -> 26
|
|
! WRITE(200,*) '"<v>f dot <N>f","<v``dot N``>f","<N>f dot <N>f","<N``dot N``>f"' ! 4 -> 30
|
|
! WRITE(200,*) '"|<N>f|","<vn+Sd>f","d<Nx>f/dx"' ! 3 -> 33
|
|
! WRITE(200,*) '"-<N dot (grad.(FSD`))/FSD`>f"' ! 1 -> 34
|
|
|
|
! IRE3
|
|
WRITE(200,*) '"<u>b","<v>b","<w>b","<rho>b","<u>b_g","<v>b_g","<w>b_g","<rho>b_g"' ! 8 -> 8
|
|
! IRE4
|
|
WRITE(200,*) '"<u>u","<v>u","<w>u","<rho>u","<u>u_g","<v>u_g","<w>u_g","<rho>u_g"' ! 8 -> 8
|
|
! IRE5
|
|
WRITE(200,*) '"RMS(u`)","RMS(v`)","RMS(w`)","<k>"' ! 4 -> 4
|
|
WRITE(200,*) '"RMS(u`)b","RMS(v`)b","RMS(w`)b","<k>b"' ! 4 -> 8
|
|
WRITE(200,*) '"RMS(u`)u","RMS(v`)u","RMS(w`)u","<k>u"' ! 4 -> 12
|
|
!WRITE(200,*) '"<u`c`>","<v`c`>","<w`c`>"' ! 3 -> 15
|
|
WRITE(200,*) '"<u`c`>","RMS(ux`)","RMS(uy`)"' ! 3 -> 15
|
|
WRITE(200,*) '"RMS(uz`)"' ! 1 -> 16
|
|
WRITE(200,*) '"RMS( U`)","RMS(ux`)_b","RMS(uy`)_b","RMS(uz`)_b"' ! 4 -> 20
|
|
WRITE(200,*) '"RMS(U`)_b","RMS(ux`)_u","RMS(uy`)_u","RMS(uz`)_u"' ! 4 -> 24
|
|
WRITE(200,*) '"RMS(U`)_u"' ! 1 -> 25
|
|
!
|
|
! IRE6
|
|
WRITE(200,*) '"RMS(u`)_g","RMS(v`)_g","RMS(w`)_g","<k>_g"' ! 4 -> 4
|
|
WRITE(200,*) '"RMS(u`)b_g","RMS(v`)b_g","RMS(w`)b_g","<k>b_g"' ! 4 -> 8
|
|
WRITE(200,*) '"RMS(u`)u_g","RMS(v`)u_g","RMS(w`)u_g","<k>u_g"' ! 4 -> 12
|
|
WRITE(200,*) '"<u`c`>_g","<v`c`>_g","<w`c`>_g"' ! 3 -> 15
|
|
!! IRE7
|
|
! WRITE(200,*) '"<vn>k<d2c/dn2>","<sd>k<d2c/dn2>","<vn+sd>k<d2c/dn2>"' ! 3 -> 3
|
|
! WRITE(200,*) '"<vn>k","<sd>k","<vn+sd>k"' ! 3 -> 6
|
|
!! IRE8
|
|
! WRITE(200,*) '"(1/<c>)*(d<c>/dx)","(1/<FSD>)*(d<FSD>/dx)","Dt_x","Dt_x_g"' ! 4 -> 4
|
|
! WRITE(200,*) '"1/Lw_3_High_Turb"' ! 1 -> 5
|
|
! WRITE(200,*) '"ST1","ST2","ST3","ST4","ST5","Lm*_x","Lm*_n","Lw","Lw_3"' ! 9 -> 14
|
|
! WRITE(200,*) '"d<FSD>/dx","Dts"' ! 2 -> 16
|
|
! WRITE(200,*) '"(1/(1-<c>))*(d(1-<c>)/dx)"' ! 1 -> 17
|
|
! WRITE(200,*) '"1/L_LE_3","1/L_LE_4","ST6","ST7","ST8"' ! 5 -> 22
|
|
! WRITE(200,*) '"1/L_LE_5=1/Lm-<DivN>f","ST_9"' ! 2 -> 24
|
|
! WRITE(200,*) '"1/L_LE_6=1/(d<c>/dx)*(d2<c>/dx2)","ST_10"' ! 2 -> 26
|
|
!
|
|
!! IRE9
|
|
! WRITE(200,*) '"<v>dotGrad.<c>","Dts*Lap.<c>","<Sd>f*<FSD>"' ! 3 -> 3
|
|
! WRITE(200,*) '"Dts*<FSD>*Div.(<n>f)","<V`` dot N``>f*<FSD>","cEqnBalance"' ! 3 -> 6
|
|
! WRITE(200,*) '"<Nx>K","<Nx>f/<Nx>K","d<Dt>/dx"' ! 3 -> 9
|
|
! WRITE(200,*) '"d<Dm+Dt>/dx","1/(d<1-c>/dx)*(d2<1-c>/dx2)"' ! 2 -> 11
|
|
! WRITE(200,*) '"1/(1-<c>)*(d(1-<c>)/dx)","d(L_LE)/dx","d(L_TE)/dx"' ! 3 -> 14
|
|
|
|
|
|
DO i=1,nxp
|
|
! WRITE(200,'(156e20.10)') REAL(i)*hxp,IRE1(1:28,i),IRE2(1:34,i),IRE3(1:8,i),& ! 1+28+34+8 = 71
|
|
! IRE4(1:8,i),IRE5(1:16,i),IRE6(1:15,i),IRE7(1:6,i),IRE8(1:26,i),IRE9(1:14,i)
|
|
! ! 8+16+15+6+26+14 = 85 -> 156
|
|
WRITE(200,'(64e20.10)') REAL(i)*hxp,IRE1(1:7,i),IRE3(1:8,i),& ! 1+7+8 =16
|
|
IRE4(1:8,i),IRE5(1:25,i),IRE6(1:15,i)
|
|
! 8+25+15=48 -> 64
|
|
ENDDO
|
|
|
|
|
|
CLOSE(200)
|
|
!CLOSE(210)
|
|
|
|
END SUBROUTINE SAVE_AVG_RESULTS
|
|
|
|
SUBROUTINE ALLOCATE_ARRAYS
|
|
|
|
INTEGER :: ierr
|
|
|
|
ALLOCATE(u(nxp,nyp,nzp),STAT=ierr) ; u=0. ! Main variables
|
|
ALLOCATE(v(nxp,nyp,nzp),STAT=ierr) ; v=0. ! Main variables
|
|
ALLOCATE(w(nxp,nyp,nzp),STAT=ierr) ; w=0. ! Main variables
|
|
|
|
ALLOCATE(old_scalar(2,nxp,nyp,nzp),STAT=ierr) ; old_scalar=0. ! Main variables
|
|
ALLOCATE(new_scalar(nxp,nyp,nzp,2),STAT=ierr) ; new_scalar=0. ! Main variables
|
|
|
|
ALLOCATE(m_v(6,nxp,nyp,nzp),STAT=ierr) ; m_v=0. ! Main variables
|
|
ALLOCATE(m_v_new(nxp,nyp,nzp,2),STAT=ierr) ; m_v_new=0. ! Main variables
|
|
!ALLOCATE(G_C(3,nxp,nyp,nzp),STAT=ierr) ; G_C=0.
|
|
!ALLOCATE(cgm(nxp,nyp,nzp),STAT=ierr) ; cgm=0.
|
|
ALLOCATE(NV(3,nxp,nyp,nzp),STAT=ierr) ; NV=0.
|
|
!ALLOCATE(sd(7,nxp,nyp,nzp),STAT=ierr) ; sd=0.
|
|
ALLOCATE(c(nxp,nyp,nzp),STAT=ierr) ; c=0.
|
|
ALLOCATE(Wc(nxp,nyp,nzp),STAT=ierr) ; Wc=0.
|
|
ALLOCATE(u_dot(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; u_dot=0.
|
|
ALLOCATE(y(nxp,nyp,nzp),STAT=ierr) ; y=0.
|
|
!ALLOCATE(DivN(nxp,nyp,nzp),STAT=ierr) ; DivN=0.
|
|
!ALLOCATE(G_V(3,nxp,nyp,nzp),STAT=ierr) ; G_V=0.
|
|
!ALLOCATE(Div_V(nxp,nyp,nzp),STAT=ierr) ; Div_V=0
|
|
!ALLOCATE(G_FSD(3,nxp,nyp,nzp),STAT=ierr) ; G_FSD=0.
|
|
!ALLOCATE(vn(nxp,nyp,nzp),STAT=ierr) ; vn=0.
|
|
!ALLOCATE(G_vn(3,nxp,nyp,nzp),STAT=ierr) ;G_vn=0.
|
|
!ALLOCATE(G_sd(3,nxp,nyp,nzp),STAT=ierr) ;G_sd=0.
|
|
!ALLOCATE(G_Y(3,nxp,nyp,nzp),STAT=ierr) ; G_Y=0.
|
|
!ALLOCATE(GN_vn(nxp,nyp,nzp),STAT=ierr) ;GN_vn=0.
|
|
!ALLOCATE(GN_sd(nxp,nyp,nzp),STAT=ierr) ;GN_sd=0.
|
|
!ALLOCATE(G2N_c(nxp,nyp,nzp),STAT=ierr) ;G2N_c=0.
|
|
!ALLOCATE(G2_Y(nxp,nyp,nzp),STAT=ierr) ;G2_Y=0.
|
|
!ALLOCATE(d2gc(4,nxp,nyp,nzp),STAT=ierr) ;d2gc=0.
|
|
|
|
ALLOCATE(ub_dot(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; ub_dot=0.
|
|
ALLOCATE(uu_dot(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; uu_dot=0.
|
|
ALLOCATE(c_dot(nxp,(eyp-syp+1),nzp),STAT=ierr) ; c_dot=0.
|
|
ALLOCATE(FSD_dot(nxp,(eyp-syp+1),nzp),STAT=ierr) ; FSD_dot=0. !dhkim
|
|
ALLOCATE(ub_dot_g(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; ub_dot_g=0.
|
|
ALLOCATE(uu_dot_g(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; uu_dot_g=0.
|
|
ALLOCATE(c_dot_g(nxp,(eyp-syp+1),nzp),STAT=ierr) ; c_dot_g=0.
|
|
ALLOCATE(c_g(nxp,(eyp-syp+1),nzp),STAT=ierr) ; c_g=0.
|
|
|
|
! Arrays for local sum
|
|
ALLOCATE(u_g(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; u_g=0.
|
|
ALLOCATE(u_dot_g(3,nxp,(eyp-syp+1),nzp),STAT=ierr) ; u_dot_g=0.
|
|
ALLOCATE(RMS_g(4,nxp),STAT=ierr) ; RMS_g=0.
|
|
ALLOCATE(RMS_b_g(4,nxp),STAT=ierr) ; RMS_b_g=0.
|
|
ALLOCATE(RMS_u_g(4,nxp),STAT=ierr) ; RMS_u_g=0.
|
|
ALLOCATE(CM_b_g(4,nxp,(eyp-syp+1),nzp),STAT=ierr) ; CM_b_g=0.
|
|
ALLOCATE(CM_u_g(4,nxp,(eyp-syp+1),nzp),STAT=ierr) ; CM_u_g=0.
|
|
ALLOCATE(COV_g(3,nxp),STAT=ierr) ; COV_g=0.
|
|
|
|
! Arrays for total sum
|
|
!ALLOCATE(SMF(23,nxp),STAT=ierr) ; SMF=0.
|
|
ALLOCATE(SMF(3,nxp),STAT=ierr) ; SMF=0.
|
|
ALLOCATE(SMC(4,nxp),STAT=ierr) ; SMC=0.
|
|
!ALLOCATE(FMS(26,nxp),STAT=ierr) ; FMS=0.
|
|
!ALLOCATE(KMS(2,nxp),STAT=ierr) ; KMS=0.
|
|
|
|
ALLOCATE(CM_b(4,nxp),STAT=ierr) ; CM_b=0.
|
|
ALLOCATE(CM_u(4,nxp),STAT=ierr) ; CM_u=0.
|
|
ALLOCATE(RMS(8,nxp),STAT=ierr) ; RMS=0.
|
|
ALLOCATE(RMS_b(4,nxp),STAT=ierr) ; RMS_b=0.
|
|
ALLOCATE(RMS_u(4,nxp),STAT=ierr) ; RMS_u=0.
|
|
ALLOCATE(COV(4,nxp),STAT=ierr) ; COV=0.
|
|
|
|
ALLOCATE(favg_ndata(nxp),STAT=ierr) ; favg_ndata=0.
|
|
ALLOCATE(RMS_u_2(4,nxp),STAT=ierr) ; RMS_u_2=0.
|
|
ALLOCATE(RMS_b_2(4,nxp),STAT=ierr) ; RMS_b_2=0.
|
|
|
|
! Arrays for final averages
|
|
!ALLOCATE(IRE1(28,nxp),STAT=ierr) ; IRE1=0.
|
|
ALLOCATE(IRE1(7,nxp),STAT=ierr) ; IRE1=0.
|
|
!ALLOCATE(IRE2(34,nxp),STAT=ierr) ; IRE2=0.
|
|
ALLOCATE(IRE3(8,nxp),STAT=ierr) ; IRE3=0.
|
|
ALLOCATE(IRE4(8,nxp),STAT=ierr) ; IRE4=0.
|
|
ALLOCATE(IRE5(25,nxp),STAT=ierr) ; IRE5=0.
|
|
ALLOCATE(IRE6(15,nxp),STAT=ierr) ; IRE6=0.
|
|
!ALLOCATE(IRE7(6,nxp),STAT=ierr) ; IRE7=0.
|
|
!ALLOCATE(IRE8(26,nxp),STAT=ierr) ; IRE8=0.
|
|
!ALLOCATE(IRE9(14,nxp),STAT=ierr) ; IRE9=0.
|
|
|
|
|
|
! hyp=l_0*pi/REAL(nyp)
|
|
hyp=l_0*pi/REAL(nyp-1) ! kwon
|
|
hxp=hyp; hzp=hyp
|
|
|
|
WRITE(*,'(a6,i3,a8,i3,a8,i3)') ' NX = ',nxp,' / NY = ',nyp,' / NZ = ',nzp
|
|
WRITE(*,*) ' Preparing memory space for COMPACT SCHEME'
|
|
CALL ludcmp(nxp,nyp,nzp,1,0,0) ! 1,1,0
|
|
|
|
WRITE(*,'(a22,i3,a3,i3,a4,i3)') ' Grid number range : ',syp,' ~ ',eyp,' of ',nyp
|
|
WRITE(*,*)
|
|
|
|
END SUBROUTINE ALLOCATE_ARRAYS
|
|
|
|
SUBROUTINE DEALLOCATES_CLOSE
|
|
DEALLOCATE(u)
|
|
DEALLOCATE(v)
|
|
DEALLOCATE(w)
|
|
DEALLOCATE(old_scalar)
|
|
DEALLOCATE(new_scalar)
|
|
DEALLOCATE(m_v)
|
|
!DEALLOCATE(G_C)
|
|
!DEALLOCATE(cgm)
|
|
DEALLOCATE(NV)
|
|
!DEALLOCATE(sd)
|
|
DEALLOCATE(c)
|
|
DEALLOCATE(Wc)
|
|
DEALLOCATE(u_dot)
|
|
DEALLOCATE(y)
|
|
!DEALLOCATE(DivN)
|
|
! DEALLOCATE(G_V)
|
|
! DEALLOCATE(Div_V)
|
|
! DEALLOCATE(G_FSD)
|
|
! DEALLOCATE(vn)
|
|
! DEALLOCATE(G_vn)
|
|
! DEALLOCATE(G_sd)
|
|
! DEALLOCATE(G_Y)
|
|
! DEALLOCATE(GN_vn)
|
|
! DEALLOCATE(GN_sd)
|
|
! DEALLOCATE(G2N_c)
|
|
! DEALLOCATE(G2_Y)
|
|
! DEALLOCATE(d2gc)
|
|
DEALLOCATE(m_v_new)
|
|
|
|
DEALLOCATE(uu_dot)
|
|
DEALLOCATE(ub_dot)
|
|
DEALLOCATE(c_dot)
|
|
DEALLOCATE(FSD_dot)
|
|
DEALLOCATE(ub_dot_g)
|
|
DEALLOCATE(uu_dot_g)
|
|
DEALLOCATE(c_dot_g)
|
|
DEALLOCATE(c_g)
|
|
|
|
DEALLOCATE(u_g)
|
|
DEALLOCATE(u_dot_g)
|
|
DEALLOCATE(RMS_g)
|
|
DEALLOCATE(RMS_b_g)
|
|
DEALLOCATE(RMS_u_g)
|
|
DEALLOCATE(CM_b_g)
|
|
DEALLOCATE(CM_u_g)
|
|
DEALLOCATE(COV_g)
|
|
|
|
DEALLOCATE(SMF)
|
|
DEALLOCATE(SMC)
|
|
!DEALLOCATE(FMS)
|
|
|
|
DEALLOCATE(CM_b)
|
|
DEALLOCATE(CM_u)
|
|
DEALLOCATE(RMS)
|
|
DEALLOCATE(RMS_b)
|
|
DEALLOCATE(RMS_u)
|
|
DEALLOCATE(COV)
|
|
DEALLOCATE(favg_ndata)
|
|
DEALLOCATE(RMS_b_2)
|
|
DEALLOCATE(RMS_u_2)
|
|
|
|
DEALLOCATE(IRE1)
|
|
!DEALLOCATE(IRE2)
|
|
DEALLOCATE(IRE3)
|
|
DEALLOCATE(IRE4)
|
|
DEALLOCATE(IRE5)
|
|
DEALLOCATE(IRE6)
|
|
!DEALLOCATE(IRE7)
|
|
!DEALLOCATE(IRE8)
|
|
!DEALLOCATE(IRE9)
|
|
|
|
IF(omitnum.gt.0) DEALLOCATE(omit_t)
|
|
END SUBROUTINE DEALLOCATES_CLOSE
|
|
|
|
END MODULE post
|