From b426be89adfbaff0ba5ae619b25cfc671abb8f0b Mon Sep 17 00:00:00 2001 From: ignis Date: Fri, 10 Nov 2017 16:37:35 +0900 Subject: [PATCH] transpose subroutine convention change and minor mods --- m_fdm_calc.f90 | 48 +++++++++++++++++++++++++++--------------------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/m_fdm_calc.f90 b/m_fdm_calc.f90 index 6d63cfd..24b236d 100644 --- a/m_fdm_calc.f90 +++ b/m_fdm_calc.f90 @@ -489,6 +489,9 @@ module m_fdm_calc real*8 :: uuz(zz) real*8 :: y + + f_(:,:,:,1) = 0.0 ! continuity + ! reaction source term DO k=1,zz DO j=1,yy @@ -503,7 +506,6 @@ module m_fdm_calc min_wr-refwr*exp(prof_wr*(c_cut-c_ref)))/(1.-exp(prof_wr*(c_cut-c_ref))) ENDIF - f_(i,j,k,1) = 0.0 ! continuity f_(i,j,k,2) = - wrate ! species conservation ENDDO @@ -512,6 +514,8 @@ module m_fdm_calc !! z-direction + +! sequential run if (numprocs.eq.1) then DO j=1,yy @@ -539,6 +543,7 @@ module m_fdm_calc ENDDO ENDDO +! parallel run else ! -( d(rho*w)/dz ) @@ -624,29 +629,29 @@ module m_fdm_calc DO k=1,zz ! -( d(rho*u)/dx ) - CALL tp2mul (r1_(:,:,k,1), uu_(:,:,k), yxbuf1, xx, yy) ! rho*u + CALL tp2mul (yxbuf1, r1_(:,:,k,1), uu_(:,:,k), yy, xx) ! rho*u CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u) - CALL tp2sub (f_(:,:,k,1), yxbuf2, xx, yy) ! continuity + CALL tp2subasgn (f_(:,:,k,1), yxbuf2, xx, yy) ! continuity ! -( d(rho*u*Yr)/dx ) + d(rho*D* d(Yr)/dx)/dx ! = -( d(rho*u*Yr)/dx ) ! + D* (rho* d2(Yr)/dx2 + d(rho)/dx * d(Yr)/dx ) - CALL tp2mul (r1_(:,:,k,2), uu_(:,:,k), yxbuf1, xx, yy) ! rho*u*Y + CALL tp2mul (yxbuf1, r1_(:,:,k,2), uu_(:,:,k), yy, xx) ! rho*u*Y CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(rho*u*Y) - CALL tp2sub (f_(:,:,k,2), yxbuf2, xx, yy) ! species conservation + CALL tp2subasgn (f_(:,:,k,2), yxbuf2, xx, yy) ! species conservation - CALL tp2div (r1_(:,:,k,2), r1_(:,:,k,1), yxbuf1, xx, yy)! Y + CALL tp2div (yxbuf1, r1_(:,:,k,2), r1_(:,:,k,1), yy, xx)! Y CALL dfnonp(xx,hx,yxbuf1,yxbuf2,yy,1) ! d/dx(Y) CALL d2fnonp(xx,hx,yxbuf1,yxbuf3,yy,1) ! d2/dx2(Y) - CALL tp2 (r1_(:,:,k,1), yxbuf1, xx, yy) ! rho + CALL tp2 (yxbuf1, r1_(:,:,k,1), yy, xx) ! rho CALL dfnonp(xx,hx,yxbuf1,yxbuf4,yy,1) ! d/dx(rho) yxbuf2 = yxbuf2 * yxbuf4 ! d/dx(Y) * d/dx(rho) yxbuf1 = -diff*(yxbuf2 + yxbuf1 * yxbuf3) ! -D( ... + rho * d2/dx2(Y)) - CALL tp2sub (f_(:,:,k,2), yxbuf1, xx, yy) ! species conservation + CALL tp2subasgn (f_(:,:,k,2), yxbuf1, xx, yy) ! species conservation ENDDO @@ -794,12 +799,13 @@ module m_fdm_calc END SUBROUTINE substep subroutine tp2 (a, b, n1, n2) + ! a = transpose(b) implicit none integer,intent(in) :: n1, n2 - real*8,intent(in) :: a(n1,n2) - real*8,intent(out) :: b(n2,n1) + real*8,intent(out) :: a(n1,n2) + real*8,intent(in) :: b(n2,n1) integer :: i,j,ii,jj integer,parameter :: nb = 16 @@ -808,7 +814,7 @@ module m_fdm_calc DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 - b(j,i) = a(i,j) + a(i,j) = b(j,i) ENDDO ENDDO @@ -819,13 +825,13 @@ module m_fdm_calc subroutine tp2div (a, b, c, n1, n2) - ! c = a / b + ! a = transpose(b/c) implicit none integer,intent(in) :: n1, n2 - real*8,intent(in) :: a(n1,n2), b(n1,n2) - real*8,intent(out) :: c(n2,n1) + real*8,intent(in) :: b(n2,n1), c(n2,n1) + real*8,intent(out) :: a(n1,n2) integer :: i,j,ii,jj integer,parameter :: nb = 16 @@ -834,7 +840,7 @@ module m_fdm_calc DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 - c(j,i) = a(i,j) / b(i,j) + a(i,j) = b(j,i) / c(j,i) ENDDO ENDDO @@ -846,13 +852,13 @@ module m_fdm_calc subroutine tp2mul (a, b, c, n1, n2) - ! c = a * b + ! a = transpose(b*c) implicit none integer,intent(in) :: n1, n2 - real*8,intent(in) :: a(n1,n2), b(n1,n2) - real*8,intent(out) :: c(n2,n1) + real*8,intent(in) :: b(n2,n1), c(n2,n1) + real*8,intent(out) :: a(n1,n2) integer :: i,j,ii,jj integer,parameter :: nb = 16 @@ -861,7 +867,7 @@ module m_fdm_calc DO j=jj,jj+nb-1 DO i=ii,ii+nb-1 - c(j,i) = a(i,j) * b(i,j) + a(i,j) = b(j,i) * c(j,i) ENDDO ENDDO @@ -871,7 +877,7 @@ module m_fdm_calc end subroutine tp2mul - subroutine tp2sub (a, b, n1, n2) + subroutine tp2subasgn (a, b, n1, n2) ! a = a - transpose(b) implicit none @@ -894,7 +900,7 @@ module m_fdm_calc ENDDO ENDDO - end subroutine tp2sub + end subroutine tp2subasgn end module m_fdm_calc