diff --git a/src/BC-Channel-flow.f90 b/src/BC-Channel-flow.f90 index a9fdfc3a2..747c132c6 100644 --- a/src/BC-Channel-flow.f90 +++ b/src/BC-Channel-flow.f90 @@ -265,12 +265,15 @@ subroutine boundary_conditions_channel (ux,uy,uz,phi) use var, only : di2 use variables use decomp_2d + use mhd, only : Bm, mhd_active, mhd_equation implicit none real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: ux,uy,uz real(mytype),dimension(xsize(1),xsize(2),xsize(3),numscalar) :: phi + integer j + if (.not. cpg ) then ! if not constant pressure gradient if (idir_stream == 1) then call channel_cfr(ux,two/three) @@ -300,6 +303,15 @@ subroutine boundary_conditions_channel (ux,uy,uz,phi) !if (nclySn.eq.2) g_sc(:,2) = zero endif endif + + if( mhd_active .and. iimplicit==0 .and. mhd_equation ) then + Bm(:,1,:,2) = one + Bm(:,ny,:,2) = one + Bm(:,1,:,1) = zero + Bm(:,ny,:,1) = zero + Bm(:,1,:,3) = zero + Bm(:,ny,:,3) = zero + endif end subroutine boundary_conditions_channel !############################################################################ @@ -394,7 +406,7 @@ subroutine visu_channel(ux1, uy1, uz1, pp3, phi1, ep1, num) use var, only : ta2,tb2,tc2,td2,te2,tf2,di2,ta3,tb3,tc3,td3,te3,tf3,di3 use var, ONLY : nzmsize use visu, only : write_field - use mhd, only : mhd_active,Je + use mhd, only : mhd_active,Je, Bm use ibm_param, only : ubcx,ubcy,ubcz @@ -458,6 +470,10 @@ subroutine visu_channel(ux1, uy1, uz1, pp3, phi1, ep1, num) call write_field(Je(:,:,:,1), ".", "J_x", num, flush = .true.) call write_field(Je(:,:,:,2), ".", "J_y", num, flush = .true.) call write_field(Je(:,:,:,3), ".", "J_z", num, flush = .true.) + call write_field(Bm(:,:,:,1), ".", "B_x", num, flush = .true.) + call write_field(Bm(:,:,:,2), ".", "B_y", num, flush = .true.) + call write_field(Bm(:,:,:,3), ".", "B_z", num, flush = .true.) + endif end subroutine visu_channel diff --git a/src/implicit.f90 b/src/implicit.f90 index ec5320f3d..39f7aaad4 100644 --- a/src/implicit.f90 +++ b/src/implicit.f90 @@ -23,7 +23,7 @@ module ludecomp module procedure ludecomp7_12 module procedure ludecomp7_0 end interface ludecomp7 - + interface ludecomp9 module procedure ludecomp9_12 module procedure ludecomp9_0 @@ -31,15 +31,15 @@ module ludecomp contains -!******************************************************************* -!Decomposition LU matrix septadiag + !******************************************************************* + !Decomposition LU matrix septadiag subroutine ludecomp7_12(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ny) -! -!******************************************************************* + ! + !******************************************************************* use decomp_2d, only : mytype USE param implicit none - + integer :: i,j,k integer, intent(in) :: ny real(mytype),dimension(ny), intent(in) :: aam,bbm,ccm,ddm,eem,rrm,qqm @@ -78,14 +78,14 @@ subroutine ludecomp7_12(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ny) end subroutine ludecomp7_12 -!******************************************************************* -!Decomposition LU matrix cyclic septadiag + !******************************************************************* + !Decomposition LU matrix cyclic septadiag subroutine ludecomp7_0(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3m,ny) -! -!******************************************************************* + ! + !******************************************************************* use decomp_2d, only : mytype USE param - + implicit none integer :: i,j,k @@ -202,11 +202,11 @@ subroutine ludecomp7_0(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l end subroutine ludecomp7_0 -!******************************************************************* -!Decomposition LU matrix nonadiag + !******************************************************************* + !Decomposition LU matrix nonadiag subroutine ludecomp9_12(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ttm,uum,sssm,zzzm,ny) -! -!******************************************************************* + ! + !******************************************************************* use decomp_2d, only : mytype use param @@ -263,11 +263,11 @@ subroutine ludecomp9_12(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ttm, enddo end subroutine ludecomp9_12 -!******************************************************************* -!Decomposition LU matrix cyclic nonadiag + !******************************************************************* + !Decomposition LU matrix cyclic nonadiag subroutine ludecomp9_0(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3m,ny) -! -!******************************************************************* + ! + !******************************************************************* use decomp_2d, only : mytype use param USE MPI @@ -327,8 +327,8 @@ module matinv contains -!******************************************************************** -!INVERSE SEPTADIAG MATRIX + !******************************************************************** + !INVERSE SEPTADIAG MATRIX subroutine septinv_12(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,nx,ny,nz) ! !******************************************************************** @@ -342,7 +342,7 @@ subroutine septinv_12(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,nx,ny,nz) real(mytype),dimension(nx,ny,nz),intent(in) :: bbb real(mytype),dimension(ny),intent(in) :: ggm,hhm,ssm,rrm real(mytype),dimension(ny),intent(in) :: vvm,wwm,zzm - + do k=1,nz !going down do i=1,nx @@ -385,11 +385,11 @@ subroutine septinv_12(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,nx,ny,nz) end subroutine septinv_12 -!******************************************************************** -!INVERSE SEPTADIAG CYCLIC + !******************************************************************** + !INVERSE SEPTADIAG CYCLIC subroutine septinv_0(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3m,nx,ny,nz) -! -!******************************************************************** + ! + !******************************************************************** use decomp_2d, only : mytype implicit none @@ -451,12 +451,12 @@ subroutine septinv_0(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3 enddo end subroutine septinv_0 - -!******************************************************************** -!INVERSE NONADIAG + + !******************************************************************** + !INVERSE NONADIAG subroutine nonainv_12(xSol,bbb,ggm,hhm,ssm,sssm,ttm,zzzm,zzm,wwm,vvm,nx,ny,nz) -! -!******************************************************************** + ! + !******************************************************************** use decomp_2d, only : mytype implicit none @@ -500,14 +500,14 @@ subroutine nonainv_12(xSol,bbb,ggm,hhm,ssm,sssm,ttm,zzzm,zzm,wwm,vvm,nx,ny,nz) end subroutine nonainv_12 -!******************************************************************** -!INVERSE NONADIAG CYCLIC + !******************************************************************** + !INVERSE NONADIAG CYCLIC subroutine nonainv_0(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3m,nx,ny,nz) ! !******************************************************************** use decomp_2d, only : mytype USE MPI - + implicit none integer :: i,j,k,kk @@ -524,7 +524,7 @@ subroutine nonainv_0(xsol,bbb,ggm,hhm,ssm,rrm,vvm,wwm,zzm,l1m,l2m,l3m,u1m,u2m,u3 call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop end subroutine nonainv_0 - + end module matinv module ydiff_implicit @@ -534,1710 +534,1854 @@ module ydiff_implicit private public :: inttimp, init_implicit, implicit_schemes - contains -! -! Time integration, (semi)implicit Y diffusion -! var1, input : variable at time n -! output : variable at time n+1 -! dvar1 : r.h.s. of transport equation -! npaire : odd / even variable, important when ncly*=1 -! isc : 0 for momentum, id of the scalar otherwise -! forcing1 : r.h.s. term not present in dvar1 (pressure gradient) -! -subroutine inttimp (var1,dvar1,npaire,isc,forcing1) +contains + ! + ! Time integration, (semi)implicit Y diffusion + ! var1, input : variable at time n + ! output : variable at time n+1 + ! dvar1 : r.h.s. of transport equation + ! npaire : odd / even variable, important when ncly*=1 + ! isc : 0 for momentum, id of the scalar otherwise + ! forcing1 : r.h.s. term not present in dvar1 (pressure gradient) + ! + subroutine inttimp (var1,dvar1,npaire,isc,forcing1,mhdvar) - USE MPI - USE param - USE variables - USE var, ONLY: ta1, ta2, tb2, tc2, td2 - USE decomp_2d - use derivY - use matinv + USE MPI + USE param + USE variables + USE var, ONLY: ta1, ta2, tb2, tc2, td2 + USE decomp_2d + use derivY + use matinv - implicit none + implicit none - integer :: i,j,k,code,ierror + integer :: i,j,k,code,ierror - !! IN - real(mytype),dimension(xsize(1),xsize(2),xsize(3)), optional, intent(in) :: forcing1 - integer, intent(in) :: npaire, isc + !! IN + real(mytype),dimension(xsize(1),xsize(2),xsize(3)), optional, intent(in) :: forcing1 + integer, intent(in) :: npaire, isc + integer, optional, intent(in) :: mhdvar !if integrating an MHD variable: 1=bx, 2=by, 3=bz - !! IN/OUT - real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: var1 - real(mytype),dimension(xsize(1),xsize(2),xsize(3),ntime) :: dvar1 + !! IN/OUT + real(mytype),dimension(xsize(1),xsize(2),xsize(3)) :: var1 + real(mytype),dimension(xsize(1),xsize(2),xsize(3),ntime) :: dvar1 - !! LOCAL - real(mytype),dimension(ysize(1),ysize(3)) :: bctop, bcbot + !! LOCAL + real(mytype),dimension(ysize(1),ysize(3)) :: bctop, bcbot - if (itimescheme.eq.1) then - !>>> Explicit Euler - ta1(:,:,:) = gdt(1) * dvar1(:,:,:,1) + if (itimescheme.eq.1) then + !>>> Explicit Euler + ta1(:,:,:) = gdt(1) * dvar1(:,:,:,1) - elseif (itimescheme.eq.2) then - !>>> AB2 - if ((itime.eq.1).and.(irestart.eq.0)) then - !>>> Start with explicit Euler - ta1(:,:,:) = dt*dvar1(:,:,:,1) + elseif (itimescheme.eq.2) then + !>>> AB2 + if ((itime.eq.1).and.(irestart.eq.0)) then + !>>> Start with explicit Euler + ta1(:,:,:) = dt*dvar1(:,:,:,1) - else - !>>> Then AB2 - ta1(:,:,:) = adt(1)*dvar1(:,:,:,1) + bdt(1)*dvar1(:,:,:,2) + else + !>>> Then AB2 + ta1(:,:,:) = adt(1)*dvar1(:,:,:,1) + bdt(1)*dvar1(:,:,:,2) - endif + endif - !save (N+Lxz)(n-1) - dvar1(:,:,:,2) = dvar1(:,:,:,1) + !save (N+Lxz)(n-1) + dvar1(:,:,:,2) = dvar1(:,:,:,1) - else if (itimescheme.eq.3) then - !>>> AB3 - if ((itime.eq.1).and.(irestart.eq.0)) then - !>>> Start with explicit Euler - ta1(:,:,:) = dt*dvar1(:,:,:,1) + else if (itimescheme.eq.3) then + !>>> AB3 + if ((itime.eq.1).and.(irestart.eq.0)) then + !>>> Start with explicit Euler + ta1(:,:,:) = dt*dvar1(:,:,:,1) - else if ((itime.eq.2).and.(irestart.eq.0)) then - !>>> Then AB2 - ta1(:,:,:) = onepfive*dt*dvar1(:,:,:,1) - half*dt*dvar1(:,:,:,2) + else if ((itime.eq.2).and.(irestart.eq.0)) then + !>>> Then AB2 + ta1(:,:,:) = onepfive*dt*dvar1(:,:,:,1) - half*dt*dvar1(:,:,:,2) - else - !>>> Then AB3 - ta1(:,:,:) = adt(1)*dvar1(:,:,:,1) + bdt(1)*dvar1(:,:,:,2) & - + cdt(1)*dvar1(:,:,:,3) + else + !>>> Then AB3 + ta1(:,:,:) = adt(1)*dvar1(:,:,:,1) + bdt(1)*dvar1(:,:,:,2) & + + cdt(1)*dvar1(:,:,:,3) - endif + endif - !save (N+Lxz)(n-1) - dvar1(:,:,:,3)=dvar1(:,:,:,2) - dvar1(:,:,:,2)=dvar1(:,:,:,1) + !save (N+Lxz)(n-1) + dvar1(:,:,:,3)=dvar1(:,:,:,2) + dvar1(:,:,:,2)=dvar1(:,:,:,1) - elseif (itimescheme.eq.4) then - !>>> AB4 - if (nrank.eq.0) then - print *, "AB4 not implemented!" - STOP - endif + elseif (itimescheme.eq.4) then + !>>> AB4 + if (nrank.eq.0) then + print *, "AB4 not implemented!" + STOP + endif - else - !>>> We should not be here - if (nrank.eq.0) then - print *, "Unrecognised implicit itimescheme: ", itimescheme - endif - call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop + else + !>>> We should not be here + if (nrank.eq.0) then + print *, "Unrecognised implicit itimescheme: ", itimescheme + endif + call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop - endif + endif - if (present(forcing1)) then - if ( (irestart.eq.1).or.(itime.gt.1) ) then - ta1(:,:,:) = ta1(:,:,:) - forcing1(:,:,:) - endif - endif + if (present(forcing1)) then + if ( (irestart.eq.1).or.(itime.gt.1) ) then + ta1(:,:,:) = ta1(:,:,:) - forcing1(:,:,:) + endif + endif - !Y-PENCIL FOR MATRIX INVERSION - call transpose_x_to_y(var1,tb2) + !Y-PENCIL FOR MATRIX INVERSION + call transpose_x_to_y(var1,tb2) - call transpose_x_to_y(ta1,ta2) + call transpose_x_to_y(ta1,ta2) - ! - ! Prepare boundary conditions - ! Used only when ncly*=2 - ! velocity - ! isc = 0 - ! Dirichlet BC - ! scalars - ! 1 <= isc <= numscalar - ! Flexible BC, alpha_sc T +- beta_sc dT/dy = g_sc - ! - ! Specific cases first - ! This is the location for exotic / nonhomogeneous boundary conditions - ! - if (itype.eq.itype_tbl .and. isc.eq.0) then - bcbot(:,:) = zero - bctop(:,:) = tb2(:,ny-1,:) - !in order to mimick a Neumann BC at the top of the domain for the TBL - ! - ! Generic homogeneous cases after - ! - else if (isc.ne.0) then - bcbot(:,:) = g_sc(isc, 1) - bctop(:,:) = g_sc(isc, 2) - else - bcbot(:,:) = zero - bctop(:,:) = zero - endif - - !ta2: A.uhat - !td2:(A+xcstB).un - !if isecondder=5, we need nona inversion - !id isecondder is not 5, we need septa inversion - - if (isecondder.ne.5) then - if (isc.eq.0) then - call multmatrix7(td2,ta2,tb2,npaire,ncly1,nclyn,xcst) - else - call multmatrix7(td2,ta2,tb2,npaire,nclyS1,nclySn,xcst_sc(isc)) - endif - else if (isecondder.eq.5) then - !TO BE DONE: Different types of BC - if ((ncly1.eq.0).and.(nclyn.eq.0)) then - !NOT READY - elseif ((ncly1.eq.1).and.(nclyn.eq.1)) then - !NOT READY - elseif ((ncly1.eq.2).and.(nclyn.eq.2)) then - call multmatrix9(td2,ta2,tb2,npaire) - endif - endif - - !Full second member BBB=A uhat+(A+xcst.B)u^n - ta2(:,:,:)=td2(:,:,:)+ta2(:,:,:) + ! + ! Prepare boundary conditions + ! Used only when ncly*=2 + ! velocity + ! isc = 0 + ! Dirichlet BC + ! scalars + ! 1 <= isc <= numscalar + ! Flexible BC, alpha_sc T +- beta_sc dT/dy = g_sc + ! + ! Specific cases first + ! This is the location for exotic / nonhomogeneous boundary conditions + ! + if (itype.eq.itype_tbl .and. isc.eq.0) then + bcbot(:,:) = zero + bctop(:,:) = tb2(:,ny-1,:) + !in order to mimick a Neumann BC at the top of the domain for the TBL + ! + ! Generic homogeneous cases after + ! + else if (isc.gt.0) then + bcbot(:,:) = g_sc(isc, 1) + bctop(:,:) = g_sc(isc, 2) + elseif(present(mhdvar)) then + if(mhdvar.eq.1) then + bcbot(:,:) = zero + bctop(:,:) = zero + elseif(mhdvar.eq.2) then + bcbot(:,:) = one + bctop(:,:) = one + elseif(mhdvar.eq.3) then + bcbot(:,:) = zero + bctop(:,:) = zero + endif + else + bcbot(:,:) = zero + bctop(:,:) = zero + endif + + !ta2: A.uhat + !td2:(A+xcstB).un + !if isecondder=5, we need nona inversion + !id isecondder is not 5, we need septa inversion + + if (isecondder.ne.5) then + if (isc.eq.0) then + call multmatrix7(td2,ta2,tb2,npaire,ncly1,nclyn,xcst) + elseif(present(mhdvar)) then + call multmatrix7(td2,ta2,tb2,npaire,nclyB1(mhdvar),nclyBn(mhdvar),xcstB) + else + call multmatrix7(td2,ta2,tb2,npaire,nclyS1,nclySn,xcst_sc(isc)) + endif + else if (isecondder.eq.5) then + !TO BE DONE: Different types of BC + if ((ncly1.eq.0).and.(nclyn.eq.0)) then + !NOT READY + elseif ((ncly1.eq.1).and.(nclyn.eq.1)) then + !NOT READY + elseif ((ncly1.eq.2).and.(nclyn.eq.2)) then + call multmatrix9(td2,ta2,tb2,npaire) + endif + endif + + !Full second member BBB=A uhat+(A+xcst.B)u^n + ta2(:,:,:)=td2(:,:,:)+ta2(:,:,:) + + ! + ! Apply boundary conditions + ! + if ((isc.eq.0.and.ncly1.eq.2).or.(isc.gt.0.and.nclyS1.eq.2)) then + ta2(:,1,:) = bcbot(:,:) + endif + if ((isc.eq.0.and.nclyn.eq.2).or.(isc.gt.0.and.nclySn.eq.2)) then + ta2(:,ny,:) = bctop(:,:) + endif + + if(present(mhdvar)) then + if ( nclyB1(mhdvar).eq.2 ) then + ta2(:,1,:) = bcbot(:,:) + endif + if ( nclyBn(mhdvar).eq.2 ) then + ta2(:,ny,:) = bctop(:,:) + endif + endif + + !Inversion of the linear system Mx=b: (A-xcst.B)u^n+1=uhat+(A+xcst.B)u^n + !if secondder=5, we need nona inversion + !if isecondder is not 5, we need septa inversion + + if (isecondder.ne.5) then + if ((isc.eq.0).and.(ncly1.eq.0).and.(nclyn.eq.0)) then + gg=>ggm0; hh=>hhm0; ss=>ssm0; rr=>rrm0; vv=>vvm0; ww=>wwm0; zz=>zzm0 + lo1=>l1m; lo2=>l2m; lo3=>l3m; up1=>u1m; up2=>u2m; up3=>u3m + elseif ((isc.gt.0).and.(nclyS1.eq.0).and.(nclySn.eq.0)) then + gg=>ggm0t(:,isc); hh=>hhm0t(:,isc); ss=>ssm0t(:,isc); rr=>rrm0t(:,isc); vv=>vvm0t(:,isc); ww=>wwm0t(:,isc); zz=>zzm0t(:,isc) + lo1=>l1mt(:,isc); lo2=>l2mt(:,isc); lo3=>l3mt(:,isc); up1=>u1mt(:,isc); up2=>u2mt(:,isc); up3=>u3mt(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.1).and.(npaire.eq.0)) then + gg=>ggm10; hh=>hhm10; ss=>ssm10; rr=>rrm10; vv=>vvm10; ww=>wwm10; zz=>zzm10 + elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.1).and.(npaire.eq.0)) then + gg=>ggm10t(:,isc); hh=>hhm10t(:,isc); ss=>ssm10t(:,isc); rr=>rrm10t(:,isc); vv=>vvm10t(:,isc); ww=>wwm10t(:,isc); zz=>zzm10t(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.1).and.(npaire.eq.1)) then + gg=>ggm11; hh=>hhm11; ss=>ssm11; rr=>rrm11; vv=>vvm11; ww=>wwm11; zz=>zzm11 + elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.1).and.(npaire.eq.1)) then + gg=>ggm11t(:,isc); hh=>hhm11t(:,isc); ss=>ssm11t(:,isc); rr=>rrm11t(:,isc); vv=>vvm11t(:,isc); ww=>wwm11t(:,isc); zz=>zzm11t(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.2)) then + gg=>ggm; hh=>hhm; ss=>ssm; rr=>rrm; vv=>vvm; ww=>wwm; zz=>zzm + elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.2)) then + gg=>ggmt(:,isc); hh=>hhmt(:,isc); ss=>ssmt(:,isc); rr=>rrmt(:,isc); vv=>vvmt(:,isc); ww=>wwmt(:,isc); zz=>zzmt(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.2).and.(npaire.eq.0)) then + gg=>ggm120; hh=>hhm120; ss=>ssm120; rr=>rrm120; vv=>vvm120; ww=>wwm120; zz=>zzm120 + elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.2).and.(npaire.eq.0)) then + gg=>ggm120t(:,isc); hh=>hhm120t(:,isc); ss=>ssm120t(:,isc); rr=>rrm120t(:,isc); vv=>vvm120t(:,isc); ww=>wwm120t(:,isc); zz=>zzm120t(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.2).and.(npaire.eq.1)) then + gg=>ggm121; hh=>hhm121; ss=>ssm121; rr=>rrm121; vv=>vvm121; ww=>wwm121; zz=>zzm121 + elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.2).and.(npaire.eq.1)) then + gg=>ggm121t(:,isc); hh=>hhm121t(:,isc); ss=>ssm121t(:,isc); rr=>rrm121t(:,isc); vv=>vvm121t(:,isc); ww=>wwm121t(:,isc); zz=>zzm121t(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.1).and.(npaire.eq.0)) then + gg=>ggm210; hh=>hhm210; ss=>ssm210; rr=>rrm210; vv=>vvm210; ww=>wwm210; zz=>zzm210 + elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.1).and.(npaire.eq.0)) then + gg=>ggm210t(:,isc); hh=>hhm210t(:,isc); ss=>ssm210t(:,isc); rr=>rrm210t(:,isc); vv=>vvm210t(:,isc); ww=>wwm210t(:,isc); zz=>zzm210t(:,isc) + elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.1).and.(npaire.eq.1)) then + gg=>ggm211; hh=>hhm211; ss=>ssm211; rr=>rrm211; vv=>vvm211; ww=>wwm211; zz=>zzm211 + elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.1).and.(npaire.eq.1)) then + gg=>ggm211t(:,isc); hh=>hhm211t(:,isc); ss=>ssm211t(:,isc); rr=>rrm211t(:,isc); vv=>vvm211t(:,isc); ww=>wwm211t(:,isc); zz=>zzm211t(:,isc) + elseif(present(mhdvar)) then + if (nclyB1(mhdvar).eq.2 .and. nclyBn(mhdvar).eq.2) then + gg=>ggmB(:,mhdvar); hh=>hhmB(:,mhdvar); ss=>ssmB(:,mhdvar); rr=>rrmB(:,mhdvar); vv=>vvmB(:,mhdvar); ww=>wwmB(:,mhdvar); + zz=>zzmB(:,mhdvar) + endif + else + ! We should not be here + if (nrank == 0) then + write(*,*) "Error for time-implicit Y diffusion." + if (isc == 0) write(*,*) " Wrong combination for ncly1, nclyn and npaire", ncly1, nclyn, npaire + if (isc > 0) write(*,*) " Wrong combination for nclyS1, nclySn and npaire", nclyS1, nclySn, npaire + endif + call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop + endif + tb2=0.; + if ((isc.eq.0.and.ncly1.eq.0.and.nclyn.eq.0).or.(isc.gt.0.and.nclyS1.eq.0.and.nclySn.eq.0)) then + call septinv(tb2,ta2,gg,hh,ss,rr,vv,ww,zz,lo1,lo2,lo3,up1,up2,up3,ysize(1),ysize(2),ysize(3)) + nullify(lo1,lo2,lo3,up1,up2,up3) + else + call septinv(tb2,ta2,gg,hh,ss,rr,vv,ww,zz,ysize(1),ysize(2),ysize(3)) + endif + nullify(gg,hh,ss,rr,vv,ww,zz) + else if (isecondder.eq.5) then + tb2=0.; + !TO BE DONE: Different types of BC + if ((ncly1.eq.0).and.(nclyn.eq.0)) then + !NOT READY + elseif ((ncly1.eq.1).and.(nclyn.eq.1)) then + !NOT READY + elseif ((ncly1.eq.2).and.(nclyn.eq.2)) then + call nonainv(tb2,ta2,ggm,hhm,ssm,sssm,ttm,zzzm,zzm,wwm,vvm,ysize(1),ysize(2),ysize(3)) + endif + endif + + call transpose_y_to_x(tb2,var1) + + if (present(forcing1)) then + if ( (irestart.eq.1).or.(itime.gt.1) ) then + var1(:,:,:)=var1(:,:,:)+forcing1(:,:,:) + endif + endif + + + return + end subroutine inttimp - ! - ! Apply boundary conditions - ! - if ((isc.eq.0.and.ncly1.eq.2).or.(isc.gt.0.and.nclyS1.eq.2)) then - ta2(:,1,:) = bcbot(:,:) - endif - if ((isc.eq.0.and.nclyn.eq.2).or.(isc.gt.0.and.nclySn.eq.2)) then - ta2(:,ny,:) = bctop(:,:) - endif - - !Inversion of the linear system Mx=b: (A-xcst.B)u^n+1=uhat+(A+xcst.B)u^n - !if secondder=5, we need nona inversion - !if isecondder is not 5, we need septa inversion - - if (isecondder.ne.5) then - if ((isc.eq.0).and.(ncly1.eq.0).and.(nclyn.eq.0)) then - gg=>ggm0; hh=>hhm0; ss=>ssm0; rr=>rrm0; vv=>vvm0; ww=>wwm0; zz=>zzm0 - lo1=>l1m; lo2=>l2m; lo3=>l3m; up1=>u1m; up2=>u2m; up3=>u3m - elseif ((isc.gt.0).and.(nclyS1.eq.0).and.(nclySn.eq.0)) then - gg=>ggm0t(:,isc); hh=>hhm0t(:,isc); ss=>ssm0t(:,isc); rr=>rrm0t(:,isc); vv=>vvm0t(:,isc); ww=>wwm0t(:,isc); zz=>zzm0t(:,isc) - lo1=>l1mt(:,isc); lo2=>l2mt(:,isc); lo3=>l3mt(:,isc); up1=>u1mt(:,isc); up2=>u2mt(:,isc); up3=>u3mt(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.1).and.(npaire.eq.0)) then - gg=>ggm10; hh=>hhm10; ss=>ssm10; rr=>rrm10; vv=>vvm10; ww=>wwm10; zz=>zzm10 - elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.1).and.(npaire.eq.0)) then - gg=>ggm10t(:,isc); hh=>hhm10t(:,isc); ss=>ssm10t(:,isc); rr=>rrm10t(:,isc); vv=>vvm10t(:,isc); ww=>wwm10t(:,isc); zz=>zzm10t(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.1).and.(npaire.eq.1)) then - gg=>ggm11; hh=>hhm11; ss=>ssm11; rr=>rrm11; vv=>vvm11; ww=>wwm11; zz=>zzm11 - elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.1).and.(npaire.eq.1)) then - gg=>ggm11t(:,isc); hh=>hhm11t(:,isc); ss=>ssm11t(:,isc); rr=>rrm11t(:,isc); vv=>vvm11t(:,isc); ww=>wwm11t(:,isc); zz=>zzm11t(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.2)) then - gg=>ggm; hh=>hhm; ss=>ssm; rr=>rrm; vv=>vvm; ww=>wwm; zz=>zzm - elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.2)) then - gg=>ggmt(:,isc); hh=>hhmt(:,isc); ss=>ssmt(:,isc); rr=>rrmt(:,isc); vv=>vvmt(:,isc); ww=>wwmt(:,isc); zz=>zzmt(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.2).and.(npaire.eq.0)) then - gg=>ggm120; hh=>hhm120; ss=>ssm120; rr=>rrm120; vv=>vvm120; ww=>wwm120; zz=>zzm120 - elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.2).and.(npaire.eq.0)) then - gg=>ggm120t(:,isc); hh=>hhm120t(:,isc); ss=>ssm120t(:,isc); rr=>rrm120t(:,isc); vv=>vvm120t(:,isc); ww=>wwm120t(:,isc); zz=>zzm120t(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.1).and.(nclyn.eq.2).and.(npaire.eq.1)) then - gg=>ggm121; hh=>hhm121; ss=>ssm121; rr=>rrm121; vv=>vvm121; ww=>wwm121; zz=>zzm121 - elseif ((isc.gt.0).and.(nclyS1.eq.1).and.(nclySn.eq.2).and.(npaire.eq.1)) then - gg=>ggm121t(:,isc); hh=>hhm121t(:,isc); ss=>ssm121t(:,isc); rr=>rrm121t(:,isc); vv=>vvm121t(:,isc); ww=>wwm121t(:,isc); zz=>zzm121t(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.1).and.(npaire.eq.0)) then - gg=>ggm210; hh=>hhm210; ss=>ssm210; rr=>rrm210; vv=>vvm210; ww=>wwm210; zz=>zzm210 - elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.1).and.(npaire.eq.0)) then - gg=>ggm210t(:,isc); hh=>hhm210t(:,isc); ss=>ssm210t(:,isc); rr=>rrm210t(:,isc); vv=>vvm210t(:,isc); ww=>wwm210t(:,isc); zz=>zzm210t(:,isc) - elseif ((isc.eq.0).and.(ncly1.eq.2).and.(nclyn.eq.1).and.(npaire.eq.1)) then - gg=>ggm211; hh=>hhm211; ss=>ssm211; rr=>rrm211; vv=>vvm211; ww=>wwm211; zz=>zzm211 - elseif ((isc.gt.0).and.(nclyS1.eq.2).and.(nclySn.eq.1).and.(npaire.eq.1)) then - gg=>ggm211t(:,isc); hh=>hhm211t(:,isc); ss=>ssm211t(:,isc); rr=>rrm211t(:,isc); vv=>vvm211t(:,isc); ww=>wwm211t(:,isc); zz=>zzm211t(:,isc) - else - ! We should not be here - if (nrank == 0) then - write(*,*) "Error for time-implicit Y diffusion." - if (isc == 0) write(*,*) " Wrong combination for ncly1, nclyn and npaire", ncly1, nclyn, npaire - if (isc /= 0) write(*,*) " Wrong combination for nclyS1, nclySn and npaire", nclyS1, nclySn, npaire - endif - call MPI_ABORT(MPI_COMM_WORLD,code,ierror); stop - endif - tb2=0.; - if ((isc.eq.0.and.ncly1.eq.0.and.nclyn.eq.0).or.(isc.gt.0.and.nclyS1.eq.0.and.nclySn.eq.0)) then - call septinv(tb2,ta2,gg,hh,ss,rr,vv,ww,zz,lo1,lo2,lo3,up1,up2,up3,ysize(1),ysize(2),ysize(3)) - nullify(lo1,lo2,lo3,up1,up2,up3) - else - call septinv(tb2,ta2,gg,hh,ss,rr,vv,ww,zz,ysize(1),ysize(2),ysize(3)) - endif - nullify(gg,hh,ss,rr,vv,ww,zz) - else if (isecondder.eq.5) then - tb2=0.; - !TO BE DONE: Different types of BC - if ((ncly1.eq.0).and.(nclyn.eq.0)) then - !NOT READY - elseif ((ncly1.eq.1).and.(nclyn.eq.1)) then - !NOT READY - elseif ((ncly1.eq.2).and.(nclyn.eq.2)) then - call nonainv(tb2,ta2,ggm,hhm,ssm,sssm,ttm,zzzm,zzm,wwm,vvm,ysize(1),ysize(2),ysize(3)) - endif - endif - - call transpose_y_to_x(tb2,var1) - - if (present(forcing1)) then - if ( (irestart.eq.1).or.(itime.gt.1) ) then - var1(:,:,:)=var1(:,:,:)+forcing1(:,:,:) - endif - endif - - - return -end subroutine inttimp - -!******************************************************************** -!PREMULT FOR INTTIMP WITH SEPTADIAG -subroutine multmatrix7(td2,ta2,ux2,npaire,cly1,clyn,xcst) -! -!******************************************************************** - USE param, only : iimplicit, istret, zero, ncly1, nclyn, two - USE variables - USE derivY - USE decomp_2d - USE ibm_param, only : ubcx,ubcy,ubcz - - implicit none - - integer, intent(in) :: npaire, cly1, clyn - real(mytype), intent(in) :: xcst - integer :: i,j,k,code - real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: ux2 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: td2,ta2 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: di2 - - ! Compute A.ta2, store it in ta2 - if (istret.ne.0) then - do j=1,ysize(2) - ta2(:,j,:)=ta2(:,j,:)/pp2y(j) - enddo - endif - - ! j=1,2,3 - if (cly1.eq.0) then - td2(:,1,:) = alsajy*ta2(:,2,:) + ta2(:,1,:) + alsajy*ta2(:,ysize(2),:) - do j=2,3 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - else if (cly1.eq.1) then - if (npaire.eq.0) then - td2(:,1,:) = ta2(:,1,:) - else - td2(:,1,:) = two*alsajy*ta2(:,2,:) + ta2(:,1,:) - endif - do j=2,3 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - else - td2(:,1,:) = zero - td2(:,2,:) = alsa2y*ta2(:,1,:) + ta2(:,2,:) + alsa2y*ta2(:,3,:) - td2(:,3,:) = alsa3y*ta2(:,2,:) + ta2(:,3,:) + alsa3y*ta2(:,4,:) - endif - - do j=4,ysize(2)-3 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - - ! j=ny-2,ny-1,ny - if (clyn.eq.0) then - do j=ysize(2)-2,ysize(2)-1 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - td2(:,ysize(2),:) = alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) + alsajy*ta2(:,1,:) - elseif (clyn.eq.1) then - do j=ysize(2)-2,ysize(2)-1 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - if (npaire.eq.0) then - td2(:,ysize(2),:) = ta2(:,ysize(2),:) - else - td2(:,ysize(2),:) = two*alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) - endif - else - td2(:,ysize(2)-2,:) = alsaty*ta2(:,ysize(2)-3,:) + ta2(:,ysize(2)-2,:) + alsaty*ta2(:,ysize(2)-1,:) - td2(:,ysize(2)-1,:) = alsamy*ta2(:,ysize(2)-2,:) + ta2(:,ysize(2)-1,:) + alsamy*ta2(:,ysize(2),:) - td2(:,ysize(2),:) = zero - endif - ta2 = td2 - -! Compute (A+nu*dt.B).un -! nu*dt*B.un first, needed for CN, not needed for backward Euler - if (iimplicit.eq.1) then - - td2(:,:,:) = zero - - else if ((cly1.eq.1.or.clyn.eq.1) .and. npaire.eq.0) then - - ! Check if we are solving momentum or scalars - if (cly1.eq.ncly1 .and. clyn.eq.nclyn) then - call deryy(td2,ux2,di2,sy,sfy,ssy,swy,ysize(1),ysize(2),ysize(3),0, ubcx) - else - call deryyS(td2,ux2,di2,sy,sfyS,ssyS,swyS,ysize(1),ysize(2),ysize(3),0, ubcx) - endif - else - ! Check if we are solving momentum or scalars - if (cly1.eq.ncly1 .and. clyn.eq.nclyn) then - call deryy(td2,ux2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, ubcx) - else - call deryyS(td2,ux2,di2,sy,sfypS,ssypS,swypS,ysize(1),ysize(2),ysize(3),1, ubcx) - endif - - endif - td2(:,:,:) = xcst * td2(:,:,:) - - if (istret.ne.0) then - do j=1,ysize(2) - ux2(:,j,:)=ux2(:,j,:)/pp2y(j) - enddo - endif - - ! j=1,2,3 - if (cly1.eq.0) then - td2(:,1,:) = alsajy*ux2(:,2,:) + ux2(:,1,:) + alsajy*ux2(:,ysize(2),:) & - + td2(:,1,:) - do j=2,3 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - else if (cly1.eq.1) then - if (npaire.eq.0) then - td2(:,1,:) = ux2(:,1,:) + td2(:,1,:) - else - td2(:,1,:) = two*alsajy*ux2(:,2,:) + ux2(:,1,:) & - + td2(:,1,:) - endif - do j=2,3 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - else - td2(:,1,:) = zero - td2(:,2,:) = alsa2y*ux2(:,1,:) + ux2(:,2,:) + alsa2y*ux2(:,3,:) & - + td2(:,2,:) - td2(:,3,:) = alsa3y*ux2(:,2,:) + ux2(:,3,:) + alsa3y*ux2(:,4,:) & - + td2(:,3,:) - endif - - do j=4,ysize(2)-3 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - - ! j=ny-2,ny-1,ny - if (clyn.eq.0) then - do j=ysize(2)-2,ysize(2)-1 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - td2(:,ysize(2),:) = alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) + alsajy*ux2(:,1,:) & - + td2(:,ysize(2),:) - elseif (clyn.eq.1) then - do j=ysize(2)-2,ysize(2)-1 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - if (npaire.eq.0) then - td2(:,ysize(2),:) = ux2(:,ysize(2),:) + td2(:,ysize(2),:) - else - td2(:,ysize(2),:) = two*alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) & - + td2(:,ysize(2),:) - endif - else - td2(:,ysize(2)-2,:) = alsaty*ux2(:,ysize(2)-3,:) + ux2(:,ysize(2)-2,:) + alsaty*ux2(:,ysize(2)-1,:) & - + td2(:,ysize(2)-2,:) - td2(:,ysize(2)-1,:) = alsamy*ux2(:,ysize(2)-2,:) + ux2(:,ysize(2)-1,:) + alsamy*ux2(:,ysize(2),:) & - + td2(:,ysize(2)-1,:) - td2(:,ysize(2),:) = zero - endif - -end subroutine multmatrix7 - -!******************************************************************** - -subroutine multmatrix9(td2,ta2,ux2,npaire) - ! !******************************************************************** - USE param - USE variables - USE derivY - USE decomp_2d - USE ibm_param, only : ubcx,ubcy,ubcz - - implicit none + !PREMULT FOR INTTIMP WITH SEPTADIAG + subroutine multmatrix7(td2,ta2,ux2,npaire,cly1,clyn,xcst) + ! + !******************************************************************** + USE param, only : iimplicit, istret, zero, ncly1, nclyn, two + USE variables + USE derivY + USE decomp_2d + USE ibm_param, only : ubcx,ubcy,ubcz - integer, intent(in) :: npaire - integer :: i,j,k,code - real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: ux2 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: td2,ta2 - real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: di2 + implicit none + + integer, intent(in) :: npaire, cly1, clyn + real(mytype), intent(in) :: xcst + integer :: i,j,k,code + real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: ux2 + real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: td2,ta2 + real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: di2 + + ! Compute A.ta2, store it in ta2 + if (istret.ne.0) then + do j=1,ysize(2) + ta2(:,j,:)=ta2(:,j,:)/pp2y(j) + enddo + endif + + ! j=1,2,3 + if (cly1.eq.0) then + td2(:,1,:) = alsajy*ta2(:,2,:) + ta2(:,1,:) + alsajy*ta2(:,ysize(2),:) + do j=2,3 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + else if (cly1.eq.1) then + if (npaire.eq.0) then + td2(:,1,:) = ta2(:,1,:) + else + td2(:,1,:) = two*alsajy*ta2(:,2,:) + ta2(:,1,:) + endif + do j=2,3 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + else + td2(:,1,:) = zero + td2(:,2,:) = alsa2y*ta2(:,1,:) + ta2(:,2,:) + alsa2y*ta2(:,3,:) + td2(:,3,:) = alsa3y*ta2(:,2,:) + ta2(:,3,:) + alsa3y*ta2(:,4,:) + endif + + do j=4,ysize(2)-3 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + + ! j=ny-2,ny-1,ny + if (clyn.eq.0) then + do j=ysize(2)-2,ysize(2)-1 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + td2(:,ysize(2),:) = alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) + alsajy*ta2(:,1,:) + elseif (clyn.eq.1) then + do j=ysize(2)-2,ysize(2)-1 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + if (npaire.eq.0) then + td2(:,ysize(2),:) = ta2(:,ysize(2),:) + else + td2(:,ysize(2),:) = two*alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) + endif + else + td2(:,ysize(2)-2,:) = alsaty*ta2(:,ysize(2)-3,:) + ta2(:,ysize(2)-2,:) + alsaty*ta2(:,ysize(2)-1,:) + td2(:,ysize(2)-1,:) = alsamy*ta2(:,ysize(2)-2,:) + ta2(:,ysize(2)-1,:) + alsamy*ta2(:,ysize(2),:) + td2(:,ysize(2),:) = zero + endif + ta2 = td2 + + ! Compute (A+nu*dt.B).un + ! nu*dt*B.un first, needed for CN, not needed for backward Euler + if (iimplicit.eq.1) then + + td2(:,:,:) = zero + + else if ((cly1.eq.1.or.clyn.eq.1) .and. npaire.eq.0) then + + ! Check if we are solving momentum or scalars + if (cly1.eq.ncly1 .and. clyn.eq.nclyn) then + call deryy(td2,ux2,di2,sy,sfy,ssy,swy,ysize(1),ysize(2),ysize(3),0, ubcx) + else + call deryyS(td2,ux2,di2,sy,sfyS,ssyS,swyS,ysize(1),ysize(2),ysize(3),0, ubcx) + endif + else + ! Check if we are solving momentum or scalars + if (cly1.eq.ncly1 .and. clyn.eq.nclyn) then + call deryy(td2,ux2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, ubcx) + else + call deryyS(td2,ux2,di2,sy,sfypS,ssypS,swypS,ysize(1),ysize(2),ysize(3),1, ubcx) + endif + + endif + td2(:,:,:) = xcst * td2(:,:,:) + + if (istret.ne.0) then + do j=1,ysize(2) + ux2(:,j,:)=ux2(:,j,:)/pp2y(j) + enddo + endif + + ! j=1,2,3 + if (cly1.eq.0) then + td2(:,1,:) = alsajy*ux2(:,2,:) + ux2(:,1,:) + alsajy*ux2(:,ysize(2),:) & + + td2(:,1,:) + do j=2,3 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + else if (cly1.eq.1) then + if (npaire.eq.0) then + td2(:,1,:) = ux2(:,1,:) + td2(:,1,:) + else + td2(:,1,:) = two*alsajy*ux2(:,2,:) + ux2(:,1,:) & + + td2(:,1,:) + endif + do j=2,3 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + else + td2(:,1,:) = zero + td2(:,2,:) = alsa2y*ux2(:,1,:) + ux2(:,2,:) + alsa2y*ux2(:,3,:) & + + td2(:,2,:) + td2(:,3,:) = alsa3y*ux2(:,2,:) + ux2(:,3,:) + alsa3y*ux2(:,4,:) & + + td2(:,3,:) + endif + + do j=4,ysize(2)-3 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + + ! j=ny-2,ny-1,ny + if (clyn.eq.0) then + do j=ysize(2)-2,ysize(2)-1 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + td2(:,ysize(2),:) = alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) + alsajy*ux2(:,1,:) & + + td2(:,ysize(2),:) + elseif (clyn.eq.1) then + do j=ysize(2)-2,ysize(2)-1 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + if (npaire.eq.0) then + td2(:,ysize(2),:) = ux2(:,ysize(2),:) + td2(:,ysize(2),:) + else + td2(:,ysize(2),:) = two*alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) & + + td2(:,ysize(2),:) + endif + else + td2(:,ysize(2)-2,:) = alsaty*ux2(:,ysize(2)-3,:) + ux2(:,ysize(2)-2,:) + alsaty*ux2(:,ysize(2)-1,:) & + + td2(:,ysize(2)-2,:) + td2(:,ysize(2)-1,:) = alsamy*ux2(:,ysize(2)-2,:) + ux2(:,ysize(2)-1,:) + alsamy*ux2(:,ysize(2),:) & + + td2(:,ysize(2)-1,:) + td2(:,ysize(2),:) = zero + endif + + end subroutine multmatrix7 + + !******************************************************************** + subroutine multmatrix9(td2,ta2,ux2,npaire) + ! + !******************************************************************** + USE param + USE variables + USE derivY + USE decomp_2d + USE ibm_param, only : ubcx,ubcy,ubcz - !A.uhat + implicit none - if (istret.ne.0) then - do j=1,ysize(2) - ta2(:,j,:)=ta2(:,j,:)/pp2y(j) - enddo - endif + integer, intent(in) :: npaire + integer :: i,j,k,code + real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: ux2 + real(mytype),dimension(ysize(1),ysize(2),ysize(3)), intent(inout) :: td2,ta2 + real(mytype),dimension(ysize(1),ysize(2),ysize(3)) :: di2 - if (ncly1.eq.0) then - td2(:,1,:) = alsajy*ta2(:,2,:) + ta2(:,1,:) + alsajy*ta2(:,ysize(2),:) - do j=2,ysize(2)-1 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - td2(:,ysize(2),:) = alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) + alsajy*ta2(:,1,:) - ta2=td2 + !A.uhat - elseif (ncly1.eq.1) then + if (istret.ne.0) then + do j=1,ysize(2) + ta2(:,j,:)=ta2(:,j,:)/pp2y(j) + enddo + endif + if (ncly1.eq.0) then - elseif (ncly1.eq.2) then + td2(:,1,:) = alsajy*ta2(:,2,:) + ta2(:,1,:) + alsajy*ta2(:,ysize(2),:) + do j=2,ysize(2)-1 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + td2(:,ysize(2),:) = alsajy*ta2(:,ysize(2)-1,:) + ta2(:,ysize(2),:) + alsajy*ta2(:,1,:) + ta2=td2 - td2(:,1,:) = zero - td2(:,2,:) = alsa2y*ta2(:,1,:) + ta2(:,2,:) + alsa2y*ta2(:,3,:) - td2(:,3,:) = alsa3y*ta2(:,2,:) + ta2(:,3,:) + alsa3y*ta2(:,4,:) - td2(:,4,:) = alsa4y*ta2(:,3,:) + ta2(:,4,:) + alsa4y*ta2(:,5,:) - do j=5,ysize(2)-4 - td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) - enddo - td2(:,ysize(2)-3,:) = alsatty*ta2(:,ysize(2)-4,:) + ta2(:,ysize(2)-3,:) + alsatty*ta2(:,ysize(2)-2,:) - td2(:,ysize(2)-2,:) = alsaty*ta2(:,ysize(2)-3,:) + ta2(:,ysize(2)-2,:) + alsaty*ta2(:,ysize(2)-1,:) - td2(:,ysize(2)-1,:) = alsamy*ta2(:,ysize(2)-2,:) + ta2(:,ysize(2)-1,:) + alsamy*ta2(:,ysize(2),:) - td2(:,ysize(2),:) = 0. - ta2=td2 + elseif (ncly1.eq.1) then - endif - !(A+nu*dt.B).un + elseif (ncly1.eq.2) then - if (iimplicit.eq.1) then + td2(:,1,:) = zero + td2(:,2,:) = alsa2y*ta2(:,1,:) + ta2(:,2,:) + alsa2y*ta2(:,3,:) + td2(:,3,:) = alsa3y*ta2(:,2,:) + ta2(:,3,:) + alsa3y*ta2(:,4,:) + td2(:,4,:) = alsa4y*ta2(:,3,:) + ta2(:,4,:) + alsa4y*ta2(:,5,:) + do j=5,ysize(2)-4 + td2(:,j,:) = alsajy*ta2(:,j-1,:) + ta2(:,j,:) + alsajy*ta2(:,j+1,:) + enddo + td2(:,ysize(2)-3,:) = alsatty*ta2(:,ysize(2)-4,:) + ta2(:,ysize(2)-3,:) + alsatty*ta2(:,ysize(2)-2,:) + td2(:,ysize(2)-2,:) = alsaty*ta2(:,ysize(2)-3,:) + ta2(:,ysize(2)-2,:) + alsaty*ta2(:,ysize(2)-1,:) + td2(:,ysize(2)-1,:) = alsamy*ta2(:,ysize(2)-2,:) + ta2(:,ysize(2)-1,:) + alsamy*ta2(:,ysize(2),:) + td2(:,ysize(2),:) = 0. + ta2=td2 - td2(:,:,:) = zero + endif - elseif ((ncly1.eq.1.or.nclyn.eq.1) .and. npaire.eq.0) then + !(A+nu*dt.B).un - call deryy(td2,ux2,di2,sy,sfy,ssy,swy,ysize(1),ysize(2),ysize(3),0, ubcx) - else - call deryy(td2,ux2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, ubcx) + if (iimplicit.eq.1) then - endif - td2(:,:,:) = xcst * td2(:,:,:) + td2(:,:,:) = zero - if (istret.ne.0) then - do j=1,ysize(2) - ux2(:,j,:)=ux2(:,j,:)/pp2y(j) - enddo - endif + elseif ((ncly1.eq.1.or.nclyn.eq.1) .and. npaire.eq.0) then - if (ncly1.eq.0) then + call deryy(td2,ux2,di2,sy,sfy,ssy,swy,ysize(1),ysize(2),ysize(3),0, ubcx) + else + call deryy(td2,ux2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1, ubcx) - td2(:,1,:) = alsajy*ux2(:,2,:) + ux2(:,1,:) + alsajy*ux2(:,ysize(2),:) & - + td2(:,1,:) - do j=2,ysize(2)-1 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - td2(:,ysize(2),:) = alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) + alsajy*ux2(:,1,:) & - + td2(:,ysize(2),:) + endif + td2(:,:,:) = xcst * td2(:,:,:) - elseif (ncly1.eq.1) then + if (istret.ne.0) then + do j=1,ysize(2) + ux2(:,j,:)=ux2(:,j,:)/pp2y(j) + enddo + endif - elseif (ncly1.eq.2) then + if (ncly1.eq.0) then - td2(:,1,:) = zero - td2(:,2,:) = alsa2y*ux2(:,1,:) + ux2(:,2,:) + alsa2y*ux2(:,3,:) & - + td2(:,2,:) - td2(:,3,:) = alsa3y*ux2(:,2,:) + ux2(:,3,:) + alsa3y*ux2(:,4,:) & - + td2(:,3,:) - td2(:,4,:) = alsa4y*ux2(:,3,:) + ux2(:,4,:) + alsa4y*ux2(:,5,:) & - + td2(:,4,:) - do j=5,ysize(2)-4 - td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & - + td2(:,j,:) - enddo - td2(:,ysize(2)-3,:) = alsatty*ux2(:,ysize(2)-4,:) + ux2(:,ysize(2)-3,:) + alsatty*ux2(:,ysize(2)-2,:) & - + td2(:,ysize(2)-3,:) - td2(:,ysize(2)-2,:) = alsaty*ux2(:,ysize(2)-3,:) + ux2(:,ysize(2)-2,:) + alsaty*ux2(:,ysize(2)-1,:) & - + td2(:,ysize(2)-2,:) - td2(:,ysize(2)-1,:) = alsamy*ux2(:,ysize(2)-2,:) + ux2(:,ysize(2)-1,:) + alsamy*ux2(:,ysize(2),:) & - + td2(:,ysize(2)-1,:) - td2(:,ysize(2),:) = zero + td2(:,1,:) = alsajy*ux2(:,2,:) + ux2(:,1,:) + alsajy*ux2(:,ysize(2),:) & + + td2(:,1,:) + do j=2,ysize(2)-1 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + td2(:,ysize(2),:) = alsajy*ux2(:,ysize(2)-1,:) + ux2(:,ysize(2),:) + alsajy*ux2(:,1,:) & + + td2(:,ysize(2),:) + + elseif (ncly1.eq.1) then + + elseif (ncly1.eq.2) then + + td2(:,1,:) = zero + td2(:,2,:) = alsa2y*ux2(:,1,:) + ux2(:,2,:) + alsa2y*ux2(:,3,:) & + + td2(:,2,:) + td2(:,3,:) = alsa3y*ux2(:,2,:) + ux2(:,3,:) + alsa3y*ux2(:,4,:) & + + td2(:,3,:) + td2(:,4,:) = alsa4y*ux2(:,3,:) + ux2(:,4,:) + alsa4y*ux2(:,5,:) & + + td2(:,4,:) + do j=5,ysize(2)-4 + td2(:,j,:) = alsajy*ux2(:,j-1,:) + ux2(:,j,:) + alsajy*ux2(:,j+1,:) & + + td2(:,j,:) + enddo + td2(:,ysize(2)-3,:) = alsatty*ux2(:,ysize(2)-4,:) + ux2(:,ysize(2)-3,:) + alsatty*ux2(:,ysize(2)-2,:) & + + td2(:,ysize(2)-3,:) + td2(:,ysize(2)-2,:) = alsaty*ux2(:,ysize(2)-3,:) + ux2(:,ysize(2)-2,:) + alsaty*ux2(:,ysize(2)-1,:) & + + td2(:,ysize(2)-2,:) + td2(:,ysize(2)-1,:) = alsamy*ux2(:,ysize(2)-2,:) + ux2(:,ysize(2)-1,:) + alsamy*ux2(:,ysize(2),:) & + + td2(:,ysize(2)-1,:) + td2(:,ysize(2),:) = zero - endif + endif -end subroutine multmatrix9 + end subroutine multmatrix9 -! -! Compute 1D arrays containing LU decomposition -! -subroutine implicit_schemes() + ! + ! Compute 1D arrays containing LU decomposition + ! + subroutine implicit_schemes() - USE param - USE derivY - USE variables - USE var - USE param - use ludecomp + USE param + USE derivY + USE variables + USE var + USE param + use ludecomp - implicit none + implicit none - integer :: i,j,k,is + integer :: i,j,k,is !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !! MATRIX M2 TIME IMPLICIT !! + !! MATRIX M2 TIME IMPLICIT !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! -!FOR VELOCITY !!!!!!! + !FOR VELOCITY !!!!!!! !!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!! -if (isecondder.ne.5) then + if (isecondder.ne.5) then !!!!!!!!!!!!!!!!!!!!!! - !7-DIAG !! + !7-DIAG !! !!!!!!!!!!!!!!!!!!!!!! !!! NCL = 2, dirichlet BC - ! - !DIAG - aam(1 )=as1y - aam(ny )=asny - aam(2 )=-two*as2y - aam(ny-1 )=-two*asmy - aam(3 )=-two*(as3y+bs3y) - aam(ny-2 )=-two*(asty+bsty) - aam(4:ny-3)=-two*(asjy+bsjy+csjy) - call init_implicit_coef(aam, aamt) - if (istret==0) then - aam = one-xcst*aam - do is = 1, numscalar - aamt(:,is) = one-xcst_sc(is)*aamt(:,is) - enddo - else - aam = one/pp2y -xcst*aam - do is = 1, numscalar - aamt(:,is) = one/pp2y -xcst_sc(is)*aamt(:,is) - enddo - endif - ! BC on aam, dirichlet - aam(1 )=one - aam(ny)=one - ! BC on aamt - if (istret.eq.0) then - do is = 1, numscalar - aamt(1 ,is) = alpha_sc(is,1) + beta_sc(is,1)*(11.d0/6.d0/dy) - aamt(ny,is) = alpha_sc(is,2) + beta_sc(is,2)*(11.d0/6.d0/dy) - enddo - else - do is = 1, numscalar - aamt(1 ,is) = alpha_sc(is,1) + beta_sc(is,1)*ppy(1 )*(11.d0/6.d0/dy) - aamt(ny,is) = alpha_sc(is,2) + beta_sc(is,2)*ppy(ny)*(11.d0/6.d0/dy) - enddo - endif - !DIAG SUP 1 - bbm(1 )=bs1y - bbm(ny )=bsny - bbm(2 )=as2y - bbm(ny-1 )=asmy - bbm(3 )=as3y - bbm(ny-2 )=asty - bbm(4:ny-3)=asjy - call init_implicit_coef(bbm, bbmt) - bbm = -xcst*bbm - do is = 1, numscalar - bbmt(:,is) = -xcst_sc(is)*bbmt(:,is) - enddo - if (istret==0) then - bbm(2 )=bbm(2 )+alsa2y - bbm(ny-1 )=bbm(ny-1 )+alsamy - bbm(3 )=bbm(3 )+alsa3y - bbm(ny-2 )=bbm(ny-2 )+alsaty - bbm(4:ny-3)=bbm(4:ny-3)+alsajy - do is = 1, numscalar - bbmt(2 ,is)=bbmt(2 ,is)+alsa2y - bbmt(ny-1 ,is)=bbmt(ny-1 ,is)+alsamy - bbmt(3 ,is)=bbmt(3 ,is)+alsa3y - bbmt(ny-2 ,is)=bbmt(ny-2 ,is)+alsaty - bbmt(4:ny-3,is)=bbmt(4:ny-3,is)+alsajy - enddo - else - bbm(2 )=bbm(2 )+alsa2y/pp2y(3) - bbm(ny-1 )=bbm(ny-1 )+alsamy/pp2y(ny) - bbm(3 )=bbm(3 )+alsa3y/pp2y(4) - bbm(ny-2 )=bbm(ny-2 )+alsaty/pp2y(ny-1) - bbm(4:ny-3)=bbm(4:ny-3)+alsajy/pp2y(5:ny-2) - do is = 1, numscalar - bbmt(2 ,is)=bbmt(2 ,is)+alsa2y/pp2y(3) - bbmt(ny-1 ,is)=bbmt(ny-1 ,is)+alsamy/pp2y(ny) - bbmt(3 ,is)=bbmt(3 ,is)+alsa3y/pp2y(4) - bbmt(ny-2 ,is)=bbmt(ny-2 ,is)+alsaty/pp2y(ny-1) - bbmt(4:ny-3,is)=bbmt(4:ny-3,is)+alsajy/pp2y(5:ny-2) - enddo - endif - ! BC on bbm, dirichlet - bbm(1 )=zero - bbm(ny)=zero - ! BC on bbmt - if (istret.eq.0) then - do is = 1, numscalar - bbmt(1 ,is) = beta_sc(is,1)*(-18.d0/6.d0/dy) - bbmt(ny,is) = zero - enddo - else - do is = 1, numscalar - bbmt(1 ,is) = beta_sc(is,1)*ppy(1)*(-18.d0/6.d0/dy) - bbmt(ny,is) = zero - enddo - endif - !DIAG SUP 2 - ccm(1 )=cs1y - ccm(ny )=csny - ccm(2 )=zero - ccm(ny-1 )=zero - ccm(3 )=bs3y - ccm(ny-2 )=bsty - ccm(4:ny-3)=bsjy - call init_implicit_coef(ccm, ccmt) - ccm = -xcst*ccm - do is = 1, numscalar - ccmt(:,is) = -xcst_sc(is)*ccmt(:,is) - enddo - ! BC on ccm - ccm(1 )=zero - ccm(ny)=zero - ! BC on ccmt - if (istret.eq.0) then - do is = 1, numscalar - ccmt(1 ,is) = beta_sc(is,1)*(nine/six/dy) - ccmt(ny,is) = zero - enddo - else - do is = 1, numscalar - ccmt(1 ,is) = beta_sc(is,1)*ppy(1)*(nine/six/dy) - ccmt(ny,is) = zero - enddo - endif - !DIAG SUP 3 - rrm(1 )=ds1y - rrm(ny )=dsny - rrm(2 )=zero - rrm(ny-1 )=zero - rrm(3 )=zero - rrm(ny-2 )=zero - rrm(4:ny-3)=csjy - call init_implicit_coef(rrm, rrmt) - rrm = -xcst*rrm - do is = 1, numscalar - rrmt(:,is) = -xcst_sc(is)*rrmt(:,is) - enddo - ! BC on rrm - rrm(1 )=zero - rrm(ny)=zero - ! BC on rrmt - if (istret.eq.0) then - do is = 1, numscalar - rrmt(1 ,is) = beta_sc(is,1)*(-two/six/dy) - rrmt(ny,is) = zero - enddo - else - do is = 1, numscalar - rrmt(1 ,is) = beta_sc(is,1)*ppy(1)*(-two/six/dy) - rrmt(ny,is) = zero - enddo - endif - !DIAG INF 1 - if (istret==0) then - ddm=bbm - ddmt = bbmt - else - ddm(1 )=bs1y - ddm(ny )=bsny - ddm(2 )=as2y - ddm(ny-1 )=asmy - ddm(3 )=as3y - ddm(ny-2 )=asty - ddm(4:ny-3)=asjy - call init_implicit_coef(ddm, ddmt) - ddm = -xcst*ddm - do is = 1, numscalar - ddmt(:,is) = -xcst_sc(is)*ddmt(:,is) - enddo - ddm(2 )=ddm(2 )+alsa2y/pp2y(1) - ddm(ny-1 )=ddm(ny-1 )+alsamy/pp2y(ny-2) - ddm(3 )=ddm(3 )+alsa3y/pp2y(2) - ddm(ny-2 )=ddm(ny-2 )+alsaty/pp2y(ny-3) - ddm(4:ny-3)=ddm(4:ny-3)+alsajy/pp2y(3:ny-4) - do is = 1, numscalar - ddmt(2 ,is)=ddmt(2 ,is)+alsa2y/pp2y(1) - ddmt(ny-1 ,is)=ddmt(ny-1 ,is)+alsamy/pp2y(ny-2) - ddmt(3 ,is)=ddmt(3 ,is)+alsa3y/pp2y(2) - ddmt(ny-2 ,is)=ddmt(ny-2 ,is)+alsaty/pp2y(ny-3) - ddmt(4:ny-3,is)=ddmt(4:ny-3,is)+alsajy/pp2y(3:ny-4) - enddo - ! BC on ddm - ddm(1 )=zero - ddm(ny)=zero - endif - ! BC on ddmt - if (istret.eq.0) then - do is = 1, numscalar - ddmt(1 ,is) = zero - ddmt(ny,is) = beta_sc(is,2)*(-eighteen/six/dy) - enddo - else - do is = 1, numscalar - ddmt(1 ,is) = zero - ddmt(ny,is) = beta_sc(is,2)*ppy(ny)*(-eighteen/six/dy) - enddo - endif - ! - !DIAG INF 2 - eem=ccm - eemt = ccmt - ! BC on eemt - if (istret.eq.0) then - do is = 1, numscalar - eemt(1 ,is) = zero - eemt(ny,is) = beta_sc(is,2)*(nine/six/dy) - enddo - else - do is = 1, numscalar - eemt(1 ,is) = zero - eemt(ny,is) = beta_sc(is,2)*ppy(ny)*(nine/six/dy) - enddo - endif - ! - !DIAG INF 3 - qqm=rrm - qqmt = rrmt - ! BC on rrmt - if (istret.eq.0) then - do is = 1, numscalar - qqmt(1 ,is) = zero - qqmt(ny,is) = beta_sc(is,2)*(-two/six/dy) - enddo - else - do is = 1, numscalar - qqmt(1 ,is) = zero - qqmt(ny,is) = beta_sc(is,2)*ppy(ny)*(-two/six/dy) - enddo - endif + ! + !DIAG + aam(1 )=as1y + aam(ny )=asny + aam(2 )=-two*as2y + aam(ny-1 )=-two*asmy + aam(3 )=-two*(as3y+bs3y) + aam(ny-2 )=-two*(asty+bsty) + aam(4:ny-3)=-two*(asjy+bsjy+csjy) + call init_implicit_coef(aam, aamt) + + do is=1, 3 + aamB(:,is) = aam + enddo + + if (istret==0) then + aam = one-xcst*aam + do is = 1, numscalar + aamt(:,is) = one-xcst_sc(is)*aamt(:,is) + enddo + do is = 1, 3 + aamB(:,is) = one-xcstB*aamB(:,is) + enddo + else + aam = one/pp2y -xcst*aam + do is = 1, numscalar + aamt(:,is) = one/pp2y -xcst_sc(is)*aamt(:,is) + enddo + do is = 1, 3 + aamB(:,is) = one/pp2y -xcstB*aamB(:,is) + enddo + endif + ! BC on aam, dirichlet + aam(1 )=one + aam(ny)=one + ! BC on aamt + if (istret.eq.0) then + do is = 1, numscalar + aamt(1 ,is) = alpha_sc(is,1) + beta_sc(is,1)*(11.d0/6.d0/dy) + aamt(ny,is) = alpha_sc(is,2) + beta_sc(is,2)*(11.d0/6.d0/dy) + enddo + else + do is = 1, numscalar + aamt(1 ,is) = alpha_sc(is,1) + beta_sc(is,1)*ppy(1 )*(11.d0/6.d0/dy) + aamt(ny,is) = alpha_sc(is,2) + beta_sc(is,2)*ppy(ny)*(11.d0/6.d0/dy) + enddo + endif + do is = 1, 3 + aamB(1,is)=one !todo -- not this for partally conducting? + aamB(ny,is)=one + enddo + + !DIAG SUP 1 + bbm(1 )=bs1y + bbm(ny )=bsny + bbm(2 )=as2y + bbm(ny-1 )=asmy + bbm(3 )=as3y + bbm(ny-2 )=asty + bbm(4:ny-3)=asjy + call init_implicit_coef(bbm, bbmt) + do is=1, 3 + bbmB(:,is) = bbm + enddo + + bbm = -xcst*bbm + do is = 1, numscalar + bbmt(:,is) = -xcst_sc(is)*bbmt(:,is) + enddo + do is = 1, 3 + bbmB(:,is) = -xcstB*bbmB(:,is) + enddo + + if (istret==0) then + bbm(2 )=bbm(2 )+alsa2y + bbm(ny-1 )=bbm(ny-1 )+alsamy + bbm(3 )=bbm(3 )+alsa3y + bbm(ny-2 )=bbm(ny-2 )+alsaty + bbm(4:ny-3)=bbm(4:ny-3)+alsajy + do is = 1, numscalar + bbmt(2 ,is)=bbmt(2 ,is)+alsa2y + bbmt(ny-1 ,is)=bbmt(ny-1 ,is)+alsamy + bbmt(3 ,is)=bbmt(3 ,is)+alsa3y + bbmt(ny-2 ,is)=bbmt(ny-2 ,is)+alsaty + bbmt(4:ny-3,is)=bbmt(4:ny-3,is)+alsajy + enddo + do is = 1, 3 + bbmB(2 ,is)=bbmB(2 ,is)+alsa2y + bbmB(ny-1 ,is)=bbmB(ny-1 ,is)+alsamy + bbmB(3 ,is)=bbmB(3 ,is)+alsa3y + bbmB(ny-2 ,is)=bbmB(ny-2 ,is)+alsaty + bbmB(4:ny-3,is)=bbmB(4:ny-3,is)+alsajy + enddo + + else + bbm(2 )=bbm(2 )+alsa2y/pp2y(3) + bbm(ny-1 )=bbm(ny-1 )+alsamy/pp2y(ny) + bbm(3 )=bbm(3 )+alsa3y/pp2y(4) + bbm(ny-2 )=bbm(ny-2 )+alsaty/pp2y(ny-1) + bbm(4:ny-3)=bbm(4:ny-3)+alsajy/pp2y(5:ny-2) + do is = 1, numscalar + bbmt(2 ,is)=bbmt(2 ,is)+alsa2y/pp2y(3) + bbmt(ny-1 ,is)=bbmt(ny-1 ,is)+alsamy/pp2y(ny) + bbmt(3 ,is)=bbmt(3 ,is)+alsa3y/pp2y(4) + bbmt(ny-2 ,is)=bbmt(ny-2 ,is)+alsaty/pp2y(ny-1) + bbmt(4:ny-3,is)=bbmt(4:ny-3,is)+alsajy/pp2y(5:ny-2) + enddo + do is = 1, 3 + bbmB(2 ,is)=bbmB(2 ,is)+alsa2y/pp2y(3) + bbmB(ny-1 ,is)=bbmB(ny-1 ,is)+alsamy/pp2y(ny) + bbmB(3 ,is)=bbmB(3 ,is)+alsa3y/pp2y(4) + bbmB(ny-2 ,is)=bbmB(ny-2 ,is)+alsaty/pp2y(ny-1) + bbmB(4:ny-3,is)=bbmB(4:ny-3,is)+alsajy/pp2y(5:ny-2) + enddo + + endif + ! BC on bbm, dirichlet + bbm(1 )=zero + bbm(ny)=zero + ! BC on bbmt + if (istret.eq.0) then + do is = 1, numscalar + bbmt(1 ,is) = beta_sc(is,1)*(-18.d0/6.d0/dy) + bbmt(ny,is) = zero + enddo + else + do is = 1, numscalar + bbmt(1 ,is) = beta_sc(is,1)*ppy(1)*(-18.d0/6.d0/dy) + bbmt(ny,is) = zero + enddo + endif + do is = 1, 3 + bbmB(1 ,is) = zero + bbmB(ny,is) = zero + enddo + + !DIAG SUP 2 + ccm(1 )=cs1y + ccm(ny )=csny + ccm(2 )=zero + ccm(ny-1 )=zero + ccm(3 )=bs3y + ccm(ny-2 )=bsty + ccm(4:ny-3)=bsjy + call init_implicit_coef(ccm, ccmt) + + do is=1, 3 + ccmB(:,is) = ccm + enddo + + ccm = -xcst*ccm + do is = 1, numscalar + ccmt(:,is) = -xcst_sc(is)*ccmt(:,is) + enddo + + do is=1, 3 + ccmB(:,is) = -xcstB*ccmB(:,is) + enddo + + ! BC on ccm + ccm(1 )=zero + ccm(ny)=zero + ! BC on ccmt + if (istret.eq.0) then + do is = 1, numscalar + ccmt(1 ,is) = beta_sc(is,1)*(nine/six/dy) + ccmt(ny,is) = zero + enddo + else + do is = 1, numscalar + ccmt(1 ,is) = beta_sc(is,1)*ppy(1)*(nine/six/dy) + ccmt(ny,is) = zero + enddo + endif + + do is=1, 3 + ccmB(1,is) = zero + ccmB(ny,is) = zero + enddo + + !DIAG SUP 3 + rrm(1 )=ds1y + rrm(ny )=dsny + rrm(2 )=zero + rrm(ny-1 )=zero + rrm(3 )=zero + rrm(ny-2 )=zero + rrm(4:ny-3)=csjy + call init_implicit_coef(rrm, rrmt) + + do is=1, 3 + rrmB(:,is) = rrm + enddo + + rrm = -xcst*rrm + do is = 1, numscalar + rrmt(:,is) = -xcst_sc(is)*rrmt(:,is) + enddo + do is=1, 3 + rrmB(:,is) = -xcstB*rrmB(:,is) + enddo + + ! BC on rrm + rrm(1 )=zero + rrm(ny)=zero + ! BC on rrmt + if (istret.eq.0) then + do is = 1, numscalar + rrmt(1 ,is) = beta_sc(is,1)*(-two/six/dy) + rrmt(ny,is) = zero + enddo + else + do is = 1, numscalar + rrmt(1 ,is) = beta_sc(is,1)*ppy(1)*(-two/six/dy) + rrmt(ny,is) = zero + enddo + endif + do is=1, 3 + rrmB(1,is) = zero + rrmB(ny,is) = zero + enddo + + + !DIAG INF 1 + if (istret==0) then + ddm=bbm + ddmt = bbmt + ddmB = bbmB + else + ddm(1 )=bs1y + ddm(ny )=bsny + ddm(2 )=as2y + ddm(ny-1 )=asmy + ddm(3 )=as3y + ddm(ny-2 )=asty + ddm(4:ny-3)=asjy + call init_implicit_coef(ddm, ddmt) + do is=1, 3 + ddmB(:,is) = ddm + enddo + + ddm = -xcst*ddm + do is = 1, numscalar + ddmt(:,is) = -xcst_sc(is)*ddmt(:,is) + enddo + do is=1, 3 + ddmB(:,is) = -xcstB*ddmB(:,is) + enddo + + ddm(2 )=ddm(2 )+alsa2y/pp2y(1) + ddm(ny-1 )=ddm(ny-1 )+alsamy/pp2y(ny-2) + ddm(3 )=ddm(3 )+alsa3y/pp2y(2) + ddm(ny-2 )=ddm(ny-2 )+alsaty/pp2y(ny-3) + ddm(4:ny-3)=ddm(4:ny-3)+alsajy/pp2y(3:ny-4) + do is = 1, numscalar + ddmt(2 ,is)=ddmt(2 ,is)+alsa2y/pp2y(1) + ddmt(ny-1 ,is)=ddmt(ny-1 ,is)+alsamy/pp2y(ny-2) + ddmt(3 ,is)=ddmt(3 ,is)+alsa3y/pp2y(2) + ddmt(ny-2 ,is)=ddmt(ny-2 ,is)+alsaty/pp2y(ny-3) + ddmt(4:ny-3,is)=ddmt(4:ny-3,is)+alsajy/pp2y(3:ny-4) + enddo + do is = 1, 3 + ddmB(2 ,is)=ddmB(2 ,is)+alsa2y/pp2y(1) + ddmB(ny-1 ,is)=ddmB(ny-1 ,is)+alsamy/pp2y(ny-2) + ddmB(3 ,is)=ddmB(3 ,is)+alsa3y/pp2y(2) + ddmB(ny-2 ,is)=ddmB(ny-2 ,is)+alsaty/pp2y(ny-3) + ddmB(4:ny-3,is)=ddmB(4:ny-3,is)+alsajy/pp2y(3:ny-4) + enddo + + ! BC on ddm + ddm(1 )=zero + ddm(ny)=zero + endif + ! BC on ddmt + if (istret.eq.0) then + do is = 1, numscalar + ddmt(1 ,is) = zero + ddmt(ny,is) = beta_sc(is,2)*(-eighteen/six/dy) + enddo + else + do is = 1, numscalar + ddmt(1 ,is) = zero + ddmt(ny,is) = beta_sc(is,2)*ppy(ny)*(-eighteen/six/dy) + enddo + endif + ddmB(1,:) = zero + ddmB(ny,:) = zero + + ! + !DIAG INF 2 + eem=ccm + eemt = ccmt + eemB = ccmB + ! BC on eemt + if (istret.eq.0) then + do is = 1, numscalar + eemt(1 ,is) = zero + eemt(ny,is) = beta_sc(is,2)*(nine/six/dy) + enddo + else + do is = 1, numscalar + eemt(1 ,is) = zero + eemt(ny,is) = beta_sc(is,2)*ppy(ny)*(nine/six/dy) + enddo + endif + do is = 1, 3 + eemB(1,is) = zero + eemB(ny,is) = zero + enddo + ! + !DIAG INF 3 + qqm=rrm + qqmt = rrmt + qqmB = rrmB + ! BC on rrmt + if (istret.eq.0) then + do is = 1, numscalar + qqmt(1 ,is) = zero + qqmt(ny,is) = beta_sc(is,2)*(-two/six/dy) + enddo + else + do is = 1, numscalar + qqmt(1 ,is) = zero + qqmt(ny,is) = beta_sc(is,2)*ppy(ny)*(-two/six/dy) + enddo + endif + do is = 1, 3 + rrmB(1,is) = zero + rrmB(ny,is) = zero + enddo !!! NCL = 1, npaire=0, neumann BC, odd function - ! - ! DIAG - aam10(1 )=zero - aam10(ny )=zero - aam10(2 )=-two*asjy-three*bsjy-two*csjy - aam10(ny-1 )=aam10(2) - aam10(3:ny-2)=-two*(asjy+bsjy+csjy) - call init_implicit_coef(aam10, aam10t) - if (istret==0) then - aam10 = one - xcst*aam10 - do is = 1, numscalar - aam10t(:,is) = one - xcst_sc(is)*aam10t(:,is) - enddo - else - aam10 = one/pp2y - xcst*aam10 - do is = 1, numscalar - aam10t(:,is) = one/pp2y - xcst_sc(is)*aam10t(:,is) - enddo - endif - ! - !DIAG SUP 1 - bbm10(1 )=zero - bbm10(ny )=zero - bbm10(2 )=asjy-csjy - bbm10(ny-1 )=asjy - bbm10(3 )=asjy - bbm10(ny-2 )=asjy-csjy - bbm10(4:ny-3)=asjy - call init_implicit_coef(bbm10, bbm10t) - if (istret==0) then - bbm10 = alsajy - xcst*bbm10 - do is = 1, numscalar - bbm10t(:,is) = alsajy - xcst_sc(is)*bbm10t(:,is) - enddo - else - bbm10(2:ny-1) = alsajy/pp2y(3:ny) - xcst*bbm10(2:ny-1) - do is = 1, numscalar - bbm10t(2:ny-1,is) = alsajy/pp2y(3:ny) - xcst_sc(is)*bbm10t(2:ny-1,is) - enddo - endif - !BC on bbm10 - bbm10(1 )=zero - bbm10(ny)=zero - ! BC on bbm10t - bbm10t(1 ,:)=zero - bbm10t(ny,:)=zero - ! - !DIAG SUP 2 - ccm10(1 )=zero - ccm10(ny )=zero - ccm10(2 )=bsjy - ccm10(ny-1 )=zero - ccm10(3 )=bsjy - ccm10(ny-2 )=bsjy - ccm10(4:ny-3)=bsjy - call init_implicit_coef(ccm10, ccm10t) - ccm10 = -xcst*ccm10 - do is = 1, numscalar - ccm10t(:,is) = -xcst_sc(is)*ccm10t(:,is) - enddo - ! - !DIAG SUP 3 - rrm10(1 )=zero - rrm10(ny )=zero - rrm10(2 )=csjy - rrm10(ny-1 )=zero - rrm10(3 )=csjy - rrm10(ny-2 )=zero - rrm10(4:ny-3)=csjy - call init_implicit_coef(rrm10, rrm10t) - rrm10 = -xcst*rrm10 - do is = 1, numscalar - rrm10t(:,is) = -xcst_sc(is)*rrm10t(:,is) - enddo - ! - !DIAG INF 1 - ddm10(1 )=zero - ddm10(ny )=zero - ddm10(2 )=asjy - ddm10(ny-1 )=asjy-csjy - ddm10(3 )=asjy-csjy - ddm10(ny-2 )=asjy - ddm10(4:ny-3)=asjy - call init_implicit_coef(ddm10, ddm10t) - if (istret==0) then - ddm10 = alsajy - xcst*ddm10 - do is = 1, numscalar - ddm10t(:,is) = alsajy - xcst_sc(is)*ddm10t(:,is) - enddo - else - ddm10(2:ny-1) = alsajy/pp2y(1:ny-2) - xcst*ddm10(2:ny-1) - do is = 1, numscalar - ddm10t(2:ny-1,is) = alsajy/pp2y(1:ny-2) - xcst_sc(is)*ddm10t(2:ny-1,is) - enddo - endif - !BC on ddm10 - ddm10(1 )=zero - ddm10(ny)=zero - ! BC on ddm10t - ddm10t(1 ,:)=zero - ddm10t(ny,:)=zero - ! - !DIAG INF 2 - eem10(1 )=zero - eem10(ny )=zero - eem10(2 )=zero - eem10(ny-1 )=bsjy - eem10(3 )=bsjy - eem10(ny-2 )=bsjy - eem10(4:ny-3)=bsjy - call init_implicit_coef(eem10, eem10t) - eem10 = -xcst*eem10 - do is = 1, numscalar - eem10t(:,is) = -xcst_sc(is)*eem10t(:,is) - enddo - ! - !DIAG INF 3 - qqm10(1 )=zero - qqm10(ny )=zero - qqm10(2 )=zero - qqm10(ny-1 )=csjy - qqm10(3 )=zero - qqm10(ny-2 )=csjy - qqm10(4:ny-3)=csjy - call init_implicit_coef(qqm10, qqm10t) - qqm10 = -xcst*qqm10 - do is = 1, numscalar - qqm10t(:,is) = -xcst_sc(is)*qqm10t(:,is) - enddo + ! + ! DIAG + aam10(1 )=zero + aam10(ny )=zero + aam10(2 )=-two*asjy-three*bsjy-two*csjy + aam10(ny-1 )=aam10(2) + aam10(3:ny-2)=-two*(asjy+bsjy+csjy) + call init_implicit_coef(aam10, aam10t) + if (istret==0) then + aam10 = one - xcst*aam10 + do is = 1, numscalar + aam10t(:,is) = one - xcst_sc(is)*aam10t(:,is) + enddo + else + aam10 = one/pp2y - xcst*aam10 + do is = 1, numscalar + aam10t(:,is) = one/pp2y - xcst_sc(is)*aam10t(:,is) + enddo + endif + ! + !DIAG SUP 1 + bbm10(1 )=zero + bbm10(ny )=zero + bbm10(2 )=asjy-csjy + bbm10(ny-1 )=asjy + bbm10(3 )=asjy + bbm10(ny-2 )=asjy-csjy + bbm10(4:ny-3)=asjy + call init_implicit_coef(bbm10, bbm10t) + if (istret==0) then + bbm10 = alsajy - xcst*bbm10 + do is = 1, numscalar + bbm10t(:,is) = alsajy - xcst_sc(is)*bbm10t(:,is) + enddo + else + bbm10(2:ny-1) = alsajy/pp2y(3:ny) - xcst*bbm10(2:ny-1) + do is = 1, numscalar + bbm10t(2:ny-1,is) = alsajy/pp2y(3:ny) - xcst_sc(is)*bbm10t(2:ny-1,is) + enddo + endif + !BC on bbm10 + bbm10(1 )=zero + bbm10(ny)=zero + ! BC on bbm10t + bbm10t(1 ,:)=zero + bbm10t(ny,:)=zero + ! + !DIAG SUP 2 + ccm10(1 )=zero + ccm10(ny )=zero + ccm10(2 )=bsjy + ccm10(ny-1 )=zero + ccm10(3 )=bsjy + ccm10(ny-2 )=bsjy + ccm10(4:ny-3)=bsjy + call init_implicit_coef(ccm10, ccm10t) + ccm10 = -xcst*ccm10 + do is = 1, numscalar + ccm10t(:,is) = -xcst_sc(is)*ccm10t(:,is) + enddo + ! + !DIAG SUP 3 + rrm10(1 )=zero + rrm10(ny )=zero + rrm10(2 )=csjy + rrm10(ny-1 )=zero + rrm10(3 )=csjy + rrm10(ny-2 )=zero + rrm10(4:ny-3)=csjy + call init_implicit_coef(rrm10, rrm10t) + rrm10 = -xcst*rrm10 + do is = 1, numscalar + rrm10t(:,is) = -xcst_sc(is)*rrm10t(:,is) + enddo + ! + !DIAG INF 1 + ddm10(1 )=zero + ddm10(ny )=zero + ddm10(2 )=asjy + ddm10(ny-1 )=asjy-csjy + ddm10(3 )=asjy-csjy + ddm10(ny-2 )=asjy + ddm10(4:ny-3)=asjy + call init_implicit_coef(ddm10, ddm10t) + if (istret==0) then + ddm10 = alsajy - xcst*ddm10 + do is = 1, numscalar + ddm10t(:,is) = alsajy - xcst_sc(is)*ddm10t(:,is) + enddo + else + ddm10(2:ny-1) = alsajy/pp2y(1:ny-2) - xcst*ddm10(2:ny-1) + do is = 1, numscalar + ddm10t(2:ny-1,is) = alsajy/pp2y(1:ny-2) - xcst_sc(is)*ddm10t(2:ny-1,is) + enddo + endif + !BC on ddm10 + ddm10(1 )=zero + ddm10(ny)=zero + ! BC on ddm10t + ddm10t(1 ,:)=zero + ddm10t(ny,:)=zero + ! + !DIAG INF 2 + eem10(1 )=zero + eem10(ny )=zero + eem10(2 )=zero + eem10(ny-1 )=bsjy + eem10(3 )=bsjy + eem10(ny-2 )=bsjy + eem10(4:ny-3)=bsjy + call init_implicit_coef(eem10, eem10t) + eem10 = -xcst*eem10 + do is = 1, numscalar + eem10t(:,is) = -xcst_sc(is)*eem10t(:,is) + enddo + ! + !DIAG INF 3 + qqm10(1 )=zero + qqm10(ny )=zero + qqm10(2 )=zero + qqm10(ny-1 )=csjy + qqm10(3 )=zero + qqm10(ny-2 )=csjy + qqm10(4:ny-3)=csjy + call init_implicit_coef(qqm10, qqm10t) + qqm10 = -xcst*qqm10 + do is = 1, numscalar + qqm10t(:,is) = -xcst_sc(is)*qqm10t(:,is) + enddo !!! NCL = 1, npaire=1, neumann BC, even function - ! - ! DIAG - aam11(1 )=-two*(asjy+bsjy+csjy) - aam11(ny )=aam11(1) - aam11(2 )=-two*asjy-bsjy-two*csjy - aam11(ny-1 )=aam11(2) - aam11(3:ny-2)=-two*(asjy+bsjy+csjy) - call init_implicit_coef(aam11, aam11t) - if (istret==0) then - aam11 = one - xcst*aam11 - do is = 1, numscalar - aam11t(:,is) = one - xcst_sc(is)*aam11t(:,is) - enddo - else - aam11 = one/pp2y - xcst*aam11 - do is = 1, numscalar - aam11t(:,is) = one/pp2y - xcst_sc(is)*aam11t(:,is) - enddo - endif - ! - !DIAG SUP 1 - bbm11(1 )=two*asjy - bbm11(ny )=zero - bbm11(2 )=asjy+csjy - bbm11(ny-1 )=asjy - bbm11(3 )=asjy - bbm11(ny-2 )=asjy+csjy - bbm11(4:ny-3)=asjy - call init_implicit_coef(bbm11, bbm11t) - if (istret==0) then - bbm11 = alsajy - xcst*bbm11 - do is = 1, numscalar - bbm11t(:,is) = alsajy - xcst_sc(is)*bbm11t(:,is) - enddo - else - bbm11(1:ny-1) = alsajy/pp2y(2:ny) - xcst*bbm11(1:ny-1) - do is = 1, numscalar - bbm11t(1:ny-1,is) = alsajy/pp2y(2:ny) - xcst_sc(is)*bbm11t(1:ny-1,is) - enddo - endif - !BC on bbm11 and bbm11t - if (istret==0) then - bbm11(1 )=bbm11(1)+alsajy - bbm11t(1 ,:)=bbm11t(1,:)+alsajy - else - bbm11(1 )=bbm11(1)+alsajy/pp2y(2) - bbm11t(1 ,:)=bbm11t(1,:)+alsajy/pp2y(2) - endif - bbm11(ny)=zero - bbm11t(ny,:)=zero - ! - !DIAG SUP 2 - ccm11(1 )=two*bsjy - ccm11(ny )=zero - ccm11(2 )=bsjy - ccm11(ny-1 )=zero - ccm11(3 )=bsjy - ccm11(ny-2 )=bsjy - ccm11(4:ny-3)=bsjy - call init_implicit_coef(ccm11, ccm11t) - ccm11 = -xcst*ccm11 - do is = 1, numscalar - ccm11t(:,is) = -xcst_sc(is)*ccm11t(:,is) - enddo - ! - !DIAG SUP 3 - rrm11(1 )=two*csjy - rrm11(ny )=zero - rrm11(2 )=csjy - rrm11(ny-1 )=zero - rrm11(3 )=csjy - rrm11(ny-2 )=zero - rrm11(4:ny-3)=csjy - call init_implicit_coef(rrm11, rrm11t) - rrm11 = -xcst*rrm11 - do is = 1, numscalar - rrm11t(:,is) = -xcst_sc(is)*rrm11t(:,is) - enddo - ! - !DIAG INF 1 - ddm11(1 )=zero - ddm11(ny )=two*asjy - ddm11(2 )=asjy - ddm11(ny-1 )=asjy+csjy - ddm11(3 )=asjy+csjy - ddm11(ny-2 )=asjy - ddm11(4:ny-3)=asjy - call init_implicit_coef(ddm11, ddm11t) - if (istret==0) then - ddm11 = alsajy - xcst*ddm11 - do is = 1, numscalar - ddm11t(:,is) = alsajy - xcst_sc(is)*ddm11t(:,is) - enddo - else - ddm11(2:ny) = alsajy/pp2y(1:ny-1) - xcst*ddm11(2:ny) - do is = 1, numscalar - ddm11t(2:ny,is) = alsajy/pp2y(1:ny-1) - xcst_sc(is)*ddm11t(2:ny,is) - enddo - endif - !BC on ddm11 and ddm11t - ddm11(1 )=zero - ddm11t(1 ,:)=zero - if (istret==0) then - ddm11(ny)=ddm11(ny)+alsajy!a1 - ddm11t(ny,:)=ddm11t(ny,:)+alsajy!a1 - else - ddm11(ny)=ddm11(ny)+alsajy/pp2y(ny-1)!a1 - ddm11t(ny,:)=ddm11t(ny,:)+alsajy/pp2y(ny-1)!a1 - endif - ! - !DIAG INF 2 - eem11(1 )=zero - eem11(ny )=two*bsjy - eem11(2 )=zero - eem11(ny-1 )=bsjy - eem11(3 )=bsjy - eem11(ny-2 )=bsjy - eem11(4:ny-3)=bsjy - call init_implicit_coef(eem11, eem11t) - eem11 = -xcst*eem11 - do is = 1, numscalar - eem11t(:,is) = -xcst_sc(is)*eem11t(:,is) - enddo - ! - !DIAG INF 3 - qqm11(1 )=zero - qqm11(ny )=two*csjy - qqm11(2 )=zero - qqm11(ny-1 )=csjy - qqm11(3 )=zero - qqm11(ny-2 )=csjy - qqm11(4:ny-3)=csjy - call init_implicit_coef(qqm11, qqm11t) - qqm11 = -xcst*qqm11 - do is = 1, numscalar - qqm11t(:,is) = -xcst_sc(is)*qqm11t(:,is) - enddo + ! + ! DIAG + aam11(1 )=-two*(asjy+bsjy+csjy) + aam11(ny )=aam11(1) + aam11(2 )=-two*asjy-bsjy-two*csjy + aam11(ny-1 )=aam11(2) + aam11(3:ny-2)=-two*(asjy+bsjy+csjy) + call init_implicit_coef(aam11, aam11t) + if (istret==0) then + aam11 = one - xcst*aam11 + do is = 1, numscalar + aam11t(:,is) = one - xcst_sc(is)*aam11t(:,is) + enddo + else + aam11 = one/pp2y - xcst*aam11 + do is = 1, numscalar + aam11t(:,is) = one/pp2y - xcst_sc(is)*aam11t(:,is) + enddo + endif + ! + !DIAG SUP 1 + bbm11(1 )=two*asjy + bbm11(ny )=zero + bbm11(2 )=asjy+csjy + bbm11(ny-1 )=asjy + bbm11(3 )=asjy + bbm11(ny-2 )=asjy+csjy + bbm11(4:ny-3)=asjy + call init_implicit_coef(bbm11, bbm11t) + if (istret==0) then + bbm11 = alsajy - xcst*bbm11 + do is = 1, numscalar + bbm11t(:,is) = alsajy - xcst_sc(is)*bbm11t(:,is) + enddo + else + bbm11(1:ny-1) = alsajy/pp2y(2:ny) - xcst*bbm11(1:ny-1) + do is = 1, numscalar + bbm11t(1:ny-1,is) = alsajy/pp2y(2:ny) - xcst_sc(is)*bbm11t(1:ny-1,is) + enddo + endif + !BC on bbm11 and bbm11t + if (istret==0) then + bbm11(1 )=bbm11(1)+alsajy + bbm11t(1 ,:)=bbm11t(1,:)+alsajy + else + bbm11(1 )=bbm11(1)+alsajy/pp2y(2) + bbm11t(1 ,:)=bbm11t(1,:)+alsajy/pp2y(2) + endif + bbm11(ny)=zero + bbm11t(ny,:)=zero + ! + !DIAG SUP 2 + ccm11(1 )=two*bsjy + ccm11(ny )=zero + ccm11(2 )=bsjy + ccm11(ny-1 )=zero + ccm11(3 )=bsjy + ccm11(ny-2 )=bsjy + ccm11(4:ny-3)=bsjy + call init_implicit_coef(ccm11, ccm11t) + ccm11 = -xcst*ccm11 + do is = 1, numscalar + ccm11t(:,is) = -xcst_sc(is)*ccm11t(:,is) + enddo + ! + !DIAG SUP 3 + rrm11(1 )=two*csjy + rrm11(ny )=zero + rrm11(2 )=csjy + rrm11(ny-1 )=zero + rrm11(3 )=csjy + rrm11(ny-2 )=zero + rrm11(4:ny-3)=csjy + call init_implicit_coef(rrm11, rrm11t) + rrm11 = -xcst*rrm11 + do is = 1, numscalar + rrm11t(:,is) = -xcst_sc(is)*rrm11t(:,is) + enddo + ! + !DIAG INF 1 + ddm11(1 )=zero + ddm11(ny )=two*asjy + ddm11(2 )=asjy + ddm11(ny-1 )=asjy+csjy + ddm11(3 )=asjy+csjy + ddm11(ny-2 )=asjy + ddm11(4:ny-3)=asjy + call init_implicit_coef(ddm11, ddm11t) + if (istret==0) then + ddm11 = alsajy - xcst*ddm11 + do is = 1, numscalar + ddm11t(:,is) = alsajy - xcst_sc(is)*ddm11t(:,is) + enddo + else + ddm11(2:ny) = alsajy/pp2y(1:ny-1) - xcst*ddm11(2:ny) + do is = 1, numscalar + ddm11t(2:ny,is) = alsajy/pp2y(1:ny-1) - xcst_sc(is)*ddm11t(2:ny,is) + enddo + endif + !BC on ddm11 and ddm11t + ddm11(1 )=zero + ddm11t(1 ,:)=zero + if (istret==0) then + ddm11(ny)=ddm11(ny)+alsajy!a1 + ddm11t(ny,:)=ddm11t(ny,:)+alsajy!a1 + else + ddm11(ny)=ddm11(ny)+alsajy/pp2y(ny-1)!a1 + ddm11t(ny,:)=ddm11t(ny,:)+alsajy/pp2y(ny-1)!a1 + endif + ! + !DIAG INF 2 + eem11(1 )=zero + eem11(ny )=two*bsjy + eem11(2 )=zero + eem11(ny-1 )=bsjy + eem11(3 )=bsjy + eem11(ny-2 )=bsjy + eem11(4:ny-3)=bsjy + call init_implicit_coef(eem11, eem11t) + eem11 = -xcst*eem11 + do is = 1, numscalar + eem11t(:,is) = -xcst_sc(is)*eem11t(:,is) + enddo + ! + !DIAG INF 3 + qqm11(1 )=zero + qqm11(ny )=two*csjy + qqm11(2 )=zero + qqm11(ny-1 )=csjy + qqm11(3 )=zero + qqm11(ny-2 )=csjy + qqm11(4:ny-3)=csjy + call init_implicit_coef(qqm11, qqm11t) + qqm11 = -xcst*qqm11 + do is = 1, numscalar + qqm11t(:,is) = -xcst_sc(is)*qqm11t(:,is) + enddo !!! NXL = 0 - !DIAG - if (istret==0) then - aam0 = one-xcst*(-two*(asjy+bsjy+csjy)) - do is = 1, numscalar - aam0t(:,is) = one-xcst_sc(is)*(-two*(asjy+bsjy+csjy)) - enddo - else - aam0 = one/pp2y-xcst*(-two*(asjy+bsjy+csjy)) - do is = 1, numscalar - aam0t(:,is) = one/pp2y-xcst_sc(is)*(-two*(asjy+bsjy+csjy)) - enddo - endif - ! - !DIAG SUP 1 - if (istret==0) then - bbm0 = alsajy-xcst*asjy - do is = 1, numscalar - bbm0t(:,is) = alsajy-xcst_sc(is)*asjy - enddo - else - bbm0(1:ny-1) = alsajy/pp2y(2:ny) -xcst*asjy - bbm0(ny) = alsajy/pp2y(1) -xcst*asjy - do is = 1, numscalar - bbm0t(1:ny-1,is) = alsajy/pp2y(2:ny) -xcst_sc(is)*asjy - bbm0t(ny,is) = alsajy/pp2y(1) -xcst_sc(is)*asjy - enddo - endif - ! - !DIAG SUP 2 - ccm0 = -xcst*bsjy - do is = 1, numscalar - ccm0t(:,is) = -xcst_sc(is)*bsjy - enddo - ! - !DIAG SUP 3 - rrm0 = -xcst*csjy - do is = 1, numscalar - rrm0t(:,is) = -xcst_sc(is)*csjy - enddo - ! - !DIAG INF 1 - if (istret==0) then - ddm0=bbm0 - ddm0t = bbm0t - else - ddm0(1) = alsajy/pp2y(ny) -xcst*asjy - ddm0(2:ny) = alsajy/pp2y(1:ny-1) -xcst*asjy - do is = 1, numscalar - ddm0t(1,is) = alsajy/pp2y(ny) -xcst_sc(is)*asjy - ddm0t(2:ny,is) = alsajy/pp2y(1:ny-1) -xcst_sc(is)*asjy - enddo - endif - ! - !DIAG INF 2 - eem0=ccm0 - eem0t=ccm0t - ! - !DIAG INF 3 - qqm0=rrm0 - qqm0t=rrm0t - - ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 - aam120=aam10; aam120(ny-2:ny)=aam(ny-2:ny) - bbm120=bbm10; bbm120(ny-2:ny)=bbm(ny-2:ny) - ccm120=ccm10; ccm120(ny-2:ny)=ccm(ny-2:ny) - ddm120=ddm10; ddm120(ny-2:ny)=ddm(ny-2:ny) - eem120=eem10; eem120(ny-2:ny)=eem(ny-2:ny) - qqm120=qqm10; qqm120(ny-2:ny)=qqm(ny-2:ny) - rrm120=rrm10; rrm120(ny-2:ny)=rrm(ny-2:ny) - ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 - aam121=aam11; aam121(ny-2:ny)=aam(ny-2:ny) - bbm121=bbm11; bbm121(ny-2:ny)=bbm(ny-2:ny) - ccm121=ccm11; ccm121(ny-2:ny)=ccm(ny-2:ny) - ddm121=ddm11; ddm121(ny-2:ny)=ddm(ny-2:ny) - eem121=eem11; eem121(ny-2:ny)=eem(ny-2:ny) - qqm121=qqm11; qqm121(ny-2:ny)=qqm(ny-2:ny) - rrm121=rrm11; rrm121(ny-2:ny)=rrm(ny-2:ny) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 - aam210=aam; aam210(ny-2:ny)=aam10(ny-2:ny) - bbm210=bbm; bbm210(ny-2:ny)=bbm10(ny-2:ny) - ccm210=ccm; ccm210(ny-2:ny)=vvm10(ny-2:ny) - ddm210=ddm; ddm210(ny-2:ny)=ddm10(ny-2:ny) - eem210=eem; eem210(ny-2:ny)=eem10(ny-2:ny) - qqm210=qqm; qqm210(ny-2:ny)=qqm10(ny-2:ny) - rrm210=rrm; rrm210(ny-2:ny)=rrm10(ny-2:ny) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 - aam211=aam; aam211(ny-2:ny)=aam11(ny-2:ny) - bbm211=bbm; bbm211(ny-2:ny)=bbm11(ny-2:ny) - ccm211=ccm; ccm211(ny-2:ny)=ccm11(ny-2:ny) - ddm211=ddm; ddm211(ny-2:ny)=ddm11(ny-2:ny) - eem211=eem; eem211(ny-2:ny)=eem11(ny-2:ny) - qqm211=qqm; qqm211(ny-2:ny)=qqm11(ny-2:ny) - rrm211=rrm; rrm211(ny-2:ny)=rrm11(ny-2:ny) - - ! scalars, ncly1 = 1, nclyn = 2, npaire = 0 - aam120t=aam10t; aam120t(ny-2:ny,:)=aamt(ny-2:ny,:) - bbm120t=bbm10t; bbm120t(ny-2:ny,:)=bbmt(ny-2:ny,:) - ccm120t=ccm10t; ccm120t(ny-2:ny,:)=ccmt(ny-2:ny,:) - ddm120t=ddm10t; ddm120t(ny-2:ny,:)=ddmt(ny-2:ny,:) - eem120t=eem10t; eem120t(ny-2:ny,:)=eemt(ny-2:ny,:) - qqm120t=qqm10t; qqm120t(ny-2:ny,:)=qqmt(ny-2:ny,:) - rrm120t=rrm10t; rrm120t(ny-2:ny,:)=rrmt(ny-2:ny,:) - ! scalars, ncly1 = 1, nclyn = 2, npaire = 1 - aam121t=aam11t; aam121t(ny-2:ny,:)=aamt(ny-2:ny,:) - bbm121t=bbm11t; bbm121t(ny-2:ny,:)=bbmt(ny-2:ny,:) - ccm121t=ccm11t; ccm121t(ny-2:ny,:)=ccmt(ny-2:ny,:) - ddm121t=ddm11t; ddm121t(ny-2:ny,:)=ddmt(ny-2:ny,:) - eem121t=eem11t; eem121t(ny-2:ny,:)=eemt(ny-2:ny,:) - qqm121t=qqm11t; qqm121t(ny-2:ny,:)=qqmt(ny-2:ny,:) - rrm121t=rrm11t; rrm121t(ny-2:ny,:)=rrmt(ny-2:ny,:) - ! scalars, ncly1 = 2, nclyn = 1, npaire = 0 - aam210t=aamt; aam210t(ny-2:ny,:)=aam10t(ny-2:ny,:) - bbm210t=bbmt; bbm210t(ny-2:ny,:)=bbm10t(ny-2:ny,:) - ccm210t=ccmt; ccm210t(ny-2:ny,:)=vvm10t(ny-2:ny,:) - ddm210t=ddmt; ddm210t(ny-2:ny,:)=ddm10t(ny-2:ny,:) - eem210t=eemt; eem210t(ny-2:ny,:)=eem10t(ny-2:ny,:) - qqm210t=qqmt; qqm210t(ny-2:ny,:)=qqm10t(ny-2:ny,:) - rrm210t=rrmt; rrm210t(ny-2:ny,:)=rrm10t(ny-2:ny,:) - ! scalars, ncly1 = 2, nclyn = 1, npaire = 1 - aam211t=aamt; aam211t(ny-2:ny,:)=aam11t(ny-2:ny,:) - bbm211t=bbmt; bbm211t(ny-2:ny,:)=bbm11t(ny-2:ny,:) - ccm211t=ccmt; ccm211t(ny-2:ny,:)=ccm11t(ny-2:ny,:) - ddm211t=ddmt; ddm211t(ny-2:ny,:)=ddm11t(ny-2:ny,:) - eem211t=eemt; eem211t(ny-2:ny,:)=eem11t(ny-2:ny,:) - qqm211t=qqmt; qqm211t(ny-2:ny,:)=qqm11t(ny-2:ny,:) - rrm211t=rrmt; rrm211t(ny-2:ny,:)=rrm11t(ny-2:ny,:) - - else + !DIAG + if (istret==0) then + aam0 = one-xcst*(-two*(asjy+bsjy+csjy)) + do is = 1, numscalar + aam0t(:,is) = one-xcst_sc(is)*(-two*(asjy+bsjy+csjy)) + enddo + else + aam0 = one/pp2y-xcst*(-two*(asjy+bsjy+csjy)) + do is = 1, numscalar + aam0t(:,is) = one/pp2y-xcst_sc(is)*(-two*(asjy+bsjy+csjy)) + enddo + endif + ! + !DIAG SUP 1 + if (istret==0) then + bbm0 = alsajy-xcst*asjy + do is = 1, numscalar + bbm0t(:,is) = alsajy-xcst_sc(is)*asjy + enddo + else + bbm0(1:ny-1) = alsajy/pp2y(2:ny) -xcst*asjy + bbm0(ny) = alsajy/pp2y(1) -xcst*asjy + do is = 1, numscalar + bbm0t(1:ny-1,is) = alsajy/pp2y(2:ny) -xcst_sc(is)*asjy + bbm0t(ny,is) = alsajy/pp2y(1) -xcst_sc(is)*asjy + enddo + endif + ! + !DIAG SUP 2 + ccm0 = -xcst*bsjy + do is = 1, numscalar + ccm0t(:,is) = -xcst_sc(is)*bsjy + enddo + ! + !DIAG SUP 3 + rrm0 = -xcst*csjy + do is = 1, numscalar + rrm0t(:,is) = -xcst_sc(is)*csjy + enddo + ! + !DIAG INF 1 + if (istret==0) then + ddm0=bbm0 + ddm0t = bbm0t + else + ddm0(1) = alsajy/pp2y(ny) -xcst*asjy + ddm0(2:ny) = alsajy/pp2y(1:ny-1) -xcst*asjy + do is = 1, numscalar + ddm0t(1,is) = alsajy/pp2y(ny) -xcst_sc(is)*asjy + ddm0t(2:ny,is) = alsajy/pp2y(1:ny-1) -xcst_sc(is)*asjy + enddo + endif + ! + !DIAG INF 2 + eem0=ccm0 + eem0t=ccm0t + ! + !DIAG INF 3 + qqm0=rrm0 + qqm0t=rrm0t + + ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 + aam120=aam10; aam120(ny-2:ny)=aam(ny-2:ny) + bbm120=bbm10; bbm120(ny-2:ny)=bbm(ny-2:ny) + ccm120=ccm10; ccm120(ny-2:ny)=ccm(ny-2:ny) + ddm120=ddm10; ddm120(ny-2:ny)=ddm(ny-2:ny) + eem120=eem10; eem120(ny-2:ny)=eem(ny-2:ny) + qqm120=qqm10; qqm120(ny-2:ny)=qqm(ny-2:ny) + rrm120=rrm10; rrm120(ny-2:ny)=rrm(ny-2:ny) + ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 + aam121=aam11; aam121(ny-2:ny)=aam(ny-2:ny) + bbm121=bbm11; bbm121(ny-2:ny)=bbm(ny-2:ny) + ccm121=ccm11; ccm121(ny-2:ny)=ccm(ny-2:ny) + ddm121=ddm11; ddm121(ny-2:ny)=ddm(ny-2:ny) + eem121=eem11; eem121(ny-2:ny)=eem(ny-2:ny) + qqm121=qqm11; qqm121(ny-2:ny)=qqm(ny-2:ny) + rrm121=rrm11; rrm121(ny-2:ny)=rrm(ny-2:ny) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 + aam210=aam; aam210(ny-2:ny)=aam10(ny-2:ny) + bbm210=bbm; bbm210(ny-2:ny)=bbm10(ny-2:ny) + ccm210=ccm; ccm210(ny-2:ny)=vvm10(ny-2:ny) + ddm210=ddm; ddm210(ny-2:ny)=ddm10(ny-2:ny) + eem210=eem; eem210(ny-2:ny)=eem10(ny-2:ny) + qqm210=qqm; qqm210(ny-2:ny)=qqm10(ny-2:ny) + rrm210=rrm; rrm210(ny-2:ny)=rrm10(ny-2:ny) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 + aam211=aam; aam211(ny-2:ny)=aam11(ny-2:ny) + bbm211=bbm; bbm211(ny-2:ny)=bbm11(ny-2:ny) + ccm211=ccm; ccm211(ny-2:ny)=ccm11(ny-2:ny) + ddm211=ddm; ddm211(ny-2:ny)=ddm11(ny-2:ny) + eem211=eem; eem211(ny-2:ny)=eem11(ny-2:ny) + qqm211=qqm; qqm211(ny-2:ny)=qqm11(ny-2:ny) + rrm211=rrm; rrm211(ny-2:ny)=rrm11(ny-2:ny) + + ! scalars, ncly1 = 1, nclyn = 2, npaire = 0 + aam120t=aam10t; aam120t(ny-2:ny,:)=aamt(ny-2:ny,:) + bbm120t=bbm10t; bbm120t(ny-2:ny,:)=bbmt(ny-2:ny,:) + ccm120t=ccm10t; ccm120t(ny-2:ny,:)=ccmt(ny-2:ny,:) + ddm120t=ddm10t; ddm120t(ny-2:ny,:)=ddmt(ny-2:ny,:) + eem120t=eem10t; eem120t(ny-2:ny,:)=eemt(ny-2:ny,:) + qqm120t=qqm10t; qqm120t(ny-2:ny,:)=qqmt(ny-2:ny,:) + rrm120t=rrm10t; rrm120t(ny-2:ny,:)=rrmt(ny-2:ny,:) + ! scalars, ncly1 = 1, nclyn = 2, npaire = 1 + aam121t=aam11t; aam121t(ny-2:ny,:)=aamt(ny-2:ny,:) + bbm121t=bbm11t; bbm121t(ny-2:ny,:)=bbmt(ny-2:ny,:) + ccm121t=ccm11t; ccm121t(ny-2:ny,:)=ccmt(ny-2:ny,:) + ddm121t=ddm11t; ddm121t(ny-2:ny,:)=ddmt(ny-2:ny,:) + eem121t=eem11t; eem121t(ny-2:ny,:)=eemt(ny-2:ny,:) + qqm121t=qqm11t; qqm121t(ny-2:ny,:)=qqmt(ny-2:ny,:) + rrm121t=rrm11t; rrm121t(ny-2:ny,:)=rrmt(ny-2:ny,:) + ! scalars, ncly1 = 2, nclyn = 1, npaire = 0 + aam210t=aamt; aam210t(ny-2:ny,:)=aam10t(ny-2:ny,:) + bbm210t=bbmt; bbm210t(ny-2:ny,:)=bbm10t(ny-2:ny,:) + ccm210t=ccmt; ccm210t(ny-2:ny,:)=vvm10t(ny-2:ny,:) + ddm210t=ddmt; ddm210t(ny-2:ny,:)=ddm10t(ny-2:ny,:) + eem210t=eemt; eem210t(ny-2:ny,:)=eem10t(ny-2:ny,:) + qqm210t=qqmt; qqm210t(ny-2:ny,:)=qqm10t(ny-2:ny,:) + rrm210t=rrmt; rrm210t(ny-2:ny,:)=rrm10t(ny-2:ny,:) + ! scalars, ncly1 = 2, nclyn = 1, npaire = 1 + aam211t=aamt; aam211t(ny-2:ny,:)=aam11t(ny-2:ny,:) + bbm211t=bbmt; bbm211t(ny-2:ny,:)=bbm11t(ny-2:ny,:) + ccm211t=ccmt; ccm211t(ny-2:ny,:)=ccm11t(ny-2:ny,:) + ddm211t=ddmt; ddm211t(ny-2:ny,:)=ddm11t(ny-2:ny,:) + eem211t=eemt; eem211t(ny-2:ny,:)=eem11t(ny-2:ny,:) + qqm211t=qqmt; qqm211t(ny-2:ny,:)=qqm11t(ny-2:ny,:) + rrm211t=rrmt; rrm211t(ny-2:ny,:)=rrm11t(ny-2:ny,:) + + else !!!!!!!!!!!!!!!!!!!!!! - !9-DIAG !! + !9-DIAG !! !!!!!!!!!!!!!!!!!!!!!! !!! NCL = 2, dirichlet BC - ! - !DIAG - aam(1 )=as1y - aam(ny )=asny - aam(2 )=-two*as2y - aam(ny-1 )=-two*asmy - aam(3 )=-two*(as3y+bs3y) - aam(ny-2 )=-two*(asty+bsty) - aam(4 )=-two*(as4y+bs4y+cs4y) - aam(ny-3 )=-two*(astty+bstty+cstty) - aam(5:ny-4)=-two*(asjy+bsjy+csjy+dsjy) - if (istret==0) then - aam = one-xcst*aam - else - aam = one/pp2y -xcst*aam - endif - !BC on aam - aam(1 )=one - aam(ny)=one - ! - !DIAG SUP 1 - bbm(1 )=bs1y - bbm(ny )=bsny - bbm(2 )=as2y - bbm(ny-1 )=asmy - bbm(3 )=as3y - bbm(ny-2 )=asty - bbm(4 )=as4y - bbm(ny-3 )=astty - bbm(5:ny-4)=asjy - bbm = -xcst*bbm - if (istret==0) then - bbm(2 )=bbm(2 )+alsa2y - bbm(ny-1 )=bbm(ny-1 )+alsamy - bbm(3 )=bbm(3 )+alsa3y - bbm(ny-2 )=bbm(ny-2 )+alsaty - bbm(4 )=bbm(4 )+alsa4y - bbm(ny-3 )=bbm(ny-3 )+alsatty - bbm(5:ny-4)=bbm(5:ny-4)+alsajy - else - bbm(2 )=bbm(2 )+alsa2y/pp2y(3) - bbm(ny-1 )=bbm(ny-1 )+alsamy/pp2y(ny) - bbm(3 )=bbm(3 )+alsa3y/pp2y(4) - bbm(ny-2 )=bbm(ny-2 )+alsaty/pp2y(ny-1) - bbm(4 )=bbm(4 )+alsa4y/pp2y(5) - bbm(ny-3 )=bbm(ny-3 )+alsatty/pp2y(ny-2) - bbm(5:ny-4)=bbm(5:ny-4)+alsajy/pp2y(6:ny-3) - endif - !BC on bbm - bbm(1 )=zero - bbm(ny)=zero - ! - !DIAG SUP 2 - ccm(1 )=cs1y - ccm(ny )=csny - ccm(2 )=zero!bs2y - ccm(ny-1 )=zero!bsmy - ccm(3 )=bs3y - ccm(ny-2 )=bsty - ccm(4 )=bs4y - ccm(ny-3 )=bstty - ccm(5:ny-4)=bsjy - ccm = -xcst*ccm - !BC on ccm - ccm(1 )=zero - ccm(ny)=zero - ! - !DIAG SUP 3 - rrm(1 )=ds1y - rrm(ny )=dsny - rrm(2 )=zero!cs2y - rrm(ny-1 )=zero!csmy - rrm(3 )=zero!cs3y - rrm(ny-2 )=zero!csty - rrm(4 )=cs4y - rrm(ny-3 )=cstty - rrm(5:ny-4)=csjy - rrm = -xcst*rrm - !BC on rrm - rrm(1 )=zero - rrm(ny)=zero - ! - !DIAG SUP 4 - ttm(1 )=zero - ttm(ny )=zero - ttm(2 )=zero - ttm(ny-1 )=zero - ttm(3 )=zero!ds3y - ttm(ny-2 )=zero!dsty - ttm(4 )=zero!ds4y - ttm(ny-3 )=zero!dstty - ttm(5:ny-4)=dsjy - ttm = -xcst*ttm - !BC on ttm - ttm(1 )=zero - ttm(ny)=zero - ! - !DIAG INF 1 - if (istret==0) then - ddm=bbm - else - ddm(1 )=bs1y - ddm(ny )=bsny - ddm(2 )=as2y - ddm(ny-1 )=asmy - ddm(3 )=as3y - ddm(ny-2 )=asty - ddm(4 )=as4y - ddm(ny-3 )=astty - ddm(5:ny-4)=asjy - ddm = -xcst*ddm - ddm(2 )=ddm(2 )+alsa2y/pp2y(1) - ddm(ny-1 )=ddm(ny-1 )+alsamy/pp2y(ny-2) - ddm(3 )=ddm(3 )+alsa3y/pp2y(2) - ddm(ny-2 )=ddm(ny-2 )+alsaty/pp2y(ny-3) - ddm(4 )=ddm(4 )+alsa4y/pp2y(3) - ddm(ny-3 )=ddm(ny-3 )+alsatty/pp2y(ny-4) - ddm(5:ny-4)=ddm(5:ny-4)+alsajy/pp2y(4:ny-5) - !BC on ddm - ddm(1 )=zero - ddm(ny)=zero - endif - ! - !DIAG INF 2 - eem=ccm - ! - !DIAG INF 3 - qqm=rrm - - !DIAG INF 4 - uum=ttm + ! + !DIAG + aam(1 )=as1y + aam(ny )=asny + aam(2 )=-two*as2y + aam(ny-1 )=-two*asmy + aam(3 )=-two*(as3y+bs3y) + aam(ny-2 )=-two*(asty+bsty) + aam(4 )=-two*(as4y+bs4y+cs4y) + aam(ny-3 )=-two*(astty+bstty+cstty) + aam(5:ny-4)=-two*(asjy+bsjy+csjy+dsjy) + if (istret==0) then + aam = one-xcst*aam + else + aam = one/pp2y -xcst*aam + endif + !BC on aam + aam(1 )=one + aam(ny)=one + ! + !DIAG SUP 1 + bbm(1 )=bs1y + bbm(ny )=bsny + bbm(2 )=as2y + bbm(ny-1 )=asmy + bbm(3 )=as3y + bbm(ny-2 )=asty + bbm(4 )=as4y + bbm(ny-3 )=astty + bbm(5:ny-4)=asjy + bbm = -xcst*bbm + if (istret==0) then + bbm(2 )=bbm(2 )+alsa2y + bbm(ny-1 )=bbm(ny-1 )+alsamy + bbm(3 )=bbm(3 )+alsa3y + bbm(ny-2 )=bbm(ny-2 )+alsaty + bbm(4 )=bbm(4 )+alsa4y + bbm(ny-3 )=bbm(ny-3 )+alsatty + bbm(5:ny-4)=bbm(5:ny-4)+alsajy + else + bbm(2 )=bbm(2 )+alsa2y/pp2y(3) + bbm(ny-1 )=bbm(ny-1 )+alsamy/pp2y(ny) + bbm(3 )=bbm(3 )+alsa3y/pp2y(4) + bbm(ny-2 )=bbm(ny-2 )+alsaty/pp2y(ny-1) + bbm(4 )=bbm(4 )+alsa4y/pp2y(5) + bbm(ny-3 )=bbm(ny-3 )+alsatty/pp2y(ny-2) + bbm(5:ny-4)=bbm(5:ny-4)+alsajy/pp2y(6:ny-3) + endif + !BC on bbm + bbm(1 )=zero + bbm(ny)=zero + ! + !DIAG SUP 2 + ccm(1 )=cs1y + ccm(ny )=csny + ccm(2 )=zero!bs2y + ccm(ny-1 )=zero!bsmy + ccm(3 )=bs3y + ccm(ny-2 )=bsty + ccm(4 )=bs4y + ccm(ny-3 )=bstty + ccm(5:ny-4)=bsjy + ccm = -xcst*ccm + !BC on ccm + ccm(1 )=zero + ccm(ny)=zero + ! + !DIAG SUP 3 + rrm(1 )=ds1y + rrm(ny )=dsny + rrm(2 )=zero!cs2y + rrm(ny-1 )=zero!csmy + rrm(3 )=zero!cs3y + rrm(ny-2 )=zero!csty + rrm(4 )=cs4y + rrm(ny-3 )=cstty + rrm(5:ny-4)=csjy + rrm = -xcst*rrm + !BC on rrm + rrm(1 )=zero + rrm(ny)=zero + ! + !DIAG SUP 4 + ttm(1 )=zero + ttm(ny )=zero + ttm(2 )=zero + ttm(ny-1 )=zero + ttm(3 )=zero!ds3y + ttm(ny-2 )=zero!dsty + ttm(4 )=zero!ds4y + ttm(ny-3 )=zero!dstty + ttm(5:ny-4)=dsjy + ttm = -xcst*ttm + !BC on ttm + ttm(1 )=zero + ttm(ny)=zero + ! + !DIAG INF 1 + if (istret==0) then + ddm=bbm + else + ddm(1 )=bs1y + ddm(ny )=bsny + ddm(2 )=as2y + ddm(ny-1 )=asmy + ddm(3 )=as3y + ddm(ny-2 )=asty + ddm(4 )=as4y + ddm(ny-3 )=astty + ddm(5:ny-4)=asjy + ddm = -xcst*ddm + ddm(2 )=ddm(2 )+alsa2y/pp2y(1) + ddm(ny-1 )=ddm(ny-1 )+alsamy/pp2y(ny-2) + ddm(3 )=ddm(3 )+alsa3y/pp2y(2) + ddm(ny-2 )=ddm(ny-2 )+alsaty/pp2y(ny-3) + ddm(4 )=ddm(4 )+alsa4y/pp2y(3) + ddm(ny-3 )=ddm(ny-3 )+alsatty/pp2y(ny-4) + ddm(5:ny-4)=ddm(5:ny-4)+alsajy/pp2y(4:ny-5) + !BC on ddm + ddm(1 )=zero + ddm(ny)=zero + endif + ! + !DIAG INF 2 + eem=ccm + ! + !DIAG INF 3 + qqm=rrm + + !DIAG INF 4 + uum=ttm !!! NCL = 1, npaire=0, neumann BC, odd function - ! - ! DIAG - aam10(1 )=zero - aam10(ny )=zero - aam10(2 )=-two*asjy-three*bsjy-two*csjy - aam10(ny-1 )=aam10(2) - aam10(3:ny-2)=-two*(asjy+bsjy+csjy) - if (istret==0) then - aam10 = one - xcst*aam10 - else - aam10 = one/pp2y - xcst*aam10 - endif - ! - !DIAG SUP 1 - bbm10(1 )=zero - bbm10(ny )=zero - bbm10(2 )=asjy-csjy - bbm10(ny-1 )=asjy - bbm10(3 )=asjy - bbm10(ny-2 )=asjy-csjy - bbm10(4:ny-3)=asjy - if (istret==0) then - bbm10 = alsajy - xcst*bbm10 - else - bbm10(2:ny-1) = alsajy/pp2y(3:ny) - xcst*bbm10(2:ny-1) - endif - !Bc on bbm10 - bbm10(1 )=zero - bbm10(ny)=zero - ! - !DIAG SUP 2 - ccm10(1 )=zero - ccm10(ny )=zero - ccm10(2 )=bsjy - ccm10(ny-1 )=zero - ccm10(3 )=bsjy - ccm10(ny-2 )=bsjy - ccm10(4:ny-3)=bsjy - ccm10 = -xcst*ccm10 - ! - !DIAG SUP 3 - rrm10(1 )=zero - rrm10(ny )=zero - rrm10(2 )=csjy - rrm10(ny-1 )=zero - rrm10(3 )=csjy - rrm10(ny-2 )=zero - rrm10(4:ny-3)=csjy - rrm10 = -xcst*rrm10 - ! - !DIAG INF 1 - ddm10(1 )=zero - ddm10(ny )=zero - ddm10(2 )=asjy - ddm10(ny-1 )=asjy-csjy - ddm10(3 )=asjy-csjy - ddm10(ny-2 )=asjy - ddm10(4:ny-3)=asjy - if (istret==0) then - ddm10 = alsajy - xcst*ddm10 - else - ddm10(2:ny-1) = alsajy/pp2y(1:ny-2) - xcst*ddm10(2:ny-1) - endif - !BC on ddm10 - ddm10(1 )=zero - ddm10(ny)=zero - ! - !DIAG INF 2 - eem10(1 )=zero - eem10(ny )=zero - eem10(2 )=zero - eem10(ny-1 )=bsjy - eem10(3 )=bsjy - eem10(ny-2 )=bsjy - eem10(4:ny-3)=bsjy - eem10 = -xcst*eem10 - ! - !DIAG INF 3 - qqm10(1 )=zero - qqm10(ny )=zero - qqm10(2 )=zero - qqm10(ny-1 )=csjy - qqm10(3 )=zero - qqm10(ny-2 )=csjy - qqm10(4:ny-3)=csjy - qqm10 = -xcst*qqm10 + ! + ! DIAG + aam10(1 )=zero + aam10(ny )=zero + aam10(2 )=-two*asjy-three*bsjy-two*csjy + aam10(ny-1 )=aam10(2) + aam10(3:ny-2)=-two*(asjy+bsjy+csjy) + if (istret==0) then + aam10 = one - xcst*aam10 + else + aam10 = one/pp2y - xcst*aam10 + endif + ! + !DIAG SUP 1 + bbm10(1 )=zero + bbm10(ny )=zero + bbm10(2 )=asjy-csjy + bbm10(ny-1 )=asjy + bbm10(3 )=asjy + bbm10(ny-2 )=asjy-csjy + bbm10(4:ny-3)=asjy + if (istret==0) then + bbm10 = alsajy - xcst*bbm10 + else + bbm10(2:ny-1) = alsajy/pp2y(3:ny) - xcst*bbm10(2:ny-1) + endif + !Bc on bbm10 + bbm10(1 )=zero + bbm10(ny)=zero + ! + !DIAG SUP 2 + ccm10(1 )=zero + ccm10(ny )=zero + ccm10(2 )=bsjy + ccm10(ny-1 )=zero + ccm10(3 )=bsjy + ccm10(ny-2 )=bsjy + ccm10(4:ny-3)=bsjy + ccm10 = -xcst*ccm10 + ! + !DIAG SUP 3 + rrm10(1 )=zero + rrm10(ny )=zero + rrm10(2 )=csjy + rrm10(ny-1 )=zero + rrm10(3 )=csjy + rrm10(ny-2 )=zero + rrm10(4:ny-3)=csjy + rrm10 = -xcst*rrm10 + ! + !DIAG INF 1 + ddm10(1 )=zero + ddm10(ny )=zero + ddm10(2 )=asjy + ddm10(ny-1 )=asjy-csjy + ddm10(3 )=asjy-csjy + ddm10(ny-2 )=asjy + ddm10(4:ny-3)=asjy + if (istret==0) then + ddm10 = alsajy - xcst*ddm10 + else + ddm10(2:ny-1) = alsajy/pp2y(1:ny-2) - xcst*ddm10(2:ny-1) + endif + !BC on ddm10 + ddm10(1 )=zero + ddm10(ny)=zero + ! + !DIAG INF 2 + eem10(1 )=zero + eem10(ny )=zero + eem10(2 )=zero + eem10(ny-1 )=bsjy + eem10(3 )=bsjy + eem10(ny-2 )=bsjy + eem10(4:ny-3)=bsjy + eem10 = -xcst*eem10 + ! + !DIAG INF 3 + qqm10(1 )=zero + qqm10(ny )=zero + qqm10(2 )=zero + qqm10(ny-1 )=csjy + qqm10(3 )=zero + qqm10(ny-2 )=csjy + qqm10(4:ny-3)=csjy + qqm10 = -xcst*qqm10 !!! NCL = 1, npaire=1, neumann BC, even function - ! - ! DIAG - aam11(1 )=-two*(asjy+bsjy+csjy) - aam11(ny )=aam11(1) - aam11(2 )=-two*asjy-bsjy-two*csjy - aam11(ny-1 )=aam11(2) - aam11(3:ny-2)=-two*(asjy+bsjy+csjy) - if (istret==0) then - aam11 = one - xcst*aam11 - else - aam11 = one/pp2y - xcst*aam11 - endif - ! - !DIAG SUP 1 - bbm11(1 )=two*asjy - bbm11(ny )=zero - bbm11(2 )=asjy+csjy - bbm11(ny-1 )=asjy - bbm11(3 )=asjy - bbm11(ny-2 )=asjy+csjy - bbm11(4:ny-3)=asjy - if (istret==0) then - bbm11 = alsajy - xcst*bbm11 - else - bbm11(1:ny-1) = alsajy/pp2y(2:ny) - xcst*bbm11(1:ny-1) - endif - !BC on bbm11 - if (istret==0) then - bbm11(1 )=bbm11(1)+alsajy - else - bbm11(1 )=bbm11(1)+alsajy/pp2y(2) - endif - bbm11(ny)=zero - ! - !DIAG SUP 2 - ccm11(1 )=two*bsjy - ccm11(ny )=zero - ccm11(2 )=bsjy - ccm11(ny-1 )=zero - ccm11(3 )=bsjy - ccm11(ny-2 )=bsjy - ccm11(4:ny-3)=bsjy - ccm11 = -xcst*ccm11 - ! - !DIAG SUP 3 - rrm11(1 )=two*csjy - rrm11(ny )=zero - rrm11(2 )=csjy - rrm11(ny-1 )=zero - rrm11(3 )=csjy - rrm11(ny-2 )=zero - rrm11(4:ny-3)=csjy - rrm11 = -xcst*rrm11 - ! - !DIAG INF 1 - ddm11(1 )=zero - ddm11(ny )=two*asjy - ddm11(2 )=asjy - ddm11(ny-1 )=asjy+csjy - ddm11(3 )=asjy+csjy - ddm11(ny-2 )=asjy - ddm11(4:ny-3)=asjy - if (istret==0) then - ddm11 = alsajy - xcst*ddm11 - else - ddm11(2:ny) = alsajy/pp2y(1:ny-1) - xcst*ddm11(2:ny) - endif - !BC on ddm11 - ddm11(1 )=zero - if (istret==0) then - ddm11(ny)=ddm11(ny)+alsajy!a1 - else - ddm11(ny)=ddm11(ny)+alsajy/pp2y(ny-1)!a1 - endif - ! - !DIAG INF 2 - eem11(1 )=zero - eem11(ny )=two*bsjy - eem11(2 )=zero - eem11(ny-1 )=bsjy - eem11(3 )=bsjy - eem11(ny-2 )=bsjy - eem11(4:ny-3)=bsjy - eem11 = -xcst*eem11 - ! - !DIAG INF 3 - qqm11(1 )=zero - qqm11(ny )=two*csjy - qqm11(2 )=zero - qqm11(ny-1 )=csjy - qqm11(3 )=zero - qqm11(ny-2 )=csjy - qqm11(4:ny-3)=csjy - qqm11 = -xcst*qqm11 + ! + ! DIAG + aam11(1 )=-two*(asjy+bsjy+csjy) + aam11(ny )=aam11(1) + aam11(2 )=-two*asjy-bsjy-two*csjy + aam11(ny-1 )=aam11(2) + aam11(3:ny-2)=-two*(asjy+bsjy+csjy) + if (istret==0) then + aam11 = one - xcst*aam11 + else + aam11 = one/pp2y - xcst*aam11 + endif + ! + !DIAG SUP 1 + bbm11(1 )=two*asjy + bbm11(ny )=zero + bbm11(2 )=asjy+csjy + bbm11(ny-1 )=asjy + bbm11(3 )=asjy + bbm11(ny-2 )=asjy+csjy + bbm11(4:ny-3)=asjy + if (istret==0) then + bbm11 = alsajy - xcst*bbm11 + else + bbm11(1:ny-1) = alsajy/pp2y(2:ny) - xcst*bbm11(1:ny-1) + endif + !BC on bbm11 + if (istret==0) then + bbm11(1 )=bbm11(1)+alsajy + else + bbm11(1 )=bbm11(1)+alsajy/pp2y(2) + endif + bbm11(ny)=zero + ! + !DIAG SUP 2 + ccm11(1 )=two*bsjy + ccm11(ny )=zero + ccm11(2 )=bsjy + ccm11(ny-1 )=zero + ccm11(3 )=bsjy + ccm11(ny-2 )=bsjy + ccm11(4:ny-3)=bsjy + ccm11 = -xcst*ccm11 + ! + !DIAG SUP 3 + rrm11(1 )=two*csjy + rrm11(ny )=zero + rrm11(2 )=csjy + rrm11(ny-1 )=zero + rrm11(3 )=csjy + rrm11(ny-2 )=zero + rrm11(4:ny-3)=csjy + rrm11 = -xcst*rrm11 + ! + !DIAG INF 1 + ddm11(1 )=zero + ddm11(ny )=two*asjy + ddm11(2 )=asjy + ddm11(ny-1 )=asjy+csjy + ddm11(3 )=asjy+csjy + ddm11(ny-2 )=asjy + ddm11(4:ny-3)=asjy + if (istret==0) then + ddm11 = alsajy - xcst*ddm11 + else + ddm11(2:ny) = alsajy/pp2y(1:ny-1) - xcst*ddm11(2:ny) + endif + !BC on ddm11 + ddm11(1 )=zero + if (istret==0) then + ddm11(ny)=ddm11(ny)+alsajy!a1 + else + ddm11(ny)=ddm11(ny)+alsajy/pp2y(ny-1)!a1 + endif + ! + !DIAG INF 2 + eem11(1 )=zero + eem11(ny )=two*bsjy + eem11(2 )=zero + eem11(ny-1 )=bsjy + eem11(3 )=bsjy + eem11(ny-2 )=bsjy + eem11(4:ny-3)=bsjy + eem11 = -xcst*eem11 + ! + !DIAG INF 3 + qqm11(1 )=zero + qqm11(ny )=two*csjy + qqm11(2 )=zero + qqm11(ny-1 )=csjy + qqm11(3 )=zero + qqm11(ny-2 )=csjy + qqm11(4:ny-3)=csjy + qqm11 = -xcst*qqm11 !!! NXL = 0 - !DIAG - if (istret==0) then - aam0 = one-xcst*(-two*(asjy+bsjy+csjy)) - else - aam0 = one/pp2y-xcst*(-two*(asjy+bsjy+csjy)) - endif - ! - !DIAG SUP 1 - if (istret==0) then - bbm0 = alsajy-xcst*asjy - else - bbm0(1:ny-1) = alsajy/pp2y(2:ny) -xcst*asjy - bbm0(ny) = alsajy/pp2y(1) -xcst*asjy - endif - ! - !DIAG SUP 2 - ccm0 = -xcst*bsjy - ! - !DIAG SUP 3 - rrm0 = -xcst*csjy - ! - !DIAG INF 1 - if (istret==0) then - ddm0=bbm0 - else - ddm0(1) = alsajy/pp2y(ny) -xcst*asjy - ddm0(2:ny) = alsajy/pp2y(1:ny-1) -xcst*asjy - endif - ! - !DIAG INF 2 - eem0=ccm0 - ! - !DIAG INF 3 - qqm0=rrm0 - endif - - if (isecondder.ne.5) then - ! velocity, ncly1 = 2, nclyn = 2 - call ludecomp7(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,& - vvm,wwm,zzm,ny) - ! velocity, ncly1 = 1, nclyn = 1, npaire = 0 - call ludecomp7(aam10,bbm10,ccm10,ddm10,eem10,qqm10,ggm10,hhm10,ssm10,rrm10,& - vvm10,wwm10,zzm10,ny) - ! velocity, ncly1 = 1, nclyn = 1, npaire = 1 - call ludecomp7(aam11,bbm11,ccm11,ddm11,eem11,qqm11,ggm11,hhm11,ssm11,rrm11,& - vvm11,wwm11,zzm11,ny) - ! velocity, ncly1 = 0, nclyn = 0 - call ludecomp7(aam0,bbm0,ccm0,ddm0,eem0,qqm0,ggm0,hhm0,ssm0,rrm0,& - vvm0,wwm0,zzm0,l1m,l2m,l3m,u1m,u2m,u3m,ny) - ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 - call ludecomp7(aam120,bbm120,ccm120,ddm120,eem120,qqm120,ggm120,hhm120,ssm120,rrm120,& - vvm120,wwm120,zzm120,ny) - ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 - call ludecomp7(aam121,bbm121,ccm121,ddm121,eem121,qqm121,ggm121,hhm121,ssm121,rrm121,& - vvm121,wwm121,zzm121,ny) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 - call ludecomp7(aam210,bbm210,ccm210,ddm210,eem210,qqm210,ggm210,hhm210,ssm210,rrm210,& - vvm210,wwm210,zzm210,ny) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 - call ludecomp7(aam211,bbm211,ccm211,ddm211,eem211,qqm211,ggm211,hhm211,ssm211,rrm211,& - vvm211,wwm211,zzm211,ny) - ! Scalars - do is = 1, numscalar - ! scalar, ncly1 = 2, nclyn = 2 - call ludecomp7(aamt(:,is),bbmt(:,is),ccmt(:,is),ddmt(:,is),eemt(:,is),& - qqmt(:,is),ggmt(:,is),hhmt(:,is),ssmt(:,is),rrmt(:,is),& - vvmt(:,is),wwmt(:,is),zzmt(:,is),ny) - ! scalar, ncly1 = 1, nclyn = 1, npaire = 0 - call ludecomp7(aam10t(:,is),bbm10t(:,is),ccm10t(:,is),ddm10t(:,is),eem10t(:,is),& - qqm10t(:,is),ggm10t(:,is),hhm10t(:,is),ssm10t(:,is),rrm10t(:,is),& - vvm10t(:,is),wwm10t(:,is),zzm10t(:,is),ny) - ! scalar, ncly1 = 1, nclyn = 1, npaire = 1 - call ludecomp7(aam11t(:,is),bbm11t(:,is),ccm11t(:,is),ddm11t(:,is),eem11t(:,is),& - qqm11t(:,is),ggm11t(:,is),hhm11t(:,is),ssm11t(:,is),rrm11t(:,is),& - vvm11t(:,is),wwm11t(:,is),zzm11t(:,is),ny) - ! scalar, ncly1 = 0, nclyn = 0 - call ludecomp7(aam0t(:,is),bbm0t(:,is),ccm0t(:,is),ddm0t(:,is),eem0t(:,is),& - qqm0t(:,is),ggm0t(:,is),hhm0t(:,is),ssm0t(:,is),rrm0t(:,is),& - vvm0t(:,is),wwm0t(:,is),zzm0t(:,is),l1mt(:,is),l2mt(:,is),l3mt(:,is),u1mt(:,is),u2mt(:,is),u3mt(:,is),ny) - ! scalar, ncly1 = 1, nclyn = 2, npaire = 0 - call ludecomp7(aam120t(:,is),bbm120t(:,is),ccm120t(:,is),ddm120t(:,is),eem120t(:,is),& - qqm120t(:,is),ggm120t(:,is),hhm120t(:,is),ssm120t(:,is),rrm120t(:,is),& - vvm120t(:,is),wwm120t(:,is),zzm120t(:,is),ny) - ! scalar, ncly1 = 1, nclyn = 2, npaire = 1 - call ludecomp7(aam121t(:,is),bbm121t(:,is),ccm121t(:,is),ddm121t(:,is),eem121t(:,is),& - qqm121t(:,is),ggm121t(:,is),hhm121t(:,is),ssm121t(:,is),rrm121t(:,is),& - vvm121t(:,is),wwm121t(:,is),zzm121t(:,is),ny) - ! scalar, ncly1 = 2, nclyn = 1, npaire = 0 - call ludecomp7(aam210t(:,is),bbm210t(:,is),ccm210t(:,is),ddm210t(:,is),eem210t(:,is),& - qqm210t(:,is),ggm210t(:,is),hhm210t(:,is),ssm210t(:,is),rrm210t(:,is),& - vvm210t(:,is),wwm210t(:,is),zzm210t(:,is),ny) - ! scalar, ncly1 = 2, nclyn = 1, npaire = 1 - call ludecomp7(aam211t(:,is),bbm211t(:,is),ccm211t(:,is),ddm211t(:,is),eem211t(:,is),& - qqm211t(:,is),ggm211t(:,is),hhm211t(:,is),ssm211t(:,is),rrm211t(:,is),& - vvm211t(:,is),wwm211t(:,is),zzm211t(:,is),ny) - enddo - else - call ludecomp9(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ttm,uum,sssm,zzzm,ny) - !NEED TO BE DONE: deal with other cases - endif + !DIAG + if (istret==0) then + aam0 = one-xcst*(-two*(asjy+bsjy+csjy)) + else + aam0 = one/pp2y-xcst*(-two*(asjy+bsjy+csjy)) + endif + ! + !DIAG SUP 1 + if (istret==0) then + bbm0 = alsajy-xcst*asjy + else + bbm0(1:ny-1) = alsajy/pp2y(2:ny) -xcst*asjy + bbm0(ny) = alsajy/pp2y(1) -xcst*asjy + endif + ! + !DIAG SUP 2 + ccm0 = -xcst*bsjy + ! + !DIAG SUP 3 + rrm0 = -xcst*csjy + ! + !DIAG INF 1 + if (istret==0) then + ddm0=bbm0 + else + ddm0(1) = alsajy/pp2y(ny) -xcst*asjy + ddm0(2:ny) = alsajy/pp2y(1:ny-1) -xcst*asjy + endif + ! + !DIAG INF 2 + eem0=ccm0 + ! + !DIAG INF 3 + qqm0=rrm0 + endif + + if (isecondder.ne.5) then + ! velocity, ncly1 = 2, nclyn = 2 + call ludecomp7(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,& + vvm,wwm,zzm,ny) + ! velocity, ncly1 = 1, nclyn = 1, npaire = 0 + call ludecomp7(aam10,bbm10,ccm10,ddm10,eem10,qqm10,ggm10,hhm10,ssm10,rrm10,& + vvm10,wwm10,zzm10,ny) + ! velocity, ncly1 = 1, nclyn = 1, npaire = 1 + call ludecomp7(aam11,bbm11,ccm11,ddm11,eem11,qqm11,ggm11,hhm11,ssm11,rrm11,& + vvm11,wwm11,zzm11,ny) + ! velocity, ncly1 = 0, nclyn = 0 + call ludecomp7(aam0,bbm0,ccm0,ddm0,eem0,qqm0,ggm0,hhm0,ssm0,rrm0,& + vvm0,wwm0,zzm0,l1m,l2m,l3m,u1m,u2m,u3m,ny) + ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 + call ludecomp7(aam120,bbm120,ccm120,ddm120,eem120,qqm120,ggm120,hhm120,ssm120,rrm120,& + vvm120,wwm120,zzm120,ny) + ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 + call ludecomp7(aam121,bbm121,ccm121,ddm121,eem121,qqm121,ggm121,hhm121,ssm121,rrm121,& + vvm121,wwm121,zzm121,ny) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 + call ludecomp7(aam210,bbm210,ccm210,ddm210,eem210,qqm210,ggm210,hhm210,ssm210,rrm210,& + vvm210,wwm210,zzm210,ny) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 + call ludecomp7(aam211,bbm211,ccm211,ddm211,eem211,qqm211,ggm211,hhm211,ssm211,rrm211,& + vvm211,wwm211,zzm211,ny) + ! Scalars + do is = 1, numscalar + ! scalar, ncly1 = 2, nclyn = 2 + call ludecomp7(aamt(:,is),bbmt(:,is),ccmt(:,is),ddmt(:,is),eemt(:,is),& + qqmt(:,is),ggmt(:,is),hhmt(:,is),ssmt(:,is),rrmt(:,is),& + vvmt(:,is),wwmt(:,is),zzmt(:,is),ny) + ! scalar, ncly1 = 1, nclyn = 1, npaire = 0 + call ludecomp7(aam10t(:,is),bbm10t(:,is),ccm10t(:,is),ddm10t(:,is),eem10t(:,is),& + qqm10t(:,is),ggm10t(:,is),hhm10t(:,is),ssm10t(:,is),rrm10t(:,is),& + vvm10t(:,is),wwm10t(:,is),zzm10t(:,is),ny) + ! scalar, ncly1 = 1, nclyn = 1, npaire = 1 + call ludecomp7(aam11t(:,is),bbm11t(:,is),ccm11t(:,is),ddm11t(:,is),eem11t(:,is),& + qqm11t(:,is),ggm11t(:,is),hhm11t(:,is),ssm11t(:,is),rrm11t(:,is),& + vvm11t(:,is),wwm11t(:,is),zzm11t(:,is),ny) + ! scalar, ncly1 = 0, nclyn = 0 + call ludecomp7(aam0t(:,is),bbm0t(:,is),ccm0t(:,is),ddm0t(:,is),eem0t(:,is),& + qqm0t(:,is),ggm0t(:,is),hhm0t(:,is),ssm0t(:,is),rrm0t(:,is),& + vvm0t(:,is),wwm0t(:,is),zzm0t(:,is),l1mt(:,is),l2mt(:,is),l3mt(:,is),u1mt(:,is),u2mt(:,is),u3mt(:,is),ny) + ! scalar, ncly1 = 1, nclyn = 2, npaire = 0 + call ludecomp7(aam120t(:,is),bbm120t(:,is),ccm120t(:,is),ddm120t(:,is),eem120t(:,is),& + qqm120t(:,is),ggm120t(:,is),hhm120t(:,is),ssm120t(:,is),rrm120t(:,is),& + vvm120t(:,is),wwm120t(:,is),zzm120t(:,is),ny) + ! scalar, ncly1 = 1, nclyn = 2, npaire = 1 + call ludecomp7(aam121t(:,is),bbm121t(:,is),ccm121t(:,is),ddm121t(:,is),eem121t(:,is),& + qqm121t(:,is),ggm121t(:,is),hhm121t(:,is),ssm121t(:,is),rrm121t(:,is),& + vvm121t(:,is),wwm121t(:,is),zzm121t(:,is),ny) + ! scalar, ncly1 = 2, nclyn = 1, npaire = 0 + call ludecomp7(aam210t(:,is),bbm210t(:,is),ccm210t(:,is),ddm210t(:,is),eem210t(:,is),& + qqm210t(:,is),ggm210t(:,is),hhm210t(:,is),ssm210t(:,is),rrm210t(:,is),& + vvm210t(:,is),wwm210t(:,is),zzm210t(:,is),ny) + ! scalar, ncly1 = 2, nclyn = 1, npaire = 1 + call ludecomp7(aam211t(:,is),bbm211t(:,is),ccm211t(:,is),ddm211t(:,is),eem211t(:,is),& + qqm211t(:,is),ggm211t(:,is),hhm211t(:,is),ssm211t(:,is),rrm211t(:,is),& + vvm211t(:,is),wwm211t(:,is),zzm211t(:,is),ny) + enddo + do is = 1, 3 + call ludecomp7(aamB(:,is),bbmB(:,is),ccmB(:,is),ddmB(:,is),eemB(:,is),qqmB(:,is),ggmB(:,is),hhmB(:,is),ssmB(:,is),rrmB(:,is),& + vvmB(:,is),wwmB(:,is),zzmB(:,is),ny) + enddo + else + call ludecomp9(aam,bbm,ccm,ddm,eem,qqm,ggm,hhm,ssm,rrm,vvm,wwm,zzm,ttm,uum,sssm,zzzm,ny) + !NEED TO BE DONE: deal with other cases + endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - !! END MATRIX M2 TIME IMPLICIT !! + !! END MATRIX M2 TIME IMPLICIT !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -end subroutine implicit_schemes + end subroutine implicit_schemes -! -! Allocate 1D arrays containing LU decompositions -! -subroutine init_implicit() + ! + ! Allocate 1D arrays containing LU decompositions + ! + subroutine init_implicit() - USE decomp_2d - USE param - USE variables - implicit none + USE decomp_2d + USE param + USE variables + implicit none - ! velocity, ncly1 = 2, nclyn = 2 - allocate(aam(ny),bbm(ny),ccm(ny),ddm(ny),eem(ny),ggm(ny),hhm(ny),wwm(ny),zzm(ny)) - allocate(rrm(ny),qqm(ny),vvm(ny),ssm(ny)) - allocate(sssm(ny),zzzm(ny),ttm(ny),uum(ny)) ! nona - ! velocity, ncly1 = 1, nclyn = 1, npaire = 0 - allocate(aam10(ny),bbm10(ny),ccm10(ny),ddm10(ny),eem10(ny),ggm10(ny),hhm10(ny),wwm10(ny),zzm10(ny)) - allocate(rrm10(ny),qqm10(ny),vvm10(ny),ssm10(ny)) - ! velocity, ncly1 = 1, nclyn = 1, npaire = 1 - allocate(aam11(ny),bbm11(ny),ccm11(ny),ddm11(ny),eem11(ny),ggm11(ny),hhm11(ny),wwm11(ny),zzm11(ny)) - allocate(rrm11(ny),qqm11(ny),vvm11(ny),ssm11(ny)) - ! velocity, ncly1 = 0, nclyn = 0 - allocate(aam0(ny),bbm0(ny),ccm0(ny),ddm0(ny),eem0(ny),ggm0(ny),hhm0(ny),wwm0(ny),zzm0(ny)) - allocate(rrm0(ny),qqm0(ny),vvm0(ny),ssm0(ny),l1m(ny),l2m(ny),l3m(ny),u1m(ny),u2m(ny),u3m(ny)) - ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 - allocate(aam120(ny),bbm120(ny),ccm120(ny),ddm120(ny),eem120(ny),ggm120(ny),hhm120(ny),wwm120(ny),zzm120(ny)) - allocate(rrm120(ny),qqm120(ny),vvm120(ny),ssm120(ny)) - ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 - allocate(aam121(ny),bbm121(ny),ccm121(ny),ddm121(ny),eem121(ny),ggm121(ny),hhm121(ny),wwm121(ny),zzm121(ny)) - allocate(rrm121(ny),qqm121(ny),vvm121(ny),ssm121(ny)) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 - allocate(aam210(ny),bbm210(ny),ccm210(ny),ddm210(ny),eem210(ny),ggm210(ny),hhm210(ny),wwm210(ny),zzm210(ny)) - allocate(rrm210(ny),qqm210(ny),vvm210(ny),ssm210(ny)) - ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 - allocate(aam211(ny),bbm211(ny),ccm211(ny),ddm211(ny),eem211(ny),ggm211(ny),hhm211(ny),wwm211(ny),zzm211(ny)) - allocate(rrm211(ny),qqm211(ny),vvm211(ny),ssm211(ny)) - ! scalar, ncly1 = 2, nclyn = 2 - allocate(aamt(ny,numscalar),bbmt(ny,numscalar),ccmt(ny,numscalar),ddmt(ny,numscalar),eemt(ny,numscalar)) - allocate(ggmt(ny,numscalar),hhmt(ny,numscalar),wwmt(ny,numscalar),zzmt(ny,numscalar)) - allocate(rrmt(ny,numscalar),qqmt(ny,numscalar),vvmt(ny,numscalar),ssmt(ny,numscalar)) - allocate(uumt(ny,numscalar),ttmt(ny,numscalar),sssmt(ny,numscalar),zzzmt(ny,numscalar)) ! nona - ! scalar, ncly1 = 1, nclyn = 1, npaire = 0 - allocate(aam10t(ny,numscalar),bbm10t(ny,numscalar),ccm10t(ny,numscalar),ddm10t(ny,numscalar),eem10t(ny,numscalar)) - allocate(ggm10t(ny,numscalar),hhm10t(ny,numscalar),wwm10t(ny,numscalar),zzm10t(ny,numscalar)) - allocate(rrm10t(ny,numscalar),qqm10t(ny,numscalar),vvm10t(ny,numscalar),ssm10t(ny,numscalar)) - ! scalar, ncly1 = 1, nclyn = 1, npaire = 1 - allocate(aam11t(ny,numscalar),bbm11t(ny,numscalar),ccm11t(ny,numscalar),ddm11t(ny,numscalar),eem11t(ny,numscalar)) - allocate(ggm11t(ny,numscalar),hhm11t(ny,numscalar),wwm11t(ny,numscalar),zzm11t(ny,numscalar)) - allocate(rrm11t(ny,numscalar),qqm11t(ny,numscalar),vvm11t(ny,numscalar),ssm11t(ny,numscalar)) - ! scalar, ncly1 = 0, nclyn = 0 - allocate(aam0t(ny,numscalar),bbm0t(ny,numscalar),ccm0t(ny,numscalar),ddm0t(ny,numscalar),eem0t(ny,numscalar)) - allocate(ggm0t(ny,numscalar),hhm0t(ny,numscalar),wwm0t(ny,numscalar),zzm0t(ny,numscalar)) - allocate(rrm0t(ny,numscalar),qqm0t(ny,numscalar),vvm0t(ny,numscalar),ssm0t(ny,numscalar)) - allocate(l1mt(ny,numscalar),l2mt(ny,numscalar),l3mt(ny,numscalar),u1mt(ny,numscalar),u2mt(ny,numscalar),u3mt(ny,numscalar)) - ! scalar, ncly1 = 1, nclyn = 2, npaire = 0 - allocate(aam120t(ny,numscalar),bbm120t(ny,numscalar),ccm120t(ny,numscalar),ddm120t(ny,numscalar),eem120t(ny,numscalar)) - allocate(ggm120t(ny,numscalar),hhm120t(ny,numscalar),wwm120t(ny,numscalar),zzm120t(ny,numscalar)) - allocate(rrm120t(ny,numscalar),qqm120t(ny,numscalar),vvm120t(ny,numscalar),ssm120t(ny,numscalar)) - ! scalar, ncly1 = 1, nclyn = 2, npaire = 1 - allocate(aam121t(ny,numscalar),bbm121t(ny,numscalar),ccm121t(ny,numscalar),ddm121t(ny,numscalar),eem121t(ny,numscalar)) - allocate(ggm121t(ny,numscalar),hhm121t(ny,numscalar),wwm121t(ny,numscalar),zzm121t(ny,numscalar)) - allocate(rrm121t(ny,numscalar),qqm121t(ny,numscalar),vvm121t(ny,numscalar),ssm121t(ny,numscalar)) - ! scalar, ncly1 = 2, nclyn = 1, npaire = 0 - allocate(aam210t(ny,numscalar),bbm210t(ny,numscalar),ccm210t(ny,numscalar),ddm210t(ny,numscalar),eem210t(ny,numscalar)) - allocate(ggm210t(ny,numscalar),hhm210t(ny,numscalar),wwm210t(ny,numscalar),zzm210t(ny,numscalar)) - allocate(rrm210t(ny,numscalar),qqm210t(ny,numscalar),vvm210t(ny,numscalar),ssm210t(ny,numscalar)) - ! scalar, ncly1 = 2, nclyn = 1, npaire = 1 - allocate(aam211t(ny,numscalar),bbm211t(ny,numscalar),ccm211t(ny,numscalar),ddm211t(ny,numscalar),eem211t(ny,numscalar)) - allocate(ggm211t(ny,numscalar),hhm211t(ny,numscalar),wwm211t(ny,numscalar),zzm211t(ny,numscalar)) - allocate(rrm211t(ny,numscalar),qqm211t(ny,numscalar),vvm211t(ny,numscalar),ssm211t(ny,numscalar)) - -end subroutine init_implicit + ! velocity, ncly1 = 2, nclyn = 2 + allocate(aam(ny),bbm(ny),ccm(ny),ddm(ny),eem(ny),ggm(ny),hhm(ny),wwm(ny),zzm(ny)) + allocate(rrm(ny),qqm(ny),vvm(ny),ssm(ny)) + allocate(sssm(ny),zzzm(ny),ttm(ny),uum(ny)) ! nona + ! velocity, ncly1 = 1, nclyn = 1, npaire = 0 + allocate(aam10(ny),bbm10(ny),ccm10(ny),ddm10(ny),eem10(ny),ggm10(ny),hhm10(ny),wwm10(ny),zzm10(ny)) + allocate(rrm10(ny),qqm10(ny),vvm10(ny),ssm10(ny)) + ! velocity, ncly1 = 1, nclyn = 1, npaire = 1 + allocate(aam11(ny),bbm11(ny),ccm11(ny),ddm11(ny),eem11(ny),ggm11(ny),hhm11(ny),wwm11(ny),zzm11(ny)) + allocate(rrm11(ny),qqm11(ny),vvm11(ny),ssm11(ny)) + ! velocity, ncly1 = 0, nclyn = 0 + allocate(aam0(ny),bbm0(ny),ccm0(ny),ddm0(ny),eem0(ny),ggm0(ny),hhm0(ny),wwm0(ny),zzm0(ny)) + allocate(rrm0(ny),qqm0(ny),vvm0(ny),ssm0(ny),l1m(ny),l2m(ny),l3m(ny),u1m(ny),u2m(ny),u3m(ny)) + ! velocity, ncly1 = 1, nclyn = 2, npaire = 0 + allocate(aam120(ny),bbm120(ny),ccm120(ny),ddm120(ny),eem120(ny),ggm120(ny),hhm120(ny),wwm120(ny),zzm120(ny)) + allocate(rrm120(ny),qqm120(ny),vvm120(ny),ssm120(ny)) + ! velocity, ncly1 = 1, nclyn = 2, npaire = 1 + allocate(aam121(ny),bbm121(ny),ccm121(ny),ddm121(ny),eem121(ny),ggm121(ny),hhm121(ny),wwm121(ny),zzm121(ny)) + allocate(rrm121(ny),qqm121(ny),vvm121(ny),ssm121(ny)) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 0 + allocate(aam210(ny),bbm210(ny),ccm210(ny),ddm210(ny),eem210(ny),ggm210(ny),hhm210(ny),wwm210(ny),zzm210(ny)) + allocate(rrm210(ny),qqm210(ny),vvm210(ny),ssm210(ny)) + ! velocity, ncly1 = 2, nclyn = 1, npaire = 1 + allocate(aam211(ny),bbm211(ny),ccm211(ny),ddm211(ny),eem211(ny),ggm211(ny),hhm211(ny),wwm211(ny),zzm211(ny)) + allocate(rrm211(ny),qqm211(ny),vvm211(ny),ssm211(ny)) + ! scalar, ncly1 = 2, nclyn = 2 + allocate(aamt(ny,numscalar),bbmt(ny,numscalar),ccmt(ny,numscalar),ddmt(ny,numscalar),eemt(ny,numscalar)) + allocate(ggmt(ny,numscalar),hhmt(ny,numscalar),wwmt(ny,numscalar),zzmt(ny,numscalar)) + allocate(rrmt(ny,numscalar),qqmt(ny,numscalar),vvmt(ny,numscalar),ssmt(ny,numscalar)) + allocate(uumt(ny,numscalar),ttmt(ny,numscalar),sssmt(ny,numscalar),zzzmt(ny,numscalar)) ! nona + ! scalar, ncly1 = 1, nclyn = 1, npaire = 0 + allocate(aam10t(ny,numscalar),bbm10t(ny,numscalar),ccm10t(ny,numscalar),ddm10t(ny,numscalar),eem10t(ny,numscalar)) + allocate(ggm10t(ny,numscalar),hhm10t(ny,numscalar),wwm10t(ny,numscalar),zzm10t(ny,numscalar)) + allocate(rrm10t(ny,numscalar),qqm10t(ny,numscalar),vvm10t(ny,numscalar),ssm10t(ny,numscalar)) + ! scalar, ncly1 = 1, nclyn = 1, npaire = 1 + allocate(aam11t(ny,numscalar),bbm11t(ny,numscalar),ccm11t(ny,numscalar),ddm11t(ny,numscalar),eem11t(ny,numscalar)) + allocate(ggm11t(ny,numscalar),hhm11t(ny,numscalar),wwm11t(ny,numscalar),zzm11t(ny,numscalar)) + allocate(rrm11t(ny,numscalar),qqm11t(ny,numscalar),vvm11t(ny,numscalar),ssm11t(ny,numscalar)) + ! scalar, ncly1 = 0, nclyn = 0 + allocate(aam0t(ny,numscalar),bbm0t(ny,numscalar),ccm0t(ny,numscalar),ddm0t(ny,numscalar),eem0t(ny,numscalar)) + allocate(ggm0t(ny,numscalar),hhm0t(ny,numscalar),wwm0t(ny,numscalar),zzm0t(ny,numscalar)) + allocate(rrm0t(ny,numscalar),qqm0t(ny,numscalar),vvm0t(ny,numscalar),ssm0t(ny,numscalar)) + allocate(l1mt(ny,numscalar),l2mt(ny,numscalar),l3mt(ny,numscalar),u1mt(ny,numscalar),u2mt(ny,numscalar),u3mt(ny,numscalar)) + ! scalar, ncly1 = 1, nclyn = 2, npaire = 0 + allocate(aam120t(ny,numscalar),bbm120t(ny,numscalar),ccm120t(ny,numscalar),ddm120t(ny,numscalar),eem120t(ny,numscalar)) + allocate(ggm120t(ny,numscalar),hhm120t(ny,numscalar),wwm120t(ny,numscalar),zzm120t(ny,numscalar)) + allocate(rrm120t(ny,numscalar),qqm120t(ny,numscalar),vvm120t(ny,numscalar),ssm120t(ny,numscalar)) + ! scalar, ncly1 = 1, nclyn = 2, npaire = 1 + allocate(aam121t(ny,numscalar),bbm121t(ny,numscalar),ccm121t(ny,numscalar),ddm121t(ny,numscalar),eem121t(ny,numscalar)) + allocate(ggm121t(ny,numscalar),hhm121t(ny,numscalar),wwm121t(ny,numscalar),zzm121t(ny,numscalar)) + allocate(rrm121t(ny,numscalar),qqm121t(ny,numscalar),vvm121t(ny,numscalar),ssm121t(ny,numscalar)) + ! scalar, ncly1 = 2, nclyn = 1, npaire = 0 + allocate(aam210t(ny,numscalar),bbm210t(ny,numscalar),ccm210t(ny,numscalar),ddm210t(ny,numscalar),eem210t(ny,numscalar)) + allocate(ggm210t(ny,numscalar),hhm210t(ny,numscalar),wwm210t(ny,numscalar),zzm210t(ny,numscalar)) + allocate(rrm210t(ny,numscalar),qqm210t(ny,numscalar),vvm210t(ny,numscalar),ssm210t(ny,numscalar)) + ! scalar, ncly1 = 2, nclyn = 1, npaire = 1 + allocate(aam211t(ny,numscalar),bbm211t(ny,numscalar),ccm211t(ny,numscalar),ddm211t(ny,numscalar),eem211t(ny,numscalar)) + allocate(ggm211t(ny,numscalar),hhm211t(ny,numscalar),wwm211t(ny,numscalar),zzm211t(ny,numscalar)) + allocate(rrm211t(ny,numscalar),qqm211t(ny,numscalar),vvm211t(ny,numscalar),ssm211t(ny,numscalar)) + + allocate(aamB(ny,3),bbmB(ny,3),ccmB(ny,3),ddmB(ny,3),eemB(ny,3),ggmB(ny,3),hhmB(ny,3),wwmB(ny,3),zzmB(ny,3)) + allocate(rrmB(ny,3),qqmB(ny,3),vvmB(ny,3),ssmB(ny,3)) + allocate(sssmB(ny,3),zzzmB(ny,3),ttmB(ny,3),uumB(ny,3)) ! nona + + + + end subroutine init_implicit -! -! Used to build the scalar implicit coefficients -! -subroutine init_implicit_coef(tab1d, tab2d) + ! + ! Used to build the scalar implicit coefficients + ! + subroutine init_implicit_coef(tab1d, tab2d) - use decomp_2d, only : mytype - use variables, only : ny, numscalar + use decomp_2d, only : mytype + use variables, only : ny, numscalar - implicit none + implicit none - real(mytype), dimension(ny), intent(in) :: tab1d - real(mytype), dimension(ny,numscalar), intent(out) :: tab2d + real(mytype), dimension(ny), intent(in) :: tab1d + real(mytype), dimension(ny,numscalar), intent(out) :: tab2d - integer :: is + integer :: is - do is = 1, numscalar - tab2d(:,is) = tab1d(:) - enddo + do is = 1, numscalar + tab2d(:,is) = tab1d(:) + enddo -end subroutine init_implicit_coef + end subroutine init_implicit_coef end module ydiff_implicit diff --git a/src/mhd.f90 b/src/mhd.f90 index 7f14a4fa8..c01a581f2 100644 --- a/src/mhd.f90 +++ b/src/mhd.f90 @@ -92,21 +92,45 @@ subroutine int_time_magnet USE param USE variables USE decomp_2d + use ydiff_implicit implicit none ! integer :: k - if(itimescheme.eq.5) then - ! - if(itr.eq.1) then - Bm=gdt(itr)*dBm(:,:,:,:,1)+Bm - else - Bm=adt(itr)*dBm(:,:,:,:,1)+bdt(itr)*dBm(:,:,:,:,2)+Bm + if( iimplicit == 1 ) then + call inttimp(Bm(:,:,:,1), dBm(:,:,:,1,:), 0, -1, mhdvar=1 ) + call inttimp(Bm(:,:,:,2), dBm(:,:,:,2,:), 0, -1, mhdvar=2 ) + call inttimp(Bm(:,:,:,3), dBm(:,:,:,3,:), 0, -1, mhdvar=3 ) + else + if(itimescheme.eq.3) then + !>>> Adams-Bashforth third order (AB3) + + ! Do first time step with Euler + if(itime.eq.1.and.irestart.eq.0) then + Bm=dt*dBm(:,:,:,:,1)+Bm + elseif(itime.eq.2.and.irestart.eq.0) then + ! Do second time step with AB2 + Bm=onepfive*dt*dBm(:,:,:,:,1)-half*dt*dBm(:,:,:,:,2)+Bm + dBm(:,:,:,:,3)=dBm(:,:,:,:,2) + else + ! Finally using AB3 + Bm=adt(itr)*dBm(:,:,:,:,1)+bdt(itr)*dBm(:,:,:,:,2)+cdt(itr)*dBm(:,:,:,:,3)+Bm + dBm(:,:,:,:,3)=dBm(:,:,:,:,2) + endif + dBm(:,:,:,:,2)=dBm(:,:,:,:,1) + + elseif(itimescheme.eq.5) then + ! + if(itr.eq.1) then + Bm=gdt(itr)*dBm(:,:,:,:,1)+Bm + else + Bm=adt(itr)*dBm(:,:,:,:,1)+bdt(itr)*dBm(:,:,:,:,2)+Bm + endif + dBm(:,:,:,:,2)=dBm(:,:,:,:,1) + ! endif - dBm(:,:,:,:,2)=dBm(:,:,:,:,1) - ! - endif + endif ! end subroutine int_time_magnet ! @@ -711,33 +735,12 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) rrem=1.d0/Rem - !SKEW SYMMETRIC FORM - !WORK X-PENCILS - ta1(:,:,:) = ux1(:,:,:) * B(:,:,:,1) - B(:,:,:,1) * ux1(:,:,:) - tb1(:,:,:) = ux1(:,:,:) * B(:,:,:,2) - B(:,:,:,1) * uy1(:,:,:) - tc1(:,:,:) = ux1(:,:,:) * B(:,:,:,3) - B(:,:,:,1) * uz1(:,:,:) - - call derx (td1,ta1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcx*ubcx) - call derx (te1,tb1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx*ubcy) - call derx (tf1,tc1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx*ubcz) - - call derx (ta1,ux1,di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx) - call derx (tb1,uy1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcy) - call derx (tc1,uz1,di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcz) - ! - call derx (tx1,B(:,:,:,1),di1,sx,ffx,fsx,fwx,xsize(1),xsize(2),xsize(3),0,ubcx) - call derx (ty1,B(:,:,:,2),di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcy) - call derx (tz1,B(:,:,:,3),di1,sx,ffxp,fsxp,fwxp,xsize(1),xsize(2),xsize(3),1,ubcz) + tb1(:,:,:) = ux1(:,:,:) * B(:,:,:,3) - B(:,:,:,1) * uz1(:,:,:) + tc1(:,:,:) = ux1(:,:,:) * B(:,:,:,2) - B(:,:,:,1) * uy1(:,:,:) - - ! Convective terms of x-pencil are stored in tg1,th1,ti1 - - tg1(:,:,:) = td1(:,:,:) + ux1(:,:,:) * tx1(:,:,:) - B(:,:,:,1) * ta1(:,:,:) - th1(:,:,:) = te1(:,:,:) + ux1(:,:,:) * ty1(:,:,:) - B(:,:,:,1) * tb1(:,:,:) - ti1(:,:,:) = tf1(:,:,:) + ux1(:,:,:) * tz1(:,:,:) - B(:,:,:,1) * tc1(:,:,:) - - ! TODO: save the x-convective terms already in dux1, duy1, duz1 + call derxBy (te1,tb1,di1,sx,ffxB(:,2),fsxB(:,2),fwxB(:,2),xsize(1),xsize(2),xsize(3),0,ubcx*ubcy) + call derxBz (tf1,tc1,di1,sx,ffxB(:,3),fsxB(:,3),fwxB(:,3),xsize(1),xsize(2),xsize(3),0,ubcx*ubcz) call transpose_x_to_y(ux1,ux2) call transpose_x_to_y(uy1,uy2) @@ -749,28 +752,12 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) !WORK Y-PENCILS - td2(:,:,:) = uy2(:,:,:)*bx2(:,:,:) - by2(:,:,:)*ux2(:,:,:) - te2(:,:,:) = uy2(:,:,:)*by2(:,:,:) - by2(:,:,:)*uy2(:,:,:) - tf2(:,:,:) = uy2(:,:,:)*bz2(:,:,:) - by2(:,:,:)*uz2(:,:,:) - - call dery (tg2,td2,di2,sy,ffy, fsy, fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcx*ubcy) - call dery (th2,te2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcy*ubcy) - call dery (ti2,tf2,di2,sy,ffy, fsy, fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcz*ubcy) - - call dery (td2,ux2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) - call dery (te2,uy2,di2,sy,ffy, fsy ,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy) - call dery (tf2,uz2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcz) - - call dery (tx2,bx2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) - call dery (ty2,by2,di2,sy,ffy, fsy ,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy) - call dery (tz2,bz2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcz) - + ta2(:,:,:) = uy2(:,:,:) * Bz2(:,:,:) - by2 * uz2(:,:,:) + tc2(:,:,:) = ux2(:,:,:) * by2(:,:,:) - bx2 * uy2(:,:,:) - ! Convective terms of y-pencil in tg2,th2,ti2 + call deryBx (td2,ta2,di2,sy,ffyB(:,1), fsyB(:,1), fwyB(:,1),ppy,ysize(1),ysize(2),ysize(3),0,ubcx*ubcy) + call deryBz (tf2,tc2,di2,sy,ffyB(:,3), fsyB(:,3), fwyB(:,3),ppy,ysize(1),ysize(2),ysize(3),0,ubcz*ubcy) - tg2(:,:,:) = tg2(:,:,:) + uy2(:,:,:) * tx2(:,:,:) - by2(:,:,:) * td2(:,:,:) - th2(:,:,:) = th2(:,:,:) + uy2(:,:,:) * ty2(:,:,:) - by2(:,:,:) * te2(:,:,:) - ti2(:,:,:) = ti2(:,:,:) + uy2(:,:,:) * tz2(:,:,:) - by2(:,:,:) * tf2(:,:,:) call transpose_y_to_z(ux2,ux3) call transpose_y_to_z(uy2,uy3) @@ -782,65 +769,43 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) !WORK Z-PENCILS - td3(:,:,:) = uz3(:,:,:)*bx3(:,:,:) - bz3(:,:,:)*ux3(:,:,:) - te3(:,:,:) = uz3(:,:,:)*by3(:,:,:) - bz3(:,:,:)*uy3(:,:,:) - tf3(:,:,:) = uz3(:,:,:)*bz3(:,:,:) - bz3(:,:,:)*uz3(:,:,:) - - - call derz (tg3,td3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcx*ubcz) - call derz (th3,te3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcy*ubcz) - call derz (ti3,tf3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,ubcz*ubcz) - - call derz (td3,ux3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,ubcx) - call derz (te3,uy3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,ubcy) - call derz (tf3,uz3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcz) - - call derz (tx3,bx3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,ubcx) - call derz (ty3,by3,di3,sz,ffzp,fszp,fwzp,zsize(1),zsize(2),zsize(3),1,ubcy) - call derz (tz3,bz3,di3,sz,ffz,fsz,fwz,zsize(1),zsize(2),zsize(3),0,ubcz) + ta3(:,:,:) = uy3(:,:,:)*bz3(:,:,:) - by3(:,:,:)*uz3(:,:,:) + tb3(:,:,:) = ux3(:,:,:)*bz3(:,:,:) - bx3(:,:,:)*uz3(:,:,:) + tc3(:,:,:) = ux3(:,:,:)*by3(:,:,:) - bx3(:,:,:)*uy3(:,:,:) - ! Convective terms of z-pencil in ta3,tb3,tc3 - ta3(:,:,:) = tg3(:,:,:) + uz3(:,:,:) * tx3(:,:,:) - bz3(:,:,:) * td3(:,:,:) - tb3(:,:,:) = th3(:,:,:) + uz3(:,:,:) * ty3(:,:,:) - bz3(:,:,:) * te3(:,:,:) - tc3(:,:,:) = ti3(:,:,:) + uz3(:,:,:) * tz3(:,:,:) - bz3(:,:,:) * tf3(:,:,:) - - - ! Convective terms of z-pencil are in ta3 -> td3, tb3 -> te3, tc3 -> tf3 - td3(:,:,:) = ta3(:,:,:) - te3(:,:,:) = tb3(:,:,:) - tf3(:,:,:) = tc3(:,:,:) + call derzBx (td3,ta3,di3,sz,ffzB(:,1),fszB(:,1),fwzB(:,1),zsize(1),zsize(2),zsize(3),0,ubcx*ubcz) + call derzBy (te3,tb3,di3,sz,ffzB(:,2),fszB(:,2),fwzB(:,2),zsize(1),zsize(2),zsize(3),0,ubcy*ubcz) + call derzBz (tf3,tc3,di3,sz,ffzpB(:,3),fszpB(:,3),fwzpB(:,3),zsize(1),zsize(2),zsize(3),1,ubcz*ubcz) !DIFFUSIVE TERMS IN Z - call derzz (ta3,bx3,di3,sz,sfzp,sszp,swzp,zsize(1),zsize(2),zsize(3),1,ubcx) - call derzz (tb3,by3,di3,sz,sfzp,sszp,swzp,zsize(1),zsize(2),zsize(3),1,ubcy) - call derzz (tc3,bz3,di3,sz,sfz ,ssz ,swz ,zsize(1),zsize(2),zsize(3),0,ubcz) - + call derzzBx (ta3,bx3,di3,sz,sfzpB(:,1),sszpB(:,1),swzpB(:,1),zsize(1),zsize(2),zsize(3),1,ubcx) + call derzzBy (tb3,by3,di3,sz,sfzpB(:,2),sszpB(:,2),swzpB(:,2),zsize(1),zsize(2),zsize(3),1,ubcy) + call derzzBz (tc3,bz3,di3,sz,sfzB(:,3) ,sszB(:,3) ,swzB(:,3) ,zsize(1),zsize(2),zsize(3),0,ubcz) ! Add convective and diffusive terms of z-pencil (half for skew-symmetric) - td3(:,:,:) = rrem*ta3(:,:,:) - half * td3(:,:,:) - te3(:,:,:) = rrem*tb3(:,:,:) - half * te3(:,:,:) - tf3(:,:,:) = rrem*tc3(:,:,:) - half * tf3(:,:,:) + td3(:,:,:) = rrem*ta3(:,:,:) - te3(:,:,:) + te3(:,:,:) = rrem*tb3(:,:,:) + te3(:,:,:) + tf3(:,:,:) = rrem*tc3(:,:,:) !WORK Y-PENCILS - call transpose_z_to_y(td3,td2) - call transpose_z_to_y(te3,te2) - call transpose_z_to_y(tf3,tf2) + call transpose_z_to_y(td3,ta2) + call transpose_z_to_y(te3,tb2) + call transpose_z_to_y(tf3,tc2) - ! Convective terms of y-pencil (tg2,th2,ti2) and sum of convective and diffusive terms of z-pencil (td2,te2,tf2) are now in tg2, th2, ti2 (half for skew-symmetric) - tg2(:,:,:) = td2(:,:,:) - half * tg2(:,:,:) - th2(:,:,:) = te2(:,:,:) - half * th2(:,:,:) - ti2(:,:,:) = tf2(:,:,:) - half * ti2(:,:,:) + ta2(:,:,:) = ta2(:,:,:) + tf2(:,:,:) + tb2(:,:,:) = tb2(:,:,:) + tc2(:,:,:) = tc2(:,:,:) - td2(:,:,:) !DIFFUSIVE TERMS IN Y if (iimplicit.le.0) then !-->for ux - call deryy (td2,bx2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1,ubcx) + call deryyBx (td2,bx2,di2,sy,sfypB(:,1),ssypB(:,1),swypB(:,1),ysize(1),ysize(2),ysize(3),1,ubcx) if (istret.ne.0) then - call dery (te2,bx2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) + call deryBx (te2,bx2,di2,sy,ffypB(:,1),fsypB(:,1),fwypB(:,1),ppy,ysize(1),ysize(2),ysize(3),1,ubcx) do k = 1,ysize(3) do j = 1,ysize(2) do i = 1,ysize(1) @@ -851,9 +816,9 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) endif !-->for uy - call deryy (te2,by2,di2,sy,sfy,ssy,swy,ysize(1),ysize(2),ysize(3),0,ubcy) + call deryyBy (te2,by2,di2,sy,sfyB(:,2),ssyB(:,2),swyB(:,2),ysize(1),ysize(2),ysize(3),0,ubcy) if (istret.ne.0) then - call dery (tf2,by2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy) + call deryBy (tf2,by2,di2,sy,ffyB(:,2),fsyB(:,2),fwyB(:,2),ppy,ysize(1),ysize(2),ysize(3),0,ubcy) do k = 1,ysize(3) do j = 1,ysize(2) do i = 1,ysize(1) @@ -864,9 +829,9 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) endif !-->for uz - call deryy (tf2,bz2,di2,sy,sfyp,ssyp,swyp,ysize(1),ysize(2),ysize(3),1,ubcz) + call deryyBz (tf2,bz2,di2,sy,sfypB(:,3),ssypB(:,3),swypB(:,3),ysize(1),ysize(2),ysize(3),1,ubcz) if (istret.ne.0) then - call dery (tj2,bz2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcz) + call deryBz (tj2,bz2,di2,sy,ffypB(:,3),fsypB(:,3),fwypB(:,3),ppy,ysize(1),ysize(2),ysize(3),1,ubcz) do k = 1,ysize(3) do j = 1,ysize(2) do i = 1,ysize(1) @@ -879,7 +844,7 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) if (istret.ne.0) then !-->for ux - call dery (te2,bx2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcx) + call deryBx (te2,bx2,di2,sy,ffypB(:,1),fsypB(:,1),fwypB(:,1),ppy,ysize(1),ysize(2),ysize(3),1,ubcx) do k=1,ysize(3) do j=1,ysize(2) do i=1,ysize(1) @@ -888,7 +853,7 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) enddo enddo !-->for uy - call dery (tf2,by2,di2,sy,ffy,fsy,fwy,ppy,ysize(1),ysize(2),ysize(3),0,ubcy) + call deryBy (tf2,by2,di2,sy,ffyB(:,2),fsyB(:,2),fwyB(:,2),ppy,ysize(1),ysize(2),ysize(3),0,ubcy) do k=1,ysize(3) do j=1,ysize(2) do i=1,ysize(1) @@ -897,7 +862,7 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) enddo enddo !-->for uz - call dery (tj2,bz2,di2,sy,ffyp,fsyp,fwyp,ppy,ysize(1),ysize(2),ysize(3),1,ubcz) + call deryBz (tj2,bz2,di2,sy,ffypB(:,3),fsypB(:,3),fwypB(:,3),ppy,ysize(1),ysize(2),ysize(3),1,ubcz) do k=1,ysize(3) do j=1,ysize(2) do i=1,ysize(1) @@ -916,9 +881,9 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) endif ! Add diffusive terms of y-pencil to convective and diffusive terms of y- and z-pencil - ta2(:,:,:) = rrem*td2(:,:,:) + tg2(:,:,:) - tb2(:,:,:) = rrem*te2(:,:,:) + th2(:,:,:) - tc2(:,:,:) = rrem*tf2(:,:,:) + ti2(:,:,:) + ta2(:,:,:) = rrem*td2(:,:,:) + ta2(:,:,:) + tb2(:,:,:) = rrem*te2(:,:,:) + tb2(:,:,:) + tc2(:,:,:) = rrem*tf2(:,:,:) + tc2(:,:,:) !WORK X-PENCILS call transpose_y_to_x(ta2,ta1) @@ -926,18 +891,14 @@ subroutine mhd_rhs_eq(dB,B,ux1,uy1,uz1) call transpose_y_to_x(tc2,tc1) !diff+conv. terms !DIFFUSIVE TERMS IN X - call derxx (td1,B(:,:,:,1),di1,sx,sfx ,ssx ,swx ,xsize(1),xsize(2),xsize(3),0,ubcx) - call derxx (te1,B(:,:,:,2),di1,sx,sfxp,ssxp,swxp,xsize(1),xsize(2),xsize(3),1,ubcy) - call derxx (tf1,B(:,:,:,3),di1,sx,sfxp,ssxp,swxp,xsize(1),xsize(2),xsize(3),1,ubcz) - - td1(:,:,:) = rrem * td1(:,:,:) - te1(:,:,:) = rrem * te1(:,:,:) - tf1(:,:,:) = rrem * tf1(:,:,:) + call derxxBx (td1,B(:,:,:,1),di1,sx,sfxB(:,1) ,ssxB(:,1) ,swxB(:,1) ,xsize(1),xsize(2),xsize(3),0,ubcx) + call derxxBy (te1,B(:,:,:,2),di1,sx,sfxpB(:,2),ssxpB(:,2),swxpB(:,2),xsize(1),xsize(2),xsize(3),1,ubcy) + call derxxBz (tf1,B(:,:,:,3),di1,sx,sfxpB(:,3),ssxpB(:,3),swxpB(:,3),xsize(1),xsize(2),xsize(3),1,ubcz) !FINAL SUM: DIFF TERMS + CONV TERMS - dB(:,:,:,1) = ta1(:,:,:) - half*tg1(:,:,:) + td1(:,:,:) - dB(:,:,:,2) = tb1(:,:,:) - half*th1(:,:,:) + te1(:,:,:) - dB(:,:,:,3) = tc1(:,:,:) - half*ti1(:,:,:) + tf1(:,:,:) + dB(:,:,:,1) = rrem * td1(:,:,:) + ta1(:,:,:) + dB(:,:,:,2) = rrem * te1(:,:,:) + tb1(:,:,:) -tf1(:,:,:) + dB(:,:,:,3) = rrem * te1(:,:,:) + tc1(:,:,:) +te1(:,:,:) return @@ -956,30 +917,16 @@ subroutine solve_poisson_mhd real(mytype),dimension(xsize(1),xsize(2),xsize(3),1:3) :: dphib integer :: i,j,k,nlock,poissiter - logical :: converged ! nlock=1 !! Corresponds to computing div(u*) ! - phib=divergence_scalar(Bm,nlock) - ! - converged=.false. - ! - poissiter = 0 - ! - do while(.not.converged) - ! + do poissiter = 1, 1 + phib=divergence_scalar(Bm,nlock) !todo: this will have incorrect BCs? call poisson(phib) - ! CALL gradp(dphib(:,:,:,1),dphib(:,:,:,2),dphib(:,:,:,3),phib) - ! - converged=.true. - ! - poissiter = poissiter + 1 - ! + Bm=Bm-dphib enddo ! - Bm=Bm-dphib - ! sync_Bm_needed=.true. ! end subroutine solve_poisson_mhd @@ -1264,7 +1211,6 @@ function divergence_scalar(vec,nlock) result(pp3) !WORK Y-PENCILS call interyvp(upi2,duxdxp2,dipp2,sy,cifyp6,cisyp6,ciwyp6,(ph1%yen(1)-ph1%yst(1)+1),ysize(2),nymsize,ysize(3),1) call deryvp(duydypi2,uyp2,dipp2,sy,cfy6,csy6,cwy6,ppyi,(ph1%yen(1)-ph1%yst(1)+1),ysize(2),nymsize,ysize(3),0) - !! Compute sum dudx + dvdy duydypi2(:,:,:) = duydypi2(:,:,:) + upi2(:,:,:) diff --git a/src/module_param.f90 b/src/module_param.f90 index 1ea03d2f6..2a0473643 100644 --- a/src/module_param.f90 +++ b/src/module_param.f90 @@ -70,6 +70,14 @@ module variables real(mytype),allocatable,dimension(:) :: ffzS,sfzS,fszS,fwzS,sszS,swzS real(mytype),allocatable,dimension(:) :: ffzpS,sfzpS,fszpS,fwzpS,sszpS,swzpS + real(mytype),allocatable,dimension(:,:) :: ffxB,sfxB,fsxB,fwxB,ssxB,swxB + real(mytype),allocatable,dimension(:,:) :: ffxpB,sfxpB,fsxpB,fwxpB,ssxpB,swxpB + real(mytype),allocatable,dimension(:,:) :: ffyB,sfyB,fsyB,fwyB,ssyB,swyB + real(mytype),allocatable,dimension(:,:) :: ffypB,sfypB,fsypB,fwypB,ssypB,swypB + real(mytype),allocatable,dimension(:,:) :: ffzB,sfzB,fszB,fwzB,sszB,swzB + real(mytype),allocatable,dimension(:,:) :: ffzpB,sfzpB,fszpB,fwzpB,sszpB,swzpB + + real(mytype), save, allocatable, dimension(:,:) :: sx,vx real(mytype), save, allocatable, dimension(:,:) :: sy,vy real(mytype), save, allocatable, dimension(:,:) :: sz,vz @@ -127,7 +135,10 @@ module variables ! scalar, ncly1 = 2, nclyn = 1, npaire = 1 real(mytype), allocatable, target, dimension(:,:) :: aam211t,bbm211t,ccm211t,ddm211t,eem211t,ggm211t,hhm211t,wwm211t,zzm211t real(mytype), allocatable, target, dimension(:,:) :: rrm211t,qqm211t,vvm211t,ssm211t - + ! Bx, ncly1 = 2, nclyn = 2 + real(mytype), allocatable, target, dimension(:,:) :: aamB,bbmB,ccmB,ddmB,eemB,ggmB,hhmB,wwmB,zzmB + real(mytype), allocatable, target, dimension(:,:) :: rrmB,qqmB,vvmB,ssmB + real(mytype), allocatable, target, dimension(:,:) :: sssmB, zzzmB, ttmB, uumB !!Nona ABSTRACT INTERFACE SUBROUTINE DERIVATIVE_X(t,u,r,s,ff,fs,fw,nx,ny,nz,npaire,lind) @@ -176,6 +187,22 @@ END SUBROUTINE DERIVATIVE_Z derzz_00,derzz_11,derzz_12,derzz_21,derzz_22 PROCEDURE (DERIVATIVE_Z), POINTER :: derz,derzz,derzS,derzzS + procedure (DERIVATIVE_X), pointer :: derxBx, derxxBx + procedure (DERIVATIVE_Y), pointer :: deryBx + procedure (DERIVATIVE_YY), pointer :: deryyBx + procedure (DERIVATIVE_Z), pointer :: derzBx, derzzBx + + procedure (DERIVATIVE_X), pointer :: derxBy, derxxBy + procedure (DERIVATIVE_Y), pointer :: deryBy + procedure (DERIVATIVE_YY), pointer :: deryyBy + procedure (DERIVATIVE_Z), pointer :: derzBy, derzzBy + + procedure (DERIVATIVE_X), pointer :: derxBz, derxxBz + procedure (DERIVATIVE_Y), pointer :: deryBz + procedure (DERIVATIVE_YY), pointer :: deryyBz + procedure (DERIVATIVE_Z), pointer :: derzBz, derzzBz + + !O6SVV real(mytype),allocatable,dimension(:) :: newsm,newtm,newsmt,newtmt !real(mytype),allocatable,dimension(:) :: newrm,ttm,newrmt,ttmt @@ -266,6 +293,10 @@ module param integer :: nclx1,nclxn,ncly1,nclyn,nclz1,nclzn integer :: nclxS1,nclxSn,nclyS1,nclySn,nclzS1,nclzSn + integer :: nclxBx1,nclxBxn,nclyBx1,nclyBxn,nclzBx1,nclzBxn + integer :: nclxBy1,nclxByn,nclyBy1,nclyByn,nclzBy1,nclzByn + integer :: nclxBz1,nclxBzn,nclyBz1,nclyBzn,nclzBz1,nclzBzn + integer, dimension(3) :: nclxB1,nclxBn,nclyB1,nclyBn,nclzB1,nclzBn !logical variable for boundary condition that is true in periodic case !and false otherwise @@ -322,7 +353,7 @@ module param real(mytype) :: cfl_diff_x,cfl_diff_y,cfl_diff_z,cfl_diff_sum !! - real(mytype) :: xcst + real(mytype) :: xcst, xcstB real(mytype), allocatable, dimension(:) :: xcst_sc real(mytype), allocatable, dimension(:,:) :: alpha_sc, beta_sc, g_sc real(mytype) :: g_bl_inf, f_bl_inf @@ -376,7 +407,7 @@ module param integer :: NTurbines, NActuatorlines character, dimension(100) :: TurbinesPath*80, ActuatorlinesPath*80 real(mytype) :: eps_factor ! Smoothing factor - + character :: filesauve*80, filenoise*80, & nchamp*80,filepath*80, fileturb*80, filevisu*80, datapath*80 real(mytype), dimension(5) :: adt,bdt,cdt,ddt,gdt diff --git a/src/parameters.f90 b/src/parameters.f90 index c55f9db48..aa00e02d6 100644 --- a/src/parameters.f90 +++ b/src/parameters.f90 @@ -15,7 +15,7 @@ subroutine parameter(input_i3d) use mpi - + use iso_fortran_env use param @@ -77,7 +77,11 @@ subroutine parameter(input_i3d) NAMELIST/ALMParam/iturboutput,NTurbines,TurbinesPath,NActuatorlines,ActuatorlinesPath,eps_factor,rho_air NAMELIST/ADMParam/Ndiscs,ADMcoords,C_T,aind,iturboutput,rho_air NAMELIST/PartiParam/ipartiout,numpartix,partirange,lfluidforce,lorentzforce - NAMELIST/MHDParam/mhd_active,mhd_equation,hartmann,stuart,rem + NAMELIST/MHDParam/mhd_active,mhd_equation,hartmann,stuart,rem, & + nclxBx1, nclxBxn, nclyBx1, nclyBxn, nclzBx1, nclzBxn, & + nclxBy1, nclxByn, nclyBy1, nclyByn, nclzBy1, nclzByn, & + nclxBz1, nclxBzn, nclyBz1, nclyBzn, nclzBz1, nclzBzn + #ifdef DEBG if (nrank == 0) write(*,*) '# parameter start' @@ -126,12 +130,12 @@ subroutine parameter(input_i3d) allocate(xld(nvol), xrd(nvol), yld(nvol), yud(nvol))!, zld(nvol), zrd(nvol)) read(10, nml=ForceCVs); rewind(10) endif - + !! Set Scalar BCs same as fluid (may be overridden) [DEFAULT] nclxS1 = nclx1; nclxSn = nclxn nclyS1 = ncly1; nclySn = nclyn nclzS1 = nclz1; nclzSn = nclzn - + if (numscalar.ne.0) then iscalar = 1 @@ -205,7 +209,30 @@ subroutine parameter(input_i3d) endif read(10, nml=PartiParam); rewind(10) !! read particle - read(10, nml=MHDParam); rewind(10) !! read mhd + read(10, nml=MHDParam); rewind(10) !! read mhd + nclxB1(1) = nclxBx1 + nclxB1(2) = nclxBy1 + nclxB1(3) = nclxBz1 + nclxBn(1) = nclxBxn + nclxBn(2) = nclxByn + nclxBn(3) = nclxBzn + nclyB1(1) = nclyBx1 + nclyB1(2) = nclyBy1 + nclyB1(3) = nclyBz1 + nclyBn(1) = nclyBxn + nclyBn(2) = nclyByn + nclyBn(3) = nclyBzn + nclzB1(1) = nclzBx1 + nclzB1(2) = nclzBy1 + nclzB1(3) = nclzBz1 + nclzBn(1) = nclzBxn + nclzBn(2) = nclzByn + nclzBn(3) = nclzBzn + + + + + ! !! These are the 'optional'/model parameters ! read(10, nml=ScalarParam) if(ilesmod==0) then @@ -269,10 +296,10 @@ subroutine parameter(input_i3d) xnu=one/re !! Constant pressure gradient, re = Re_tau -> use to compute Re_centerline if (cpg) then - re_cent = (re/0.116_mytype)**(1.0_mytype/0.88_mytype) - xnu = one/re_cent ! viscosity based on Re_cent to keep same scaling as CFR - ! - fcpg = two/yly * (re/re_cent)**2 + re_cent = (re/0.116_mytype)**(1.0_mytype/0.88_mytype) + xnu = one/re_cent ! viscosity based on Re_cent to keep same scaling as CFR + ! + fcpg = two/yly * (re/re_cent)**2 end if if (ilmn) then @@ -317,6 +344,7 @@ subroutine parameter(input_i3d) stop endif if (iscalar.eq.1) xcst_sc = xcst / sc + if (mhd_active) xcstB = xcst / rem * re endif if (itype==itype_tbl.and.A_tr .gt. zero.and.nrank==0) write(*,*) "TBL tripping is active" @@ -368,17 +396,17 @@ subroutine parameter(input_i3d) endif print *,'===========================================================' if (itype.eq.itype_channel) then - if (.not.cpg) then - write(*,*) 'Channel forcing with constant flow rate (CFR)' - write(*,"(' Re_cl : ',F17.3)") re - else - write(*,*) 'Channel forcing with constant pressure gradient (CPG)' - write(*,"(' Re_tau : ',F17.3)") re - write(*,"(' Re_cl (estimated) : ',F17.3)") re_cent - write(*,"(' fcpg : ',F17.8)") fcpg - end if + if (.not.cpg) then + write(*,*) 'Channel forcing with constant flow rate (CFR)' + write(*,"(' Re_cl : ',F17.3)") re + else + write(*,*) 'Channel forcing with constant pressure gradient (CPG)' + write(*,"(' Re_tau : ',F17.3)") re + write(*,"(' Re_cl (estimated) : ',F17.3)") re_cent + write(*,"(' fcpg : ',F17.8)") fcpg + end if else - write(*,"(' Reynolds number Re : ',F17.3)") re + write(*,"(' Reynolds number Re : ',F17.3)") re endif write(*,"(' xnu : ',F17.8)") xnu write(*,*) '===========================================================' @@ -387,53 +415,53 @@ subroutine parameter(input_i3d) write(*,"(' Time step dt : ',F17.8)") dt ! if (itimescheme.eq.1) then - !print *,'Temporal scheme : Forwards Euler' - write(*,"(' Temporal scheme : ',A20)") "Forwards Euler" + !print *,'Temporal scheme : Forwards Euler' + write(*,"(' Temporal scheme : ',A20)") "Forwards Euler" elseif (itimescheme.eq.2) then - !print *,'Temporal scheme : Adams-bashforth 2' - write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 2" + !print *,'Temporal scheme : Adams-bashforth 2' + write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 2" elseif (itimescheme.eq.3) then - !print *,'Temporal scheme : Adams-bashforth 3' - write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 3" + !print *,'Temporal scheme : Adams-bashforth 3' + write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 3" elseif (itimescheme.eq.4) then - !print *,'Temporal scheme : Adams-bashforth 4' - write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 4" - print *,'Error: Adams-bashforth 4 not implemented!' - stop + !print *,'Temporal scheme : Adams-bashforth 4' + write(*,"(' Temporal scheme : ',A20)") "Adams-bashforth 4" + print *,'Error: Adams-bashforth 4 not implemented!' + stop elseif (itimescheme.eq.5) then - !print *,'Temporal scheme : Runge-kutta 3' - write(*,"(' Temporal scheme : ',A20)") "Runge-kutta 3" + !print *,'Temporal scheme : Runge-kutta 3' + write(*,"(' Temporal scheme : ',A20)") "Runge-kutta 3" elseif (itimescheme.eq.6) then - !print *,'Temporal scheme : Runge-kutta 4' - write(*,"(' Temporal scheme : ',A20)") "Runge-kutta 4" - print *,'Error: Runge-kutta 4 not implemented!' - stop + !print *,'Temporal scheme : Runge-kutta 4' + write(*,"(' Temporal scheme : ',A20)") "Runge-kutta 4" + print *,'Error: Runge-kutta 4 not implemented!' + stop else - print *,'Error: itimescheme must be specified as 1-6' - stop + print *,'Error: itimescheme must be specified as 1-6' + stop endif ! if (iimplicit.ne.0) then - if (iimplicit.eq.1) then - write(*,"(' ',A40)") "With backward Euler for Y diffusion" - else if (iimplicit.eq.2) then - write(*,"(' ',A40)") "With CN for Y diffusion" - endif + if (iimplicit.eq.1) then + write(*,"(' ',A40)") "With backward Euler for Y diffusion" + else if (iimplicit.eq.2) then + write(*,"(' ',A40)") "With CN for Y diffusion" + endif endif ! if (ilesmod.ne.0) then - write(*,*) ' : DNS' + write(*,*) ' : DNS' else - if (jles==1) then - write(*,*) ' : Phys Smag' - else if (jles==2) then - write(*,*) ' : Phys WALE' - else if (jles==3) then - write(*,*) ' : Phys dyn. Smag' - else if (jles==4) then - write(*,*) ' : iSVV' - else - endif + if (jles==1) then + write(*,*) ' : Phys Smag' + else if (jles==2) then + write(*,*) ' : Phys WALE' + else if (jles==3) then + write(*,*) ' : Phys dyn. Smag' + else if (jles==4) then + write(*,*) ' : iSVV' + else + endif endif write(*,*) '===========================================================' write(*,"(' ifirst : ',I17)") ifirst @@ -456,43 +484,43 @@ subroutine parameter(input_i3d) if (iscalar==1) write(*,"(' Scalar : ',A17)") "on" write(*,"(' numscalar : ',I17)") numscalar if (iscalar.eq.1) then - do is=1, numscalar - write(*,"(' Schmidt number sc(',I2,') : ',F17.8)") is, sc(is) - write(*,"(' Richardson n. ri(',I2,') : ',F17.8)") is, ri(is) - if (scalar_lbound(is).gt.-huge(one)) then - write(*,"(' Lower bound (',I2,') : ',F17.8)") is, scalar_lbound(is) - else - ! This is the default option, no information printed in the listing - endif - if (scalar_ubound(is).lt.huge(one)) then - write(*,"(' Upper bound (',I2,') : ',F17.8)") is, scalar_ubound(is) - else - ! This is the default option, no information printed in the listing - endif - if (iscalar.eq.1) then - if (nclxS1.eq.1 .or. nclxSn.eq.1 .or. & - nclyS1.eq.1 .or. nclySn.eq.1 .or. & - nclzS1.eq.1 .or. nclzSn.eq.1) then - if (sc_even(is)) then - ! This is the default option, no information printed in the listing - else - write(*,"(' Scalar ',I2,' is odd')") is - endif - endif - if (sc_skew(is)) then - write(*,"(' Scalar ',I2,' with skew-symmetric convection')") is - else - ! This is the default option, no information printed in the listing - endif - endif - end do + do is=1, numscalar + write(*,"(' Schmidt number sc(',I2,') : ',F17.8)") is, sc(is) + write(*,"(' Richardson n. ri(',I2,') : ',F17.8)") is, ri(is) + if (scalar_lbound(is).gt.-huge(one)) then + write(*,"(' Lower bound (',I2,') : ',F17.8)") is, scalar_lbound(is) + else + ! This is the default option, no information printed in the listing + endif + if (scalar_ubound(is).lt.huge(one)) then + write(*,"(' Upper bound (',I2,') : ',F17.8)") is, scalar_ubound(is) + else + ! This is the default option, no information printed in the listing + endif + if (iscalar.eq.1) then + if (nclxS1.eq.1 .or. nclxSn.eq.1 .or. & + nclyS1.eq.1 .or. nclySn.eq.1 .or. & + nclzS1.eq.1 .or. nclzSn.eq.1) then + if (sc_even(is)) then + ! This is the default option, no information printed in the listing + else + write(*,"(' Scalar ',I2,' is odd')") is + endif + endif + if (sc_skew(is)) then + write(*,"(' Scalar ',I2,' with skew-symmetric convection')") is + else + ! This is the default option, no information printed in the listing + endif + endif + end do endif write(*,*) '===========================================================' if (lpartack) then - write(*,"(' Particle tracking : ',A17)") "on" - write(*,"(' numparticle : ',I17)") numparticle + write(*,"(' Particle tracking : ',A17)") "on" + write(*,"(' numparticle : ',I17)") numparticle else - write(*,"(' Particle tracking : ',A17)") "off" + write(*,"(' Particle tracking : ',A17)") "off" endif write(*,*) '===========================================================' write(*,"(' spinup_time : ',I17)") spinup_time @@ -500,17 +528,17 @@ subroutine parameter(input_i3d) write(*,*) '===========================================================' if (iibm==0) write(*,"(' Immersed boundary : ',A17)") "off" if (iibm.gt.1) then - write(*,"(' Immersed boundary : ',A17)") "on" - write(*,"(' iibm : ',I17)") iibm + write(*,"(' Immersed boundary : ',A17)") "on" + write(*,"(' iibm : ',I17)") iibm end if if (iibm==1) write(*,*) 'Simple immersed boundary method' if (iibm==2) then - write(*,*) 'Lagrangian polynomial reconstruction' - write(*,*) '===========================================================' - write(*,"(' npif : ',I17)") npif - write(*,"(' izap : ',I17)") izap - write(*,"(' nraf : ',I17)") nraf - write(*,"(' nobjmax : ',I17)") nobjmax + write(*,*) 'Lagrangian polynomial reconstruction' + write(*,*) '===========================================================' + write(*,"(' npif : ',I17)") npif + write(*,"(' izap : ',I17)") izap + write(*,"(' nraf : ',I17)") nraf + write(*,"(' nobjmax : ',I17)") nobjmax end if write(*,*) '===========================================================' write(*,"(' Boundary condition velocity field: ')") @@ -519,11 +547,11 @@ subroutine parameter(input_i3d) write(*,"(' nclz1, nclzn : ',I15,',',I1 )") nclz1,nclzn write(*,*) '===========================================================' if ((iscalar==1).or.(ilmn)) then - write(*,"(' Boundary condition scalar field: ')") - write(*,"(' nclxS1, nclxSn : ',I15,',',I1 )") nclxS1,nclxSn - write(*,"(' nclyS1, nclySn : ',I15,',',I1 )") nclyS1,nclySn - write(*,"(' nclzS1, nclzSn : ',I15,',',I1 )") nclzS1,nclzSn - write(*,*) '===========================================================' + write(*,"(' Boundary condition scalar field: ')") + write(*,"(' nclxS1, nclxSn : ',I15,',',I1 )") nclxS1,nclxSn + write(*,"(' nclyS1, nclySn : ',I15,',',I1 )") nclyS1,nclySn + write(*,"(' nclzS1, nclzSn : ',I15,',',I1 )") nclzS1,nclzSn + write(*,*) '===========================================================' endif #ifdef DOUBLE_PREC @@ -561,7 +589,7 @@ subroutine parameter(input_i3d) endif write(*,*) '===========================================================' endif - + if (iibm.eq.3) then ! This is only for the Cubic Spline Reconstruction npif=npif+1 endif diff --git a/src/schemes.f90 b/src/schemes.f90 index 8bfbca879..9e91aa4d6 100644 --- a/src/schemes.f90 +++ b/src/schemes.f90 @@ -15,7 +15,7 @@ subroutine schemes() USE variables USE var USE ydiff_implicit, only : init_implicit, implicit_schemes - use mhd, only: mhd_active + use mhd, only: mhd_active, mhd_equation implicit none @@ -164,6 +164,209 @@ subroutine schemes() alsakz,askz,bskz,cskz,dskz,& sfzS,sszS,swzS,sfzpS,sszpS,swzpS,dz2,nz,nclzS1,nclzSn) endif + + if( mhd_active .and. mhd_equation ) then + ! First derivative + if (nclxBx1.eq.0.and.nclxBxn.eq.0) derxBx => derx_00 + if (nclxBx1.eq.1.and.nclxBxn.eq.1) derxBx => derx_11 + if (nclxBx1.eq.1.and.nclxBxn.eq.2) derxBx => derx_12 + if (nclxBx1.eq.2.and.nclxBxn.eq.1) derxBx => derx_21 + if (nclxBx1.eq.2.and.nclxBxn.eq.2) derxBx => derx_22 + ! + if (nclyBx1.eq.0.and.nclyBxn.eq.0) deryBx => dery_00 + if (nclyBx1.eq.1.and.nclyBxn.eq.1) deryBx => dery_11 + if (nclyBx1.eq.1.and.nclyBxn.eq.2) deryBx => dery_12 + if (nclyBx1.eq.2.and.nclyBxn.eq.1) deryBx => dery_21 + if (nclyBx1.eq.2.and.nclyBxn.eq.2) deryBx => dery_22 + ! + if (nclzBx1.eq.0.and.nclzBxn.eq.0) derzBx => derz_00 + if (nclzBx1.eq.1.and.nclzBxn.eq.1) derzBx => derz_11 + if (nclzBx1.eq.1.and.nclzBxn.eq.2) derzBx => derz_12 + if (nclzBx1.eq.2.and.nclzBxn.eq.1) derzBx => derz_21 + if (nclzBx1.eq.2.and.nclzBxn.eq.2) derzBx => derz_22 + ! SecondBxderivative + if (nclxBx1.eq.0.and.nclxBxn.eq.0) derxxBx => derxx_00 + if (nclxBx1.eq.1.and.nclxBxn.eq.1) derxxBx => derxx_11 + if (nclxBx1.eq.1.and.nclxBxn.eq.2) derxxBx => derxx_12 + if (nclxBx1.eq.2.and.nclxBxn.eq.1) derxxBx => derxx_21 + if (nclxBx1.eq.2.and.nclxBxn.eq.2) derxxBx => derxx_22 + !y + if (nclyBx1.eq.0.and.nclyBxn.eq.0) deryyBx => deryy_00 + if (nclyBx1.eq.1.and.nclyBxn.eq.1) deryyBx => deryy_11 + if (nclyBx1.eq.1.and.nclyBxn.eq.2) deryyBx => deryy_12 + if (nclyBx1.eq.2.and.nclyBxn.eq.1) deryyBx => deryy_21 + if (nclyBx1.eq.2.and.nclyBxn.eq.2) deryyBx => deryy_22 + !z + if (nclzBx1.eq.0.and.nclzBxn.eq.0) derzzBx => derzz_00 + if (nclzBx1.eq.1.and.nclzBxn.eq.1) derzzBx => derzz_11 + if (nclzBx1.eq.1.and.nclzBxn.eq.2) derzzBx => derzz_12 + if (nclzBx1.eq.2.and.nclzBxn.eq.1) derzzBx => derzz_21 + if (nclzBx1.eq.2.and.nclzBxn.eq.2) derzzBx => derzz_22 + call first_derivative(alfa1x,af1x,bf1x,cf1x,df1x,alfa2x,af2x,alfanx,afnx,bfnx,& + cfnx,dfnx,alfamx,afmx,alfaix,afix,bfix,& + ffxB(:,1),fsxB(:,1),fwxB(:,1),ffxpB(:,1),fsxpB(:,1),fwxpB(:,1),dx,nx,nclxBx1,nclxBxn) + call first_derivative(alfa1y,af1y,bf1y,cf1y,df1y,alfa2y,af2y,alfany,afny,bfny,& + cfny,dfny,alfamy,afmy,alfajy,afjy,bfjy,& + ffyB(:,1),fsyB(:,1),fwyB(:,1),ffypB(:,1),fsypB(:,1),fwypB(:,1),dy,ny,nclyBx1,nclyBxn) + call first_derivative(alfa1z,af1z,bf1z,cf1z,df1z,alfa2z,af2z,alfanz,afnz,bfnz,& + cfnz,dfnz,alfamz,afmz,alfakz,afkz,bfkz,& + ffzB(:,1),fszB(:,1),fwzB(:,1),ffzpB(:,1),fszpB(:,1),fwzpB(:,1),dz,nz,nclzBx1,nclzBxn) + call second_derivative(alsa1x,as1x,bs1x,& + cs1x,ds1x,alsa2x,as2x,alsanx,asnx,bsnx,csnx,dsnx,alsamx,& + asmx,alsa3x,as3x,bs3x,alsatx,astx,bstx,& + alsa4x,as4x,bs4x,cs4x,& + alsattx,asttx,bsttx,csttx,& + alsaix,asix,bsix,csix,dsix,& + sfxB(:,1),ssxB(:,1),swxB(:,1),sfxpB(:,1),ssxpB(:,1),swxpB(:,1),dx2,nx,nclxBx1,nclxBxn) + call second_derivative(alsa1y,as1y,bs1y,& + cs1y,ds1y,alsa2y,as2y,alsany,asny,bsny,csny,dsny,alsamy,& + asmy,alsa3y,as3y,bs3y,alsaty,asty,bsty,& + alsa4y,as4y,bs4y,cs4y,& + alsatty,astty,bstty,cstty,& + alsajy,asjy,bsjy,csjy,dsjy,& + sfyB(:,1),ssyB(:,1),swyB(:,1),sfypB(:,1),ssypB(:,1),swypB(:,1),dy2,ny,nclyBx1,nclyBxn) + call second_derivative(alsa1z,as1z,bs1z,& + cs1z,ds1z,alsa2z,as2z,alsanz,asnz,bsnz,csnz,dsnz,alsamz,& + asmz,alsa3z,as3z,bs3z,alsatz,astz,bstz,& + alsa4z,as4z,bs4z,cs4z,& + alsattz,asttz,bsttz,csttz,& + alsakz,askz,bskz,cskz,dskz,& + sfzB(:,1),sszB(:,1),swzB(:,1),sfzpB(:,1),sszpB(:,1),swzpB(:,1),dz2,nz,nclzBx1,nclzBxn) + + ! First derivative + if (nclxBy1.eq.0.and.nclxByn.eq.0) derxBy => derx_00 + if (nclxBy1.eq.1.and.nclxByn.eq.1) derxBy => derx_11 + if (nclxBy1.eq.1.and.nclxByn.eq.2) derxBy => derx_12 + if (nclxBy1.eq.2.and.nclxByn.eq.1) derxBy => derx_21 + if (nclxBy1.eq.2.and.nclxByn.eq.2) derxBy => derx_22 + ! + if (nclyBy1.eq.0.and.nclyByn.eq.0) deryBy => dery_00 + if (nclyBy1.eq.1.and.nclyByn.eq.1) deryBy => dery_11 + if (nclyBy1.eq.1.and.nclyByn.eq.2) deryBy => dery_12 + if (nclyBy1.eq.2.and.nclyByn.eq.1) deryBy => dery_21 + if (nclyBy1.eq.2.and.nclyByn.eq.2) deryBy => dery_22 + ! + if (nclzBy1.eq.0.and.nclzByn.eq.0) derzBy => derz_00 + if (nclzBy1.eq.1.and.nclzByn.eq.1) derzBy => derz_11 + if (nclzBy1.eq.1.and.nclzByn.eq.2) derzBy => derz_12 + if (nclzBy1.eq.2.and.nclzByn.eq.1) derzBy => derz_21 + if (nclzBy1.eq.2.and.nclzByn.eq.2) derzBy => derz_22 + ! SecondByderivative + if (nclxBy1.eq.0.and.nclxByn.eq.0) derxxBy => derxx_00 + if (nclxBy1.eq.1.and.nclxByn.eq.1) derxxBy => derxx_11 + if (nclxBy1.eq.1.and.nclxByn.eq.2) derxxBy => derxx_12 + if (nclxBy1.eq.2.and.nclxByn.eq.1) derxxBy => derxx_21 + if (nclxBy1.eq.2.and.nclxByn.eq.2) derxxBy => derxx_22 + !y + if (nclyBy1.eq.0.and.nclyByn.eq.0) deryyBy => deryy_00 + if (nclyBy1.eq.1.and.nclyByn.eq.1) deryyBy => deryy_11 + if (nclyBy1.eq.1.and.nclyByn.eq.2) deryyBy => deryy_12 + if (nclyBy1.eq.2.and.nclyByn.eq.1) deryyBy => deryy_21 + if (nclyBy1.eq.2.and.nclyByn.eq.2) deryyBy => deryy_22 + !z + if (nclzBy1.eq.0.and.nclzByn.eq.0) derzzBy => derzz_00 + if (nclzBy1.eq.1.and.nclzByn.eq.1) derzzBy => derzz_11 + if (nclzBy1.eq.1.and.nclzByn.eq.2) derzzBy => derzz_12 + if (nclzBy1.eq.2.and.nclzByn.eq.1) derzzBy => derzz_21 + if (nclzBy1.eq.2.and.nclzByn.eq.2) derzzBy => derzz_22 + call first_derivative(alfa1x,af1x,bf1x,cf1x,df1x,alfa2x,af2x,alfanx,afnx,bfnx,& + cfnx,dfnx,alfamx,afmx,alfaix,afix,bfix,& + ffxB(:,2),fsxB(:,2),fwxB(:,2),ffxpB(:,2),fsxpB(:,2),fwxpB(:,2),dx,nx,nclxBy1,nclxByn) + call first_derivative(alfa1y,af1y,bf1y,cf1y,df1y,alfa2y,af2y,alfany,afny,bfny,& + cfny,dfny,alfamy,afmy,alfajy,afjy,bfjy,& + ffyB(:,2),fsyB(:,2),fwyB(:,2),ffypB(:,2),fsypB(:,2),fwypB(:,2),dy,ny,nclyBy1,nclyByn) + call first_derivative(alfa1z,af1z,bf1z,cf1z,df1z,alfa2z,af2z,alfanz,afnz,bfnz,& + cfnz,dfnz,alfamz,afmz,alfakz,afkz,bfkz,& + ffzB(:,2),fszB(:,2),fwzB(:,2),ffzpB(:,2),fszpB(:,2),fwzpB(:,2),dz,nz,nclzBy1,nclzByn) + call second_derivative(alsa1x,as1x,bs1x,& + cs1x,ds1x,alsa2x,as2x,alsanx,asnx,bsnx,csnx,dsnx,alsamx,& + asmx,alsa3x,as3x,bs3x,alsatx,astx,bstx,& + alsa4x,as4x,bs4x,cs4x,& + alsattx,asttx,bsttx,csttx,& + alsaix,asix,bsix,csix,dsix,& + sfxB(:,2),ssxB(:,2),swxB(:,2),sfxpB(:,2),ssxpB(:,2),swxpB(:,2),dx2,nx,nclxBy1,nclxByn) + call second_derivative(alsa1y,as1y,bs1y,& + cs1y,ds1y,alsa2y,as2y,alsany,asny,bsny,csny,dsny,alsamy,& + asmy,alsa3y,as3y,bs3y,alsaty,asty,bsty,& + alsa4y,as4y,bs4y,cs4y,& + alsatty,astty,bstty,cstty,& + alsajy,asjy,bsjy,csjy,dsjy,& + sfyB(:,2),ssyB(:,2),swyB(:,2),sfypB(:,2),ssypB(:,2),swypB(:,2),dy2,ny,nclyBy1,nclyByn) + call second_derivative(alsa1z,as1z,bs1z,& + cs1z,ds1z,alsa2z,as2z,alsanz,asnz,bsnz,csnz,dsnz,alsamz,& + asmz,alsa3z,as3z,bs3z,alsatz,astz,bstz,& + alsa4z,as4z,bs4z,cs4z,& + alsattz,asttz,bsttz,csttz,& + alsakz,askz,bskz,cskz,dskz,& + sfzB(:,2),sszB(:,2),swzB(:,2),sfzpB(:,2),sszpB(:,2),swzpB(:,2),dz2,nz,nclzBy1,nclzByn) + ! First derivative + if (nclxBz1.eq.0.and.nclxBzn.eq.0) derxBz => derx_00 + if (nclxBz1.eq.1.and.nclxBzn.eq.1) derxBz => derx_11 + if (nclxBz1.eq.1.and.nclxBzn.eq.2) derxBz => derx_12 + if (nclxBz1.eq.2.and.nclxBzn.eq.1) derxBz => derx_21 + if (nclxBz1.eq.2.and.nclxBzn.eq.2) derxBz => derx_22 + ! + if (nclyBz1.eq.0.and.nclyBzn.eq.0) deryBz => dery_00 + if (nclyBz1.eq.1.and.nclyBzn.eq.1) deryBz => dery_11 + if (nclyBz1.eq.1.and.nclyBzn.eq.2) deryBz => dery_12 + if (nclyBz1.eq.2.and.nclyBzn.eq.1) deryBz => dery_21 + if (nclyBz1.eq.2.and.nclyBzn.eq.2) deryBz => dery_22 + ! + if (nclzBz1.eq.0.and.nclzBzn.eq.0) derzBz => derz_00 + if (nclzBz1.eq.1.and.nclzBzn.eq.1) derzBz => derz_11 + if (nclzBz1.eq.1.and.nclzBzn.eq.2) derzBz => derz_12 + if (nclzBz1.eq.2.and.nclzBzn.eq.1) derzBz => derz_21 + if (nclzBz1.eq.2.and.nclzBzn.eq.2) derzBz => derz_22 + ! SecondBzderivative + if (nclxBz1.eq.0.and.nclxBzn.eq.0) derxxBz => derxx_00 + if (nclxBz1.eq.1.and.nclxBzn.eq.1) derxxBz => derxx_11 + if (nclxBz1.eq.1.and.nclxBzn.eq.2) derxxBz => derxx_12 + if (nclxBz1.eq.2.and.nclxBzn.eq.1) derxxBz => derxx_21 + if (nclxBz1.eq.2.and.nclxBzn.eq.2) derxxBz => derxx_22 + !y + if (nclyBz1.eq.0.and.nclyBzn.eq.0) deryyBz => deryy_00 + if (nclyBz1.eq.1.and.nclyBzn.eq.1) deryyBz => deryy_11 + if (nclyBz1.eq.1.and.nclyBzn.eq.2) deryyBz => deryy_12 + if (nclyBz1.eq.2.and.nclyBzn.eq.1) deryyBz => deryy_21 + if (nclyBz1.eq.2.and.nclyBzn.eq.2) deryyBz => deryy_22 + !z + if (nclzBz1.eq.0.and.nclzBzn.eq.0) derzzBz => derzz_00 + if (nclzBz1.eq.1.and.nclzBzn.eq.1) derzzBz => derzz_11 + if (nclzBz1.eq.1.and.nclzBzn.eq.2) derzzBz => derzz_12 + if (nclzBz1.eq.2.and.nclzBzn.eq.1) derzzBz => derzz_21 + if (nclzBz1.eq.2.and.nclzBzn.eq.2) derzzBz => derzz_22 + call first_derivative(alfa1x,af1x,bf1x,cf1x,df1x,alfa2x,af2x,alfanx,afnx,bfnx,& + cfnx,dfnx,alfamx,afmx,alfaix,afix,bfix,& + ffxB(:,3),fsxB(:,3),fwxB(:,3),ffxpB(:,3),fsxpB(:,3),fwxpB(:,3),dx,nx,nclxBz1,nclxBzn) + call first_derivative(alfa1y,af1y,bf1y,cf1y,df1y,alfa2y,af2y,alfany,afny,bfny,& + cfny,dfny,alfamy,afmy,alfajy,afjy,bfjy,& + ffyB(:,3),fsyB(:,3),fwyB(:,3),ffypB(:,3),fsypB(:,3),fwypB(:,3),dy,ny,nclyBz1,nclyBzn) + call first_derivative(alfa1z,af1z,bf1z,cf1z,df1z,alfa2z,af2z,alfanz,afnz,bfnz,& + cfnz,dfnz,alfamz,afmz,alfakz,afkz,bfkz,& + ffzB(:,3),fszB(:,3),fwzB(:,3),ffzpB(:,3),fszpB(:,3),fwzpB(:,3),dz,nz,nclzBz1,nclzBzn) + call second_derivative(alsa1x,as1x,bs1x,& + cs1x,ds1x,alsa2x,as2x,alsanx,asnx,bsnx,csnx,dsnx,alsamx,& + asmx,alsa3x,as3x,bs3x,alsatx,astx,bstx,& + alsa4x,as4x,bs4x,cs4x,& + alsattx,asttx,bsttx,csttx,& + alsaix,asix,bsix,csix,dsix,& + sfxB(:,3),ssxB(:,3),swxB(:,3),sfxpB(:,3),ssxpB(:,3),swxpB(:,3),dx2,nx,nclxBz1,nclxBzn) + call second_derivative(alsa1y,as1y,bs1y,& + cs1y,ds1y,alsa2y,as2y,alsany,asny,bsny,csny,dsny,alsamy,& + asmy,alsa3y,as3y,bs3y,alsaty,asty,bsty,& + alsa4y,as4y,bs4y,cs4y,& + alsatty,astty,bstty,cstty,& + alsajy,asjy,bsjy,csjy,dsjy,& + sfyB(:,3),ssyB(:,3),swyB(:,3),sfypB(:,3),ssypB(:,3),swypB(:,3),dy2,ny,nclyBz1,nclyBzn) + call second_derivative(alsa1z,as1z,bs1z,& + cs1z,ds1z,alsa2z,as2z,alsanz,asnz,bsnz,csnz,dsnz,alsamz,& + asmz,alsa3z,as3z,bs3z,alsatz,astz,bstz,& + alsa4z,as4z,bs4z,cs4z,& + alsattz,asttz,bsttz,csttz,& + alsakz,askz,bskz,cskz,dskz,& + sfzB(:,3),sszB(:,3),swzB(:,3),sfzpB(:,3),sszpB(:,3),swzpB(:,3),dz2,nz,nclzBz1,nclzBzn) + endif + call interpolation(dx,nxm,nx,nclx1,nclxn,& alcaix6,acix6,bcix6,& ailcaix6,aicix6,bicix6,cicix6,dicix6,& diff --git a/src/tools.f90 b/src/tools.f90 index 782c341e1..c755bf63f 100644 --- a/src/tools.f90 +++ b/src/tools.f90 @@ -780,6 +780,7 @@ subroutine compute_cfldiff() use param, only : cfl_diff_sum, cfl_diff_x, cfl_diff_y, cfl_diff_z use variables, only : dyp use decomp_2d, only : nrank + use mhd, only: mhd_active, mhd_equation,rem implicit none @@ -804,6 +805,29 @@ subroutine compute_cfldiff() write(*,*) '===========================================================' endif + if( mhd_active.and.mhd_equation) then + + cfl_diff_x = dt/ (dx**2) / rem + cfl_diff_z = dt/ (dz**2) / rem + + if (istret == 0) then + cfl_diff_y = dt / (dy**2) / rem + else + cfl_diff_y = dt / (minval(dyp)**2) / rem + end if + + cfl_diff_sum = cfl_diff_x + cfl_diff_y + cfl_diff_z + + if (nrank==0) then + write(*,*) '===========================================================' + write(*,*) 'Magnetic Diffusion number' + write(*,"(' B cfl_diff_x : ',F13.8)") cfl_diff_x + write(*,"(' B cfl_diff_y : ',F13.8)") cfl_diff_y + write(*,"(' B cfl_diff_z : ',F13.8)") cfl_diff_z + write(*,"(' B cfl_diff_sum : ',F13.8)") cfl_diff_sum + write(*,*) '===========================================================' + endif + endif return end subroutine compute_cfldiff !################################################################## diff --git a/src/variables.f90 b/src/variables.f90 index bb3917716..d5188b5c2 100644 --- a/src/variables.f90 +++ b/src/variables.f90 @@ -878,6 +878,86 @@ subroutine init_variables allocate(swzpS(nz)) swzpS=zero + !MHD -- todo if active + allocate(ffxB(nx,3)) + ffxB=zero + allocate(sfxB(nx,3)) + sfxB=zero + allocate(fsxB(nx,3)) + fsxB=zero + allocate(fwxB(nx,3)) + fwxB=zero + allocate(ssxB(nx,3)) + ssxB=zero + allocate(swxB(nx,3)) + swxB=zero + + allocate(ffxpB(nx,3)) + ffxpB=zero + allocate(sfxpB(nx,3)) + sfxpB=zero + allocate(fsxpB(nx,3)) + fsxpB=zero + allocate(fwxpB(nx,3)) + fwxpB=zero + allocate(ssxpB(nx,3)) + ssxpB=zero + allocate(swxpB(nx,3)) + swxpB=zero + + allocate(ffyB(ny,3)) + ffyB=zero + allocate(sfyB(ny,3)) + sfyB=zero + allocate(fsyB(ny,3)) + fsyB=zero + allocate(fwyB(ny,3)) + fwyB=zero + allocate(ssyB(ny,3)) + ssyB=zero + allocate(swyB(ny,3)) + swyB=zero + + allocate(ffypB(ny,3)) + ffypB=zero + allocate(sfypB(ny,3)) + sfypB=zero + allocate(fsypB(ny,3)) + fsypB=zero + allocate(fwypB(ny,3)) + fwypB=zero + allocate(ssypB(ny,3)) + ssypB=zero + allocate(swypB(ny,3)) + swypB=zero + + allocate(ffzB(nz,3)) + ffzB=zero + allocate(sfzB(nz,3)) + sfzB=zero + allocate(fszB(nz,3)) + fszB=zero + allocate(fwzB(nz,3)) + fwzB=zero + allocate(sszB(nz,3)) + sszB=zero + allocate(swzB(nz,3)) + swzB=zero + + allocate(ffzpB(nz,3)) + ffzpB=zero + allocate(sfzpB(nz,3)) + sfzpB=zero + allocate(fszpB(nz,3)) + fszpB=zero + allocate(fwzpB(nz,3)) + fwzpB=zero + allocate(sszpB(nz,3)) + sszpB=zero + allocate(swzpB(nz,3)) + swzpB=zero + + allocate(sx(xsize(2),xsize(3))) sx=zero allocate(vx(xsize(2),xsize(3)))