commit 7b4dec31f5b785126ebb6b60f37baffe0bf92bcf Author: ignis Date: Mon Apr 1 15:12:36 2019 +0900 Initial commit diff --git a/code/Compact.f90 b/code/Compact.f90 new file mode 100644 index 0000000..e75fcf1 --- /dev/null +++ b/code/Compact.f90 @@ -0,0 +1,386 @@ + MODULE Compact + IMPLICIT NONE + PRIVATE + REAL, DIMENSION(:), ALLOCATABLE :: lxf,lxs,wxf,wxs, & + lyf,lys,wyf,wys, & + lzf,lzs,wzf,wzs +! lyzf,lyzs,wyzf,wyzs + INTEGER :: nxc,nyc,nzc + REAL, PARAMETER :: ezero = 1.0e-14 + + PUBLIC :: ludcmp,dfnonp,d2fnonp,dfp,d2fp + + CONTAINS + + SUBROUTINE ludcmp(nx,ny,nz,xp,yp,zp) + INTEGER, INTENT(IN) :: nx,ny,nz + INTEGER, INTENT(IN) :: xp,yp,zp + INTEGER :: ierr + + nxc=nx + nyc=ny + nzc=nz +! IF(nyc /= nzc) PRINT*,'ny should be equal nz' + +! xp, yp, zp = 0 : periodic + ALLOCATE(lxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(xp.eq.0) THEN + ALLOCATE(wxf(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wxs(nxc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(1,nxc) + ELSE + CALL nonp_lud(1,nxc) + ENDIF + + ALLOCATE(lyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(yp.eq.0) THEN + ALLOCATE(wyf(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wys(nyc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(2,nyc) + ELSE + call nonp_lud(2,nyc) + ENDIF + + ALLOCATE(lzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(lzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + IF(zp.eq.0) THEN + ALLOCATE(wzf(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + ALLOCATE(wzs(nzc),STAT=ierr) + IF(ierr /= 0) PRINT*, 'work array for lud allocation failed' + CALL p_lud(3,nzc) + ELSE + call nonp_lud(3,nzc) + ENDIF + +! CALL x_lud +! CALL yz_lud + + END SUBROUTINE ludcmp + + SUBROUTINE nonp_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL, DIMENSION(xx) :: aa + aa=3. + aa(1)=0.5 ; aa(2)=4. + aa(xx-1)=4. ; aa(xx)=0.5 +! first derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxf) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lyf) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzf) ! z-direction + aa=5.5 + aa(1)=2./11. ; aa(2)=10. + aa(xx-1)=10. ; aa(xx)=2./11. +! second derivative + IF (xyz.eq.1) CALL stdlu(aa,xx,lxs) ! x-direction + IF (xyz.eq.2) CALL stdlu(aa,xx,lys) ! y-direction + IF (xyz.eq.3) CALL stdlu(aa,xx,lzs) ! z-direction + END SUBROUTINE nonp_lud + + + SUBROUTINE p_lud(xyz,xx) + INTEGER :: i,xyz,xx + REAL :: a + a=3. ! first derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxf,wxf) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lyf,wyf) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzf,wzf) ! z-direction + a=11./2. ! second derivative + IF (xyz.eq.1) CALL ptdlu(a,xx,lxs,wxs) ! x-direction + IF (xyz.eq.2) CALL ptdlu(a,xx,lys,wys) ! y-direction + IF (xyz.eq.3) CALL ptdlu(a,xx,lzs,wzs) ! z-direction + END SUBROUTINE p_lud + + SUBROUTINE stdlu(a,n,l) + INTEGER :: n + REAL :: a(n),l(n) + REAL :: d + INTEGER :: i + l(1)=1.0/a(1) + DO i=2,n + d=a(i)-l(i-1) + l(i)=1.0/d + ENDDO + END SUBROUTINE stdlu + + SUBROUTINE ptdlu(a,n,l,w) + INTEGER :: n + REAL :: a,l(n),w(n) + INTEGER :: i + REAL :: aa(n),d + + DO i=1,n-1 + aa(i)=a + ENDDO + i=n-1 + call stdlu(aa,i,l) + w(1)=1.0 + DO i=2,n-2 + w(i)=-l(i-1)*w(i-1) + ENDDO + w(n-1)=1.0-l(n-2)*w(n-2) + DO i=1,n-1 + w(i)=w(i)*l(i) + ENDDO + d=a + DO i=1,n-1 + d=d-w(i)*w(i)/l(i) + ENDDO + l(n)=1./d + END SUBROUTINE ptdlu + + SUBROUTINE dfnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL,INTENT(IN) :: h + REAL,INTENT(IN),DIMENSION(nd,n) :: x + REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL :: r1,r2,r3,a,b,c,h1,t1,t2,t3,t4 + + h1=1./h + + r1=7./3. + r2=1./12. + r3=3. + a=-1.25 + b=1. + c=0.25 + + DO j=1,nd + dx(j,n-1)=x(j,n)-x(j,n-2) + dx(j,n)=-(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)) + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)) + dx(j,2)=x(j,3)-x(j,1) + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n)=0. + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3)) dx(j,1)=0. + dx(j,n-1)=dx(j,n-1)*h1*r3 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(j,2)*h1*r3 + ENDDO + + DO i=3,n-2 + DO j=1,nd + t1=x(j,i+1)-x(j,i-1) + t2=x(j,i+2)-x(j,i-2) + dx(j,i)=h1*(r1*t1+r2*t2) + ENDDO + ENDDO + IF (dir.eq.1) CALL tdslv(dx,n,lxf,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lyf,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzf,nd) ! z-direction + END SUBROUTINE dfnonp + + SUBROUTINE dfp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL,INTENT(IN) :: h + REAL,INTENT(IN),DIMENSION(nd,n) :: x + REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL :: r1,r2,h1 + + h1=1./h + r1=7./3. + r2=1./12. + + DO j=1,nd + dx(j,n-1)=(r1*(x(j,n)-x(j,n-2))+r2*(x(j,1)-x(j,n-3))) + dx(j,n)=(r1*(x(j,1)-x(j,n-1))+r2*(x(j,2)-x(j,n-2))) + dx(j,1)=(r1*(x(j,2)-x(j,n))+r2*(x(j,3)-x(j,n-1))) + dx(j,2)=(r1*(x(j,3)-x(j,1))+r2*(x(j,4)-x(j,n))) + dx(j,n-1)=dx(j,n-1)*h1 + dx(j,n)=dx(j,n)*h1 + dx(j,1)=dx(j,1)*h1 + dx(j,2)=dx(j,2)*h1 + ENDDO + + DO i=3,n-2 + DO j=1,nd + dx(j,i)=(r1*(x(j,i+1)-x(j,i-1))+r2*(x(j,i+2)-x(j,i-2))) + dx(j,i)=dx(j,i)*h1 + ENDDO + ENDDO + + IF (dir.eq.1) CALL ptdslv(dx,n,lxf,wxf,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lyf,wyf,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzf,wzf,nd) ! z-direction + + END SUBROUTINE dfp + + SUBROUTINE ptdslv(r,n,l,w,nd) + INTEGER,INTENT(IN) :: n,nd + REAL,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL,INTENT(IN),DIMENSION(:) :: l,w + INTEGER i,j + REAL, DIMENSION(nd) :: sum + DO j=1,nd + sum(j)=w(1)*r(j,1) + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n-1 + DO j=1,nd + r(j,i)=r(j,i)-r(j,i-1) + sum(j)=sum(j)+w(i)*r(j,i) + r(j,i)=r(j,i)*l(i) + ENDDO + ENDDO + DO j=1,nd + r(j,n)=l(n)*(r(j,n)-sum(j)) + r(j,n-1)=r(j,n-1)-w(n-1)*r(j,n) + ENDDO + DO i=n-2,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1)-w(i)*r(j,n) + ENDDO + ENDDO + END SUBROUTINE ptdslv + + SUBROUTINE d2fp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL,INTENT(IN) :: h + REAL,INTENT(IN),DIMENSION(nd,n) :: x + REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL :: h2,r1,r2,t1,t2 + h2=1./(h*h) + r1=6. + r2=3./8. + DO j=1,nd + t1 = (x(j,n)-2.*x(j,n-1)+x(j,n-2)) + t2 = (x(j,1)-2.*x(j,n-1)+x(j,n-3)) + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) t1=0. + IF (x(j,1).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-3)) t2=0. + dx(j,n-1)=(r1*t1+r2*t2) + + t1 = (x(j,1)-2.*x(j,n)+x(j,n-1)) + t2 = (x(j,2)-2.*x(j,n)+x(j,n-2)) + IF (x(j,1).eq.x(j,n).and.x(j,n).eq.x(j,n-1)) t1=0. + IF (x(j,2).eq.x(j,n).and.x(j,n).eq.x(j,n-2)) t2=0. +! dx(j,n)=(r1*(x(j,1)-2.*x(j,n)+x(j,n-1)) & +! +r2*(x(j,2)-2.*x(j,n)+x(j,n-2))) + dx(j,n)=(r1*t1+r2*t2) + + t1 = (x(j,2)-2.*x(j,1)+x(j,n)) + t2 = (x(j,3)-2.*x(j,1)+x(j,n-1)) + IF (x(j,2).eq.x(j,1).and.x(j,1).eq.x(j,n)) t1=0. + IF (x(j,3).eq.x(j,1).and.x(j,1).eq.x(j,n-1)) t2=0. +! dx(j,1)=(r1*(x(j,2)-2.*x(j,1)+x(j,n)) & +! +r2*(x(j,3)-2.*x(j,1)+x(j,n-1))) + dx(j,1)=(r1*t1+r2*t2) + + t1 = (x(j,3)-2.*x(j,2)+x(j,1)) + t2 = (x(j,4)-2.*x(j,2)+x(j,n)) + IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) t1=0. + IF (x(j,4).eq.x(j,2).and.x(j,2).eq.x(j,n)) t2=0. +! dx(j,2)=(r1*(x(j,3)-2.*x(j,2)+x(j,1)) & +! +r2*(x(j,4)-2.*x(j,2)+x(j,n))) + dx(j,2)=(r1*t1+r2*t2) + + dx(j,n-1)=dx(j,n-1)*h2 + dx(j,n)=dx(j,n)*h2 + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*h2 + ENDDO + DO i=3,n-2 + DO j=1,nd + t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1)) + t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2)) + IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0. + IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0. +! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) & +! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2))) + dx(j,i)=(r1*t1+r2*t2) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + IF (dir.eq.1) CALL ptdslv(dx,n,lxs,wxs,nd) ! x-direction + IF (dir.eq.2) CALL ptdslv(dx,n,lys,wys,nd) ! y-direction + IF (dir.eq.3) CALL ptdslv(dx,n,lzs,wzs,nd) ! z-direction + END SUBROUTINE d2fp + + SUBROUTINE tdslv(r,n,l,nd) + INTEGER,INTENT(IN) :: n,nd + REAL,INTENT(INOUT),DIMENSION(nd,n) :: r + REAL,INTENT(IN),DIMENSION(:) :: l + INTEGER i,j + REAL t1 + DO j=1,nd + r(j,1)=r(j,1)*l(1) + ENDDO + DO i=2,n + DO j=1,nd + t1=r(j,i)-r(j,i-1) + r(j,i)=l(i)*t1 + ENDDO + ENDDO + DO i=n-1,1,-1 + DO j=1,nd + r(j,i)=r(j,i)-l(i)*r(j,i+1) + ENDDO + ENDDO + END SUBROUTINE tdslv + + SUBROUTINE d2fnonp(n,h,x,dx,nd,dir) + INTEGER,INTENT(IN) :: n,nd,dir + REAL,INTENT(IN) :: h + REAL,INTENT(IN),DIMENSION(nd,n) :: x + REAL,INTENT(OUT),DIMENSION(nd,n) :: dx + INTEGER :: i,j + REAL :: h2,r1,r2,r3,a,b,c,e,t1,t2 + + h2=1./(h*h) + r1=6. + r2=3./8. + r3=12. + a=13./11. + b=-27./11. + c=15./11. + e=-1./11. + + DO j=1,nd + dx(j,1)=(a*x(j,1)+b*x(j,2)+c*x(j,3)+e*x(j,4)) + dx(j,2)=(x(j,3)-2.*x(j,2)+x(j,1)) + dx(j,n-1)=(x(j,n)-2.*x(j,n-1)+x(j,n-2)) + dx(j,n)=(a*x(j,n)+b*x(j,n-1)+c*x(j,n-2)+e*x(j,n-3)) + IF (x(j,1).eq.x(j,2).and.x(j,2).eq.x(j,3).and.x(j,3).eq.x(j,4)) dx(j,1)=0. + IF (x(j,3).eq.x(j,2).and.x(j,2).eq.x(j,1)) dx(j,2)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2).and.x(j,n-2).eq.x(j,n-3)) dx(j,n)=0. + IF (x(j,n).eq.x(j,n-1).and.x(j,n-1).eq.x(j,n-2)) dx(j,n-1)=0. + dx(j,1)=dx(j,1)*h2 + dx(j,2)=dx(j,2)*h2*r3 + dx(j,n-1)=dx(j,n-1)*h2*r3 + dx(j,n)=dx(j,n)*h2 + ENDDO + DO i=3,n-2 + DO j=1,nd + t1 = (x(j,i+1)-2.*x(j,i)+x(j,i-1)) + t2 = (x(j,i+2)-2.*x(j,i)+x(j,i-2)) + IF (x(j,i+1).eq.x(j,i).and.x(j,i).eq.x(j,i-1)) t1=0. + IF (x(j,i+2).eq.x(j,i).and.x(j,i).eq.x(j,i-2)) t2=0. + +! dx(j,i)=(r1*(x(j,i+1)-2.*x(j,i)+x(j,i-1)) & +! +r2*(x(j,i+2)-2.*x(j,i)+x(j,i-2))) + dx(j,i)=(r1*t1+r2*t2) + dx(j,i)=dx(j,i)*h2 + ENDDO + ENDDO + + IF (dir.eq.1) CALL tdslv(dx,n,lxs,nd) ! x-direction + IF (dir.eq.2) CALL tdslv(dx,n,lys,nd) ! y-direction + IF (dir.eq.3) CALL tdslv(dx,n,lzs,nd) ! z-direction + + END SUBROUTINE d2fnonp + + END MODULE Compact diff --git a/code/makefile b/code/makefile new file mode 100644 index 0000000..e183b87 --- /dev/null +++ b/code/makefile @@ -0,0 +1,20 @@ +#flags = -r8 -i8 -msse3 -O3 -axT -m64 -fPIC -i_dynamic -openmp +#Neuman flags +#flags = -O3 -r8 -i8 -openmp -msse3 -axSSSE3 -m64 -fPIC -i_dynamic +#Cluster_2(16.161) flags +flags = -O3 -r8 -i8 -qopenmp -msse3 -axSSSE3 -m64 -fPIC +compiler = ifort + +Post_Cbar_FSD_Incomp : test.o post.mod Compact.mod + ${compiler} -o x-edge-cold-bc-uPrime-hybrid ${flags} test.o post.o Compact.o + +test.o : test.f90 post.mod + ${compiler} -c ${flags} test.f90 + +post.mod: post.f90 Compact.mod + ${compiler} -c ${flags} post.f90 + +Compact.mod: Compact.f90 + ${compiler} -c ${flags} Compact.f90 +clean: + rm *.o *.mod x-edge-cold-bc-uPrime-hybrid diff --git a/code/post.f90 b/code/post.f90 new file mode 100644 index 0000000..cae10c8 --- /dev/null +++ b/code/post.f90 @@ -0,0 +1,1778 @@ + MODULE post + USE Compact + IMPLICIT NONE + PRIVATE + INTEGER :: nxp,nyp,nzp,startnum,endnum,skipnum,ncyc + INTEGER :: syp,eyp,startsw,twod + INTEGER :: omitnum,countnum + INTEGER :: shiftnum ! hybrid + INTEGER :: num_, dummyu_ ! hybrid + REAL :: tnow,rod,SL_u + REAL :: hxp,hyp,hzp,l_0 + REAL :: scp,prp,lep,vis0p + REAL :: refwr,prof_wr,min_wr,min_fsd + REAL :: pre,ac,bc,c_cut,c_ref + REAL :: min_c + 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 :: 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 + + INTEGER, DIMENSION(:,:), ALLOCATABLE :: omit_t + + PUBLIC :: main + CONTAINS + + SUBROUTINE main + INTEGER :: fread,i,postsw + + CALL READ_INTRO + CALL PREPARE_SAVE + + ! OPEN (210,FILE="ui.dat") + ! WRITE(210,*) 'VARIABLES = "x","y","z","ui"' + countnum=0 + + DO fread=startnum,endnum,skipnum + + postsw=1 + IF (omitnum.gt.0) THEN + DO i=1,omitnum + IF (fread.ge.omit_t(i,1).and.fread.le.omit_t(i,2)) THEN + postsw=-1*i + ENDIF + ENDDO + ENDIF + + IF (postsw.lt.0) 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((-1*postsw),1),' <= fullsavenum <= ',omit_t((-1*postsw),2) + ELSEIF (postsw.eq.1) THEN + 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 + + write(*,*) 'After 1st loop' + CALL AVERAGING + + CALL SECOND_LOOP + + ! 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 READ_INTRO + CHARACTER(LEN=10) :: cdum + INTEGER :: ierr,i + + OPEN (100,FILE='post-edge-cold-bc-hybrid-intro') + OPEN (101,FILE='otape') + READ (100,*) cdum,l_0 + WRITE(101,*) cdum,l_0 + READ (100,*) cdum,startnum + WRITE(101,*) cdum,startnum + READ (100,*) cdum,endnum + WRITE(101,*) cdum,endnum + READ (100,*) cdum,skipnum + WRITE(101,*) cdum,skipnum + READ (100,*) cdum,shiftnum + WRITE(101,*) cdum,shiftnum + READ (100,*) cdum,vis0p + WRITE(101,*) cdum,vis0p + READ (100,*) cdum,scp + WRITE(101,*) cdum,scp + READ (100,*) cdum,lep + WRITE(101,*) cdum,lep + READ (100,*) cdum,min_wr + WRITE(101,*) cdum,min_wr + READ (100,*) cdum,prof_wr + WRITE(101,*) cdum,prof_wr + READ (100,*) cdum,min_fsd + WRITE(101,*) cdum,min_fsd + READ (100,*) cdum,min_c + WRITE(101,*) cdum,min_c + READ (100,*) cdum,pre + WRITE(101,*) cdum,pre + READ (100,*) cdum,ac + WRITE(101,*) cdum,ac + READ (100,*) cdum,bc + WRITE(101,*) cdum,bc + READ (100,*) cdum,c_cut + WRITE(101,*) cdum,c_cut + READ (100,*) cdum,c_ref + WRITE(101,*) cdum,c_ref + READ (100,*) cdum,syp + WRITE(101,*) cdum,syp + READ (100,*) cdum,eyp + WRITE(101,*) cdum,eyp + READ (100,*) cdum,SL_u + WRITE(101,*) cdum,SL_u + READ (100,*) cdum,omitnum + WRITE(101,*) cdum,omitnum + IF (omitnum.gt.0) THEN + ALLOCATE(omit_t(omitnum,2),STAT=ierr) + DO i=1,omitnum + READ (100,*) cdum,omit_t(i,1) ! begining fullsave number + WRITE(101,*) cdum,omit_t(i,1) + READ (100,*) cdum,omit_t(i,2) ! ending fullsave number + WRITE(101,*) cdum,omit_t(i,2) + ENDDO + ENDIF + + prp=scp/lep + CLOSE (100) + CLOSE (101) + + + END SUBROUTINE READ_INTRO + + SUBROUTINE PREPARE_SAVE + 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,'"' + + startsw=0. + + END SUBROUTINE PREPARE_SAVE + + SUBROUTINE READ_FILE(num) + INTEGER, INTENT(IN) :: num + INTEGER :: ierr,i,j,k + REAL, DIMENSION(2) :: tmpr + REAL :: tmp1,tmp2 + REAL :: dt,dummyu + + OPEN(num,FORM='unformatted',STATUS='unknown') + READ (num) tnow,nxp,nyp,nzp,tmp1,tmp2 + + IF (startsw.eq.0) THEN +! syp=1 +! eyp=nyp +! Arrays for local values + 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 + + CALL ludcmp(nxp,nyp,nzp,1,0,0) ! 1,1,0 + + WRITE(*,'(a6,i3,a8,i3,a8,i3)') ' NX = ',nxp,' / NY = ',nyp,' / NZ = ',nzp + WRITE(*,*) ' Preparing memory space for COMPACT SCHEME' + + + WRITE(*,'(a22,i3,a3,i3,a4,i3)') ' Grid number range : ',syp,' ~ ',eyp,' of ',nyp + WRITE(*,*) + + twod=0 + IF (nzp.eq.1) twod=1 + startsw=1 + + ENDIF + + !READ (num) ncyc,tmpr(1:2) + 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) m_v(2,:,:,:),m_v(3,:,:,:),m_v(4,:,:,:),m_v(5:6,:,:,:) + ELSE + WRITE(*,*) ' with a new fort data from Comb-Cluster' + READ (num) m_v(2,:,:,:),m_v(3,:,:,:),m_v(4,:,:,:),m_v_new(:,:,:,1:2) + m_v(2,:,:,:) = m_v(2,:,:,:) + dummyu + ENDIF + + 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 / +! / / / / / / / + +! IRE2 : 1 / 2 / 3 / +! f / f / 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 / +! b / b / b / b / + +! IRE5 : 1 / 2 / 3 / 4 / +! u / u / u / 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) ! k * +! 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) ! f +! +! IRE2(27,i)=IRE2(19,i)*IRE2(16,i)+IRE2(20,i)*IRE2(17,i)+IRE2(21,i)*IRE2(18,i) ! f dot f +! IRE2(28,i)=IRE2(22,i)-IRE2(27,i) ! f +! IRE2(29,i)=IRE2(16,i)**2.+IRE2(17,i)**2.+IRE2(18,i)**2. ! f dot f +! IRE2(30,i)=IRE2(23,i)-IRE2(29,i) ! f +! IRE2(31,i)=SQRT(IRE2(29,i)) ! |f| = abs(f) +! IRE2(32,i)=IRE2(2,i)+IRE2(22,i) ! f +! ELSE +! IRE2(2:32,i)=0. ! To avoid infinity +! ENDIF +! ENDDO + +!! df/dx +! DO i=1,nxp +! ux(1,i)=IRE2(16,i) ! f +! ENDDO +! CALL dfnonp(nxp,hxp,ux,dux,1,1) +! DO i=1,nxp +! IRE2(33,i)=dux(1,i) ! df/dx +! ENDDO + +! Conditional average, < >b, < >u + IRE3(1:4,:)=IRE3(1:4,:)/ndata ! + 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 READ_FILE2(num) + INTEGER, INTENT(IN) :: num + INTEGER, DIMENSION(4) :: itmp + REAL, DIMENSION(3) :: rtmp + REAL, DIMENSION(2) :: tmpr + REAL :: dt,dummyu + + OPEN(num,FORM='unformatted',STATUS='unknown') + READ (num) tnow,nxp,nyp,nzp,tmpr(1:2) + +! Bgn - DQ version +! READ (num) ncyc,tmpr(1:2) +! READ (num) tmpr(1:4) +! READ (num) tmpr(1:4) +! READ (num) tmpr(1:2) +! End - DQ version ------------------------- + + !READ (num) ncyc,tmpr(1:2) + READ (num) ncyc,dt,dummyu + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + READ (num) tmpr(1:2) + + WRITE(*,'(a27,f8.3,a1,i6,a4,i5,a3,i5)') ' 2nd LOOP : T & NCYC = ',tnow,',' & + ,itmp(4),' || ',(num-startnum+1),' / ',(endnum-startnum+1) + + m_v(1,:,:,:)=1. + + num_=num + IF(num.le.shiftnum) THEN + WRITE(*,*) ' Reading an old fort data from Nueman-0X' + READ (num) m_v(2,:,:,:),m_v(3,:,:,:),m_v(4,:,:,:),m_v(5:6,:,:,:) + ELSE + WRITE(*,*) ' Reading a new fort data from Comb-Cluster' + READ (num) m_v(2,:,:,:),m_v(3,:,:,:),m_v(4,:,:,:),m_v_new(:,:,:,1:2) + m_v(2,:,:,:) = m_v(2,:,:,:) + dummyu + ENDIF + CLOSE (num) + + + + + END SUBROUTINE READ_FILE2 + + SUBROUTINE SECOND_LOOP + INTEGER :: i,j,k,iv,post_sw + REAL :: ntime,ui,vi,wi +!------------------------------------------------------------- +! X_DOT_F : fluctuation of X by surface average. +!---------------------------------------------------------------- + + DO iv=startnum,endnum,skipnum + +! CALL READ_FILE2(iv) + + post_sw=1 + IF (omitnum.gt.0) THEN + DO i=1,omitnum + IF (iv.ge.omit_t(i,1).and.iv.le.omit_t(i,2)) THEN + post_sw=i*-1 + ENDIF + ENDDO + ENDIF + + IF (post_sw.lt.0) THEN + WRITE(*,'(a40,i7,a4,i4,a3,i4)') ' Current fullsavenum = ',iv,& + ' || ',(iv-startnum+1),' / ',(endnum-startnum+1) + WRITE(*,'(a12,i6,a20,i6)') & + ' Skip. ',omit_t((-1*post_sw),1),' <= fullsavenum <= ',omit_t((-1*post_sw),2) + + ELSEIF (post_sw.eq.1) THEN + + CALL READ_FILE2(iv) + CALL CAL_Yrs + !CALL CAL_CGM_N !dhkim + CALL CAL_FLUCTUATION(iv) + CALL SAVE_SUM_FLUCTUATION + ENDIF + + ENDDO ! Time-loop + + CALL FLUCTUATION_AVG + + END SUBROUTINE SECOND_LOOP + + SUBROUTINE CAL_FLUCTUATION(num) + INTEGER :: i,j,k,jj,num + + 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 ! [] +!====================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) ! [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] +! 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) ! +! 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/)*(d/dx) = 1/Lw1 +! IF(IRE1(1,i).le.0.) IRE8(1,i)=0. +! IRE8(2,i)=dux(1,i)/IRE1(2,i) ! (1/)*(d/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/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/|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-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/dx)*(d2/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) +! ! dot Grad. +! IRE9(2,i) = IRE8(16,i)*IRE1(23,i) ! Dts*Lap. +! IRE9(3,i) = IRE2(2,i)*IRE1(2,i) ! f* +! IRE9(4,i) = IRE8(16,i)*IRE1(2,i)*IRE2(33,i) ! Dts**Div.(f) +! IRE9(5,i) = IRE2(28,i)*IRE1(2,i) ! f* +! 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-)*(d(1-)/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) !
+! ux(2,i)=IRE8(3,i)+IRE1(10,i) !
+ +! 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) ! K +! if(IRE1(19,i).eq.0.) IRE9(7,i)=0. +! IRE9(8,i) = IRE2(16,i)/IRE9(7,i) ! f/K +! if(IRE9(7,i).eq.0.) IRE9(8,i)=0. +! IRE9(9,i) = dux(1,i) ! d
/dx +! IRE9(10,i) = dux(2,i) ! d/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-)*(d(1-)/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","","","","","","",""' ! 8 + !WRITE(200,*) 'VARIABLES = "X","","","","","","","",""' ! 9 +! WRITE(200,*) '"","","","","","","",""' ! 8 -> 17 +! WRITE(200,*) '"","_g",""' ! 3 -> 20 +! WRITE(200,*) '"d(1-)/dx","d(1-)/dy","d(1-)/dz"' ! 3 -> 23 +! WRITE(200,*) '"","","",""' ! 4 -> 27 +! WRITE(200,*) '"",""' ! 2 -> 29 + +! IRE2 +! WRITE(200,*) '"","f","f","f","f","f","f","f"' ! 8 -> 8 +! WRITE(200,*) '"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","f","f","f","f","f","f"' ! 7 -> 21 +! WRITE(200,*) '"f","f","f","f","f"' ! 5 -> 26 +! WRITE(200,*) '"f dot f","f","f dot f","f"' ! 4 -> 30 +! WRITE(200,*) '"|f|","f","df/dx"' ! 3 -> 33 +! WRITE(200,*) '"-f"' ! 1 -> 34 + +! IRE3 + WRITE(200,*) '"b","b","b","b","b_g","b_g","b_g","b_g"' ! 8 -> 8 +! IRE4 + WRITE(200,*) '"u","u","u","u","u_g","u_g","u_g","u_g"' ! 8 -> 8 +! IRE5 + WRITE(200,*) '"RMS(u`)","RMS(v`)","RMS(w`)",""' ! 4 -> 4 + WRITE(200,*) '"RMS(u`)b","RMS(v`)b","RMS(w`)b","b"' ! 4 -> 8 + WRITE(200,*) '"RMS(u`)u","RMS(v`)u","RMS(w`)u","u"' ! 4 -> 12 + !WRITE(200,*) '"","",""' ! 3 -> 15 + WRITE(200,*) '"","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","_g"' ! 4 -> 4 + WRITE(200,*) '"RMS(u`)b_g","RMS(v`)b_g","RMS(w`)b_g","b_g"' ! 4 -> 8 + WRITE(200,*) '"RMS(u`)u_g","RMS(v`)u_g","RMS(w`)u_g","u_g"' ! 4 -> 12 + WRITE(200,*) '"_g","_g","_g"' ! 3 -> 15 +!! IRE7 +! WRITE(200,*) '"k","k","k"' ! 3 -> 3 +! WRITE(200,*) '"k","k","k"' ! 3 -> 6 +!! IRE8 +! WRITE(200,*) '"(1/)*(d/dx)","(1/)*(d/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/dx","Dts"' ! 2 -> 16 +! WRITE(200,*) '"(1/(1-))*(d(1-)/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-f","ST_9"' ! 2 -> 24 +! WRITE(200,*) '"1/L_LE_6=1/(d/dx)*(d2/dx2)","ST_10"' ! 2 -> 26 +! +!! IRE9 +! WRITE(200,*) '"dotGrad.","Dts*Lap.","f*"' ! 3 -> 3 +! WRITE(200,*) '"Dts**Div.(f)","f*","cEqnBalance"' ! 3 -> 6 +! WRITE(200,*) '"K","f/K","d
/dx"' ! 3 -> 9 +! WRITE(200,*) '"d/dx","1/(d<1-c>/dx)*(d2<1-c>/dx2)"' ! 2 -> 11 +! WRITE(200,*) '"1/(1-)*(d(1-)/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 DEALLOCATES_CLOSE + 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 diff --git a/code/test.f90 b/code/test.f90 new file mode 100644 index 0000000..bbe2343 --- /dev/null +++ b/code/test.f90 @@ -0,0 +1,11 @@ +! This code has been written since 2011.01.08 +! This code is written by Dongkyu Lee (THDHRKDG) +! Purpose of this code is to postprocess results of stagnating DNS code + + PROGRAM test + USE post + USE Compact + + CALL main + + END PROGRAM