Commit ba5f5ea6 authored by Alexander Neukirchen's avatar Alexander Neukirchen

Experimental changes considering source free B-fields.

parent abce4e21
......@@ -106,7 +106,7 @@ CONTAINS
obsolete,cell,oneD,sliceplot,mpi,noco,den,denRot,EnergyDen,vTot,vx,results)
!ToDo, check if this is needed for more potentials as well...
CALL vgen_finalize(atoms,stars,vacuum,sym,noco,input,sphhar,vTot,vCoul,denRot)
CALL vgen_finalize(mpi,dimension,oneD,field,cell,atoms,stars,vacuum,sym,noco,input,sphhar,vTot,vCoul,denRot)
!DEALLOCATE(vcoul%pw_w)
CALL bfield(input,noco,atoms,field,vTot)
......
......@@ -510,9 +510,16 @@ MODULE m_intgr
REAL :: dr, alpha
logr=LOG(x)
dr=x(2)-x(1)
dr=logr(2)-logr(1)
rf=x*y
IF (y(1)*y(2).GT.zero) THEN
alpha = 1.0 + log(y(2)/y(1))/dr
IF (alpha.GT.zero) z(1) = x(1)*y(1)/alpha
ENDIF
z=z(1)
DO i=2, jri
z(i:)=z(i:)+(rf(i-1)+rf(i))*dr/2.0
END DO
......@@ -585,4 +592,57 @@ MODULE m_intgr
RETURN
END SUBROUTINE intgr5
SUBROUTINE intgr6(y,x,h,jri,z)
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: h
REAL, INTENT (IN) :: x(jri), y(jri)
REAL, INTENT (OUT):: z(jri)
REAL :: dr, r(7), alpha
INTEGER :: i, j
REAL :: yr(7)
z(1) = zero
IF (y(1)*y(2).GT.zero) THEN
alpha = 1.0 + log(y(2)/y(1))/h
IF (alpha.GT.zero) z(1) = x(1)*y(1)/alpha
ENDIF
z(2)=z(1)+h*(x(1)*y(1)+x(2)*y(2))/2
z(3)=z(1)+h*(x(1)*y(1)+4*x(2)*y(2)+x(3)*y(3))/3.
DO j = 4,jri
z(j) = z(j-3) + 3*h*(x(j-3)*y(j-3)+3*x(j-2)*y(j-2)+3*x(j-1)*y(j-1)+x(j)*y(j))/8.
ENDDO
RETURN
END SUBROUTINE intgr6
SUBROUTINE intgr7(y,x,h,jri,z)
INTEGER, INTENT (IN) :: jri
REAL, INTENT (IN) :: h
REAL, INTENT (IN) :: x(jri), y(jri)
REAL, INTENT (OUT):: z(jri)
REAL :: dr, r(7), alpha
INTEGER :: i, j
REAL :: yr(7)
z(1) = zero
IF (y(1)*y(2).GT.zero) THEN
alpha = 1.0 + log(y(2)/y(1))/h
IF (alpha.GT.zero) z(1) = x(1)*y(1)/alpha
ENDIF
z(2)=z(1)+h*(x(1)*y(1)+x(2)*y(2))/2
DO j = 3,jri
z(j) = z(j-2) + h*(x(j-2)*y(j-2)+4*x(j-1)*y(j-1)+x(j)*y(j))/3.
ENDDO
RETURN
END SUBROUTINE intgr7
END MODULE m_intgr
MODULE m_sfTests
IMPLICIT NONE
CONTAINS
SUBROUTINE plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, &
noco, xcB, div, phi, cvec, corrB, div2)
USE m_plot
IMPLICIT NONE
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
......@@ -53,13 +52,11 @@ CONTAINS
END SUBROUTINE plotBtest
SUBROUTINE buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,itest,Avec,icut,denMat,factor)
SUBROUTINE buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,itest,Avec,denMat,factor)
USE m_mt_tofrom_grid
USE m_pw_tofrom_grid
USE m_xcBfield
IMPLICIT NONE
TYPE(t_stars), INTENT(IN) :: stars
TYPE(t_atoms), INTENT(IN) :: atoms
TYPE(t_sphhar), INTENT(IN) :: sphhar
......@@ -70,7 +67,6 @@ CONTAINS
TYPE(t_cell), INTENT(IN) :: cell
INTEGER, INTENT(IN) :: itest
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: Avec
INTEGER, INTENT(OUT) :: icut(atoms%ntype)
TYPE(t_potden), OPTIONAL, INTENT(IN) :: denMat
REAL, OPTIONAL, INTENT(IN) :: factor
......@@ -82,14 +78,13 @@ CONTAINS
REAL :: vec1(3), vec2(3), vec3(3), zero(3), point(3)
INTEGER :: grid(3)
icut=1
IF (itest.EQ.0) THEN
RETURN
END IF
IF (PRESENT(denMat)) THEN
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denMat,factor,Avec,icut)
CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denMat,factor,Avec)
RETURN
END IF
......@@ -229,7 +224,7 @@ CONTAINS
TYPE(t_potden), DIMENSION(3) :: aVec, cvec, corrB
TYPE(t_potden) :: div, phi, checkdiv
INTEGER :: i, n, lh, l, icut(3,0:sphhar%nlhd,atoms%ntype)
INTEGER :: i, n, lh, l
REAL :: g(atoms%jmtd)
REAL :: radii(atoms%jmtd,atoms%ntype), funcsr(atoms%jmtd,atoms%ntype,3), trueder(atoms%jmtd,atoms%ntype,3), trueint(atoms%jmtd,atoms%ntype,3), &
......@@ -243,11 +238,11 @@ CONTAINS
! whether the sourcefree routine made it sourcefree.
IF (PRESENT(denMat)) THEN
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,1,aVec,icut,denMat,factor)
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,icut,div,phi,cvec,corrB,checkdiv)
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,1,aVec,denMat,factor)
CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,div,phi,cvec,corrB,checkdiv)
CALL plotBtest(stars, atoms, sphhar, vacuum, input, oneD, sym, cell, noco, aVec, div, phi, cvec, corrB, checkdiv)
ELSE
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,0,aVec,icut)
CALL buildAtest(stars,atoms,sphhar,vacuum,input,noco,sym,cell,0,aVec)
RETURN
END IF
......@@ -265,7 +260,7 @@ CONTAINS
SUBROUTINE difftester(atoms,n,l,f,g)
USE m_types
USE m_grdchlh
IMPLICIT NONE
TYPE(t_atoms), INTENT(IN) :: atoms
INTEGER, INTENT(IN) :: n,l
REAL, INTENT(IN) :: f(atoms%jri(n))
......
......@@ -5,8 +5,9 @@
!--------------------------------------------------------------------------------
MODULE m_vgen_finalize
USE m_juDFT
USE m_xcBfield
CONTAINS
SUBROUTINE vgen_finalize(atoms,stars,vacuum,sym,noco,input,sphhar,vTot,vCoul,denRot)
SUBROUTINE vgen_finalize(mpi,dimension,oneD,field,cell,atoms,stars,vacuum,sym,noco,input,sphhar,vTot,vCoul,denRot)
! ***********************************************************
! FLAPW potential generator *
! ***********************************************************
......@@ -20,6 +21,11 @@ CONTAINS
USE m_rotate_mt_den_tofrom_local
USE m_sfTests
IMPLICIT NONE
TYPE(t_mpi), INTENT(IN) :: mpi
TYPE(t_dimension), INTENT(IN) :: dimension
TYPE(t_oneD), INTENT(IN) :: oneD
TYPE(t_field), INTENT(INOUT) :: field
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_vacuum),INTENT(IN) :: vacuum
TYPE(t_noco),INTENT(IN) :: noco
TYPE(t_sym),INTENT(IN) :: sym
......@@ -32,6 +38,9 @@ CONTAINS
! .. Local Scalars ..
INTEGER i,js,n
TYPE(t_potden) :: div, phi, checkdiv
TYPE(t_potden), DIMENSION(3) :: cvec, corrB, bxc
! Rescale vTot%pw_w with number of stars
......@@ -48,6 +57,17 @@ CONTAINS
END IF
ENDIF
! Source-free testwise
!CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,inDen,1.0)
!CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,vTot,2.0)
! Once it is tested:
!IF (noco%l_mtnocoPot.AND.(.FALSE.)) THEN ! l_sf will go here
!CALL makeVectorField(stars,atoms,sphhar,vacuum,input,noco,vTot,2.0,bxc)
!CALL sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,bxc,div,phi,cvec,corrB,checkdiv)
!CALL correctPot(vTot,cvec)
!END IF
! ---> store v(l=0) component as r*v(l=0)/sqrt(4pi)
DO js = 1,input%jspins !Used input%jspins instead of SIZE(vtot%mt,4) since the off diag, elements of VTot%mt need no rescaling.
......@@ -66,25 +86,12 @@ CONTAINS
!Copy first vacuum into second vacuum if this was not calculated before
IF (vacuum%nvac==1) THEN
vTot%vacz(:,2,:) = vTot%vacz(:,1,:)
! DO i=1,3
! xcB(i)%vacz(:,2,:) = xcB(i)%vacz(:,1,:)
! ENDDO
IF (sym%invs) THEN
vTot%vacxy(:,:,2,:) = CMPLX(vTot%vacxy(:,:,1,:))
! DO i=1,3
! xcB(i)%vacxy(:,:,2,:) = CMPLX(xcB(i)%vacxy(:,:,1,:))
! ENDDO
ELSE
vTot%vacxy(:,:,2,:) = vTot%vacxy(:,:,1,:)
! DO i=1,3
! xcB(i)%vacxy(:,:,2,:) = xcB(i)%vacxy(:,:,1,:)
! ENDDO
ENDIF
ENDIF
! Source-free testwise
!CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,inDen,1.0)
!CALL sftest(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,1,vTot,2.0)
END SUBROUTINE vgen_finalize
END MODULE m_vgen_finalize
......@@ -61,7 +61,7 @@ contains
real :: green_factor, termsR
real :: green_1 (1:atoms%jmtd), green_2 (1:atoms%jmtd)
real :: integrand_1(1:atoms%jmtd), integrand_2(1:atoms%jmtd)
real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)
real :: integral_1 (1:atoms%jmtd), integral_2 (1:atoms%jmtd)!, integral_3 (1:atoms%jmtd)
real :: sbf(0:atoms%lmaxd)
real, allocatable, dimension(:,:) :: il, kl
......@@ -168,24 +168,26 @@ contains
end if
integrand_1(1:imax) = green_1(1:imax) * rho(1:imax,lh,n)
integrand_2(1:imax) = green_2(1:imax) * rho(1:imax,lh,n)
call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
if (.not.dosf) THEN
call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
! Source-free testwise
!if (dosf) then
else
!if (l==5) THEN
! integrand_2(1:300)=integrand_2(1:300)*(atoms%rmsh(1:300,n)/atoms%rmsh(300,n))**l
!end if
!call intgrt(integrand_1(1:imax),atoms%rmsh(:,n),imax,integral_1(1:imax))
!call intgrt(integrand_2(1:imax),atoms%rmsh(:,n),imax,integral_2(1:imax))
!vtl(lh,n)=(0.0,0.0)
! integrand_2(1:300)=integrand_2(1:300)*(atoms%rmsh(1:300,n)/atoms%rmsh(300,n))**(4*l)
!end if
call intgrtlog(integrand_1(1:imax),atoms%rmsh(:,n),imax,integral_1(1:imax))
call intgrtlog(integrand_2(1:imax),atoms%rmsh(:,n),imax,integral_2(1:imax))
!call intgr4(integrand_1(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_1(1:imax))
!call intgr4(integrand_2(1:imax),atoms%rmsh(:,n),atoms%dx(n),imax,integral_2(1:imax))
!call intgr2( integrand_1(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_1(1:imax) )
!call intgr2( integrand_2(1:imax), atoms%rmsh(1,n), atoms%dx(n), imax, integral_2(1:imax) )
!end if
!integral_3(1:imax)=(integral_2(imax)-integral_2(1:imax))*green_1(imax)
end if
termsR = integral_2(imax) + ( vtl(lh,n) / green_factor - integral_1(imax) * green_2(imax) ) / green_1(imax)
vr(1:imax,lh,n) = green_factor * ( green_1(1:imax) * ( termsR - integral_2(1:imax) ) &
+ green_2(1:imax) * integral_1(1:imax) )
end do
nat = nat + atoms%neq(n)
end do
......
......@@ -8,6 +8,8 @@ MODULE m_xcBfield
USE m_constants
USE m_plot
USE m_divergence
IMPLICIT NONE
!-----------------------------------------------------------------------------
! This module contains all the operations on exchange-correlation B-fields
......@@ -19,7 +21,7 @@ MODULE m_xcBfield
PUBLIC :: makeBxc, sourcefree
CONTAINS
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec,icut)
SUBROUTINE makeVectorField(stars,atoms,sphhar,vacuum,input,noco,denmat,factor,aVec)
! Contructs the exchange-correlation magnetic field from the total poten-
! tial matrix or the magnetic density for the density matrix. Only used for
......@@ -44,14 +46,12 @@ CONTAINS
TYPE(t_potden), INTENT(IN) :: denmat
REAL, INTENT(IN) :: factor
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: aVec
INTEGER, INTENT(OUT) :: icut(3,0:sphhar%nlhd,atoms%ntype)
TYPE(t_potden), DIMENSION(4) :: dummyDen
INTEGER :: i, itype, ir, lh
REAL :: r2(atoms%jmtd), fcut
fcut=1.e-12
icut=1
! Initialize and fill a dummy density array, that takes the initial result
! of matrixsplit.
......@@ -80,7 +80,6 @@ CONTAINS
r2=atoms%rmsh(:,itype)**2
END IF
aVec(i)%mt(:,lh,itype,1) = aVec(i)%mt(:,lh,itype,1)*r2
!WHERE (ABS(aVec(i)%mt(ir,0:,itype,:)/r2) < fcut) aVec(i)%mt(ir,0:,itype,:) = 0.0
END DO !lh
END DO !itype
aVec(i)%pw(1:,:) = dummyDen(i+1)%pw(1:,:)
......@@ -90,7 +89,7 @@ CONTAINS
END SUBROUTINE makeVectorField
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,icut,div,phi,cvec,corrB,checkdiv)
SUBROUTINE sourcefree(mpi,dimension,field,stars,atoms,sphhar,vacuum,input,oneD,sym,cell,noco,aVec,div,phi,cvec,corrB,checkdiv)
USE m_vgen_coulomb
USE m_gradYlm
USE m_grdchlh
......@@ -121,7 +120,6 @@ CONTAINS
TYPE(t_cell), INTENT(IN) :: cell
TYPE(t_noco), INTENT(IN) :: noco
TYPE(t_potden), DIMENSION(3), INTENT(INOUT) :: aVec
INTEGER, INTENT(IN) :: icut(3,0:sphhar%nlhd,atoms%ntype)
TYPE(t_potden), INTENT(OUT) :: div, phi, checkdiv
TYPE(t_potden), DIMENSION(3), INTENT(OUT) :: cvec, corrB
......@@ -130,34 +128,17 @@ CONTAINS
INTEGER :: n, jr, lh, lhmax, jcut, nat
REAL :: xp(3,dimension%nspd)
!DO i=1,3
! aVec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
!END DO
CALL div%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype, &
atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN, &
vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(div%pw_w,mold=div%pw)
CALL divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,aVec,div)
atloc=atoms
atloc%zatom=0.0 !Local atoms variable with no charges; needed for the potential generation.
fcut=1.e-6
!div%mt(:300,25,:,1)=0.0
!DO n=1,atoms%ntype
!lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
!DO lh=0, lhmax
!div%mt(:,lh,n,1)=div%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
!WHERE (ABS(div%mt(:,lh,n,1))<MAXVAL(ABS(div%mt(:,lh,n,1)))*fcut) div%mt(:,lh,n,1)=0.0
!END DO
!DO lh=0, lhmax
!div%mt(:,lh,n,1)=div%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
!END DO
!END DO
! Local atoms variable with no charges;
! needed for the potential generation from the divergence.
atloc=atoms
atloc%zatom=0.0
CALL phi%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_POTCOUL,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(phi%pw_w(SIZE(phi%pw,1),size(phi%pw,2)))
phi%pw_w = CMPLX(0.0,0.0)
......@@ -171,10 +152,6 @@ CONTAINS
CALL divpotgrad2(stars,atoms,sphhar,vacuum,sym,cell,noco,phi,cvec)
!DO i=1,3
! cvec(i)%mt(:,atoms%lmaxd**2:,:,:)=0.0
!END DO
DO i=1,3
CALL corrB(i)%init_potden_simple(stars%ng3,atoms%jmtd,sphhar%nlhd,atoms%ntype,atoms%n_u,1,.FALSE.,.FALSE.,POTDEN_TYPE_DEN,vacuum%nmzd,vacuum%nmzxyd,stars%ng2)
ALLOCATE(corrB(i)%pw_w,mold=corrB(i)%pw)
......@@ -188,32 +165,13 @@ CONTAINS
CALL divergence2(stars,atoms,sphhar,vacuum,sym,cell,noco,corrB,checkdiv)
nat=1
DO n=1,atoms%ntype
lhmax=sphhar%nlh(atoms%ntypsy(SUM(atoms%neq(:n - 1)) + 1))
DO lh=0, lhmax
aVec(3)%mt(:,lh,n,1)=aVec(3)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
div%mt(:,lh,n,1)=div%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(1)%mt(:,lh,n,1)=cvec(1)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(2)%mt(:,lh,n,1)=cvec(2)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
cvec(3)%mt(:,lh,n,1)=cvec(3)%mt(:,lh,n,1)/(atoms%rmsh(:, n)**2)
END DO
DO lh=0, lhmax
aVec(3)%mt(:,lh,n,1)=aVec(3)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
div%mt(:,lh,n,1)=div%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
checkdiv%mt(:,lh,n,1)=checkdiv%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(1)%mt(:,lh,n,1)=cvec(1)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(2)%mt(:,lh,n,1)=cvec(2)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
cvec(3)%mt(:,lh,n,1)=cvec(3)%mt(:,lh,n,1)*(atoms%rmsh(:, n)**2)
END DO
CALL sphpts(xp,dimension%nspd,atoms%rmt(n),atoms%pos(1,nat))
CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,dimension,atoms,sphhar,stars,sym,vacuum,cell,oneD,div)
CALL checkdop(xp,dimension%nspd,n,nat,0,-1,1,dimension,atoms,sphhar,stars,sym,vacuum,cell,oneD,phi)
nat = nat + atoms%neq(n)
END DO
END SUBROUTINE sourcefree
......@@ -233,6 +191,8 @@ CONTAINS
TYPE(t_potden), INTENT(INOUT) :: vTot
TYPE(t_potden), DIMENSION(3), INTENT(IN) :: c
REAL :: pwr(SIZE(vTot%pw(1:,3))), pwi(SIZE(vTot%pw(1:,3)))
vTot%mt(:,0:,:,1)=vTot%mt(:,0:,:,1)+c(3)%mt(:,0:,:,1)
vTot%mt(:,0:,:,2)=vTot%mt(:,0:,:,2)-c(3)%mt(:,0:,:,1)
vTot%mt(:,0:,:,3)=vTot%mt(:,0:,:,3)+c(1)%mt(:,0:,:,1)
......@@ -240,18 +200,25 @@ CONTAINS
vTot%pw(1:,1)=vTot%pw(1:,1)+c(3)%pw(1:,1)
vTot%pw(1:,2)=vTot%pw(1:,2)-c(3)%pw(1:,1)
vTot%pw(1:,3)=vTot%pw(1:,3)+c(1)%pw(1:,1)
vTot%pw(1:,4)=vTot%pw(1:,4)+c(2)%pw(1:,1)
vTot%vacz(1:,1:,1)=vTot%vacz(1:,1:,1)+c(3)%vacz(1:,1:,1)
vTot%vacz(1:,1:,2)=vTot%vacz(1:,1:,2)-c(3)%vacz(1:,1:,1)
vTot%vacz(1:,1:,3)=vTot%vacz(1:,1:,3)+c(1)%vacz(1:,1:,1)
vTot%vacz(1:,1:,4)=vTot%vacz(1:,1:,4)+c(2)%vacz(1:,1:,1)
vTot%vacxy(1:,1:,1:,1)=vTot%vacxy(1:,1:,1:,1)+c(3)%vacxy(1:,1:,1:,1)
vTot%vacxy(1:,1:,1:,2)=vTot%vacxy(1:,1:,1:,2)-c(3)%vacxy(1:,1:,1:,1)
vTot%vacxy(1:,1:,1:,3)=vTot%vacxy(1:,1:,1:,3)+c(1)%vacxy(1:,1:,1:,1)
vTot%vacxy(1:,1:,1:,4)=vTot%vacxy(1:,1:,1:,4)+c(2)%vacxy(1:,1:,1:,1)
pwr=REAL(vTot%pw(1:,3)) + REAL(c(1)%pw(1:,1)) - AIMAG(c(2)%pw(1:,1))
pwi=AIMAG(vTot%pw(1:,3)) + AIMAG(c(1)%pw(1:,1)) + REAL(c(2)%pw(1:,1))
vTot%pw(1:,3)=CMPLX(pwr,pwi)
vTot%pw_w(1:,1)=vTot%pw_w(1:,1)+c(3)%pw_w(1:,1)
vTot%pw_w(1:,2)=vTot%pw_w(1:,2)-c(3)%pw_w(1:,1)
pwr=REAL(vTot%pw_w(1:,3)) + REAL(c(1)%pw_w(1:,1)) - AIMAG(c(2)%pw_w(1:,1))
pwi=AIMAG(vTot%pw_w(1:,3)) + AIMAG(c(1)%pw_w(1:,1)) + REAL(c(2)%pw_w(1:,1))
vTot%pw_w(1:,3)=CMPLX(pwr,pwi)
!vTot%vacz(1:,1:,1)=vTot%vacz(1:,1:,1)+c(3)%vacz(1:,1:,1)
!vTot%vacz(1:,1:,2)=vTot%vacz(1:,1:,2)-c(3)%vacz(1:,1:,1)
!vTot%vacz(1:,1:,3)=vTot%vacz(1:,1:,3)+c(1)%vacz(1:,1:,1)
!vTot%vacz(1:,1:,4)=vTot%vacz(1:,1:,4)+c(2)%vacz(1:,1:,1)
!vTot%vacxy(1:,1:,1:,1)=vTot%vacxy(1:,1:,1:,1)+c(3)%vacxy(1:,1:,1:,1)
!vTot%vacxy(1:,1:,1:,2)=vTot%vacxy(1:,1:,1:,2)-c(3)%vacxy(1:,1:,1:,1)
!vTot%vacxy(1:,1:,1:,3)=vTot%vacxy(1:,1:,1:,3)+c(1)%vacxy(1:,1:,1:,1)
!vTot%vacxy(1:,1:,1:,4)=vTot%vacxy(1:,1:,1:,4)+c(2)%vacxy(1:,1:,1:,1)
END SUBROUTINE
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment