#include "macrodef.fh" !#define OLD ! 0: rk4, 1: Adams-Moulton #define RKorAM 0 function beta_rhs(xx,CJx,Kx) result(gont) implicit none double complex,intent(in) :: CJx real*8,intent(in) :: xx,Kx real*8 :: gont gont = xx*(1.d0-xx)/8.d0*(dreal(CJx*dconjg(CJx))-Kx*Kx) return end function beta_rhs function Q_rhs(xx,CJ,CJx,DCJx,KK,Ck,Ckx,Cnux,KKx,CBx,Cnu,DCJ,CB,CQ) result(gont) implicit none double complex,intent(in) :: CJ,CJx,DCJx,Ck,Ckx,Cnux,CBx,Cnu,DCJ,CB,CQ real*8,intent(in) :: xx,KK,KKx double complex :: gont gont = -KK*(Ckx+Cnux)+Cnu*KKx+CJ*dconjg(Ckx)+2.d0*CBx & +dconjg(Cnu)*CJx+dconjg(CJ)*DCJx-dconjg(Ck)*CJx & +(dconjg(Cnu)*(CJx-CJ*CJ*dconjg(CJx)) & +DCJ*(dconjg(CJx)-dconjg(CJ*CJ)*CJx)/2.d0/KK/KK) & -2.d0*(2.d0*CB+CQ)/xx/(1.d0-xx) return end function Q_rhs function U_rhs(xx,Rmin,beta,KK,CQ,CJ) result(gont) implicit none double complex,intent(in) :: CQ,CJ real*8,intent(in) :: xx,Rmin,beta,KK double complex :: gont #if 1 gont = dexp(2.d0*beta)/Rmin/xx/xx*(KK*CQ-CJ*dconjg(CQ)) #else gont = CQ/Rmin/xx/xx #endif #if 0 if(cdabs(gont)>1)then write(*,*)beta,KK,CQ,CJ stop endif #endif return end function U_rhs function W_rhs(xx,Rmin,beta,KK,DCB,CB,CJ,Cnu,Ck,W, & CQ,bDCk,bDCnu,bDCB,bDCU,bDCUx,DCJ) result(gont) implicit none double complex,intent(in) :: DCB,CB,CJ,Cnu,Ck,CQ,bDCk,bDCnu,bDCB,bDCU,bDCUx,DCJ real*8,intent(in) :: xx,Rmin,beta,KK,W real*8 :: Ric,gont Ric = dreal(2.d0*KK+bDCnu-bDCk+(DCJ*dconjg(DCJ)-Cnu*dconjg(Cnu))/4.d0/KK) gont = dreal(dexp(2.d0*beta)*(Ric/2.d0-KK*(bDCB+CB*dconjg(CB))+dconjg(CJ)*(bDCB+CB*CB) & +(Cnu-Ck)*dconjg(CB))-1.d0+2.d0*Rmin*xx/(1.d0-xx)*(bDCU-W) & +Rmin*xx*xx/2.d0*bDCUx-dexp(2.d0*beta)/4.d0* & (KK*KK-CJ*dconjg(CJ))*(KK*dconjg(CQ)-dconjg(CJ)*CQ)*CQ) gont = gont/Rmin/xx/xx return end function W_rhs function Theta_rhs(xx,Rmin,beta,betax,KK,KKx,CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,W,Wx,CJ,DCJ,CJx,CJxx, & DCJx,bDCB,Cnu,Cnux,Ck,Theta) result(gont) implicit none double complex,intent(in) :: CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,CJ,DCJ,CJx,CJxx,DCJx double complex,intent(in) :: bDCB,Cnu,Cnux,Ck,Theta real*8,intent(in) :: xx,Rmin,beta,betax,KK,KKx,W,Wx double complex :: JH,II,gont real*8 :: V,Vx,Pu II = dcmplx(0.d0,1.d0) V = xx*Rmin/(1.d0-xx)*(1.d0+xx*Rmin/(1.d0-xx)*W) Vx = Rmin/(1.d0-xx)**2+2.d0*xx*Rmin*Rmin/(1.d0-xx)**3*W+xx*xx*Rmin*Rmin/(1.d0-xx)**2*Wx Pu = 2.d0*xx*(1.d0-xx)/KK*dreal(Theta*(dconjg(CJx)*KK-dconjg(CJ)*KKx)) JH = (1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(-KK*DCJ*dconjg(CB)+ & (KK*Cnu+(KK*KK-1.d0)*DCJ-2.d0*KK*Ck)*CB+CJ* & ((2.d0*Ck-Cnu)*dconjg(CB)-2.d0*KK*(bDCB+CB*dconjg(CB))+ & 2.d0*dreal((Cnu-Ck)*dconjg(CB)+dconjg(CJ)*(DCB+CB*CB)))) & +0.5d0*Rmin*xx**3*(1.d0-xx)*dexp(-2.d0*beta)* & ((KK*CUx+CJ*dconjg(CUx))**2- & CJ*dreal(dconjg(CUx)*(KK*CUx+CJ*dconjg(CUx)))) & -0.5d0*(Cnu*(xx*(1.d0-xx)*CUx+2.d0*CU)+DCJ*(xx*(1.d0-xx)*dconjg(CUx)+ & 2.d0*dconjg(CU)))+CJ*II*dimag(xx*(1.d0-xx)*bDCUx+2.d0*bDCU) & -xx*(1.d0-xx)*CJx*dreal(bDCU) & +xx*(1.d0-xx)*(dconjg(CU)*DCJ+CU*Cnu)*II*dimag(CJ*dconjg(CJx)) & -xx*(1.d0-xx)*(dconjg(CU)*DCJx+CU*Cnux) & -2.d0*xx*(1.d0-xx)*(CJ*KKx-KK*CJx)*(dreal(dconjg(CU)*Ck)+ & II*dimag(KK*bDCU-dconjg(CJ)*DCU)) & -8.d0*CJ*((1.d0-xx)**2/Rmin+xx*(1.d0-xx)*W)*betax gont = -KK*(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(DCB+CB*CB) & -(xx*(1.d0-xx)*Wx+W)*CJ+JH+CJ*Pu-2.d0*Theta & -(1.d0-xx)*(1.d0-xx)/xx/xx/Rmin/Rmin*V*(CJ+xx*(1.d0-xx)*CJx) & +(1.d0-xx)*(1.d0-xx)*(1.d0-xx)/xx/Rmin/Rmin*Vx*(CJ+xx*(1.d0-xx)*CJx) & +(1.d0-xx)**4/xx/Rmin/Rmin*V*(2.d0*CJx+xx*CJxx) #if 0 gont = -(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*DCB & -2.d0*Theta & +(1.d0-xx)**3/Rmin*(2.d0*CJx+xx*CJxx) #endif gont = gont/2.d0/xx/(1.d0-xx) return end function Theta_rhs !/////////////////////////////////////////////////////////////////////////////////////////////////////////////// subroutine fake_Theta_rhs(lx,X,rhs,Theta) implicit none integer,intent(in) :: lx double complex,dimension(lx),intent(in) :: Theta double complex,dimension(lx),intent(out) :: rhs real*8,dimension(lx),intent(in) :: X call cderivs_x(lx,X,Theta,rhs) return end subroutine fake_Theta_rhs !/////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! try other guy's old method function Theta_rhs_o(xx,Rmin,beta,betax,KK,KKx,CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,W,Wx,CJ,DCJ,CJx,CJxx, & DCJx,bDCB,Cnu,Cnux,Ck,Theta) result(gont) implicit none double complex,intent(in) :: CU,CUx,DCUx,bDCU,bDCUx,DCU,CB,DCB,CJ,DCJ,CJx,CJxx,DCJx double complex,intent(in) :: Cnu,Cnux,Ck,bDCB,Theta real*8,intent(in) :: xx,Rmin,beta,betax,KK,KKx,W,Wx double complex :: JH,II,gont real*8 :: V,Vx,Pu II = dcmplx(0.d0,1.d0) V = xx*Rmin/(1.d0-xx)*(1.d0+xx*Rmin/(1.d0-xx)*W) Vx = Rmin/(1.d0-xx)**2+2.d0*xx*Rmin*Rmin/(1.d0-xx)**3*W+xx*xx*Rmin*Rmin/(1.d0-xx)**2*Wx Pu = 2.d0*xx*(1.d0-xx)/KK*dreal(Theta*(dconjg(CJx)*KK-dconjg(CJ)*KKx)) JH = (1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(-KK*DCJ*dconjg(CB)+ & (KK*Cnu+(KK*KK-1.d0)*DCJ-2.d0*KK*Ck)*CB+CJ* & ((2.d0*Ck-Cnu)*dconjg(CB)-2.d0*KK*(bDCB+CB*dconjg(CB))+ & 2.d0*dreal((Cnu-Ck)*dconjg(CB)+dconjg(CJ)*(DCB+CB*CB)))) & +0.5d0*Rmin*xx**3*(1.d0-xx)*dexp(-2.d0*beta)* & ((KK*CUx+CJ*dconjg(CUx))**2- & CJ*dreal(dconjg(CUx)*(KK*CUx+CJ*dconjg(CUx)))) & -0.5d0*(Cnu*(xx*(1.d0-xx)*CUx+2.d0*CU)+DCJ*(xx*(1.d0-xx)*dconjg(CUx)+ & 2.d0*dconjg(CU)))+CJ*II*dimag(xx*(1.d0-xx)*bDCUx+2.d0*bDCU) & -xx*(1.d0-xx)*CJx*dreal(bDCU) & +xx*(1.d0-xx)*(dconjg(CU)*DCJ+CU*Cnu)*II*dimag(CJ*dconjg(CJx)) & -xx*(1.d0-xx)*(dconjg(CU)*DCJx+CU*Cnux) & -2.d0*xx*(1.d0-xx)*(CJ*KKx-KK*CJx)*(dreal(dconjg(CU)*Ck)+ & II*dimag(KK*bDCU-dconjg(CJ)*DCU)) & -8.d0*CJ*((1.d0-xx)**2/Rmin+xx*(1.d0-xx)*W)*betax gont = -KK*(xx*(1-xx)*DCUx+2.d0*DCU)+2.d0*(1.d0-xx)/xx/Rmin*dexp(2.d0*beta)*(DCB+CB*CB) & -(xx*(1.d0-xx)*Wx+W)*CJ+JH+CJ*Pu & -(1.d0-xx)*(1.d0-xx)/xx/xx/Rmin/Rmin*V*(CJ+xx*(1.d0-xx)*CJx) & +(1.d0-xx)*(1.d0-xx)*(1.d0-xx)/xx/Rmin/Rmin*Vx*(CJ+xx*(1.d0-xx)*CJx) & +(1.d0-xx)**4/xx/Rmin/Rmin*V*(2.d0*CJx+xx*CJxx) return end function Theta_rhs_o #if (RKorAM == 0) !-------------------------------------------------------------------- ! this R is indeed x function NullEvol_Theta_o(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB real*8,dimension(ex(3)) :: KK,KKx,HKK,HKKx,Hbeta,betax,Hbetax,HW,Wx,HWx double complex :: Theta_rhs_o real*8 :: dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) call cget_half_x(ex(3),CB(i,j,:),HCB) call cget_half_x(ex(3),DCB(i,j,:),HDCB) call cget_half_x(ex(3),bDCB(i,j,:),HbDCB) call cget_half_x(ex(3),Cnu,HCnu) call cderivs_x(ex(3),R,Cnu,Cnux) call cget_half_x(ex(3),Cnux,HCnux) call cget_half_x(ex(3),Ck,HCk) call rget_half_x(ex(3),beta(i,j,:),Hbeta) call rderivs_x(ex(3),R,beta(i,j,:),betax) call rget_half_x(ex(3),betax,Hbetax) KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:)) call rget_half_x(ex(3),KK,HKK) call rderivs_x(ex(3),R,KK,KKx) call rget_half_x(ex(3),KKx,HKKx) call rderivs_x(ex(3),R,W,Wx) call rget_half_x(ex(3),Wx,HWx) call rget_half_x(ex(3),W(i,j,:),HW) call cget_half_x(ex(3),CU(i,j,:),HCU) call cderivs_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_x(ex(3),R,CU(i,j,:),CUx) call cget_half_x(ex(3),DCUx,HDCUx) call cget_half_x(ex(3),CUx,HCUx) call cget_half_x(ex(3),DCU(i,j,:),HDCU) call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx) call cget_half_x(ex(3),bDCUx,HbDCUx) call cget_half_x(ex(3),bDCU(i,j,:),HbDCU) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx) call cget_half_x(ex(3),CJx,HCJx) call cget_half_x(ex(3),CJxx,HCJxx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) call cget_half_x(ex(3),DCJx,HDCJx) do k=1,ex(3)-1 RHS = Theta_rhs_o(R(k)+dR/2.d0,Rmin,Hbeta(k),betax(k),HKK(k),KKx(k),HCU(k),CUx(k),DCUx(k),HbDCU(k),bDCUx(k), & HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k),HCJ(k),HDCJ(k), & CJx(k),CJxx(k),DCJx(k),HbDCB(k),HCnu(k),Cnux(k),HCk(k),CTheta0) CTheta1 = RHS-(1-2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)*CTheta0 CTheta1 = CTheta1/(1+2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR) CTheta0 = CTheta1 RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo enddo enddo gont = 0 return end function NullEvol_Theta_o !------------------------------------------------------------------------------ ! this R is indeed x function NullEvol_beta(ex,crho,sigma,R,RJ,IJ,beta,KKx,HKKx) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KKx,HKKx ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR double complex, dimension(ex(3)):: CJ,CJx,HCJx real*8 :: betah0,betah1,betah,rhs integer :: i,j,k,RK4 real*8 :: beta_rhs !!! sanity check dR = sum(RJ)+sum(IJ)+sum(beta)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_beta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_beta: find NaN in IJ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_beta: find NaN in beta" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_beta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_beta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) do j=1,ex(2) do i=1,ex(1) betah0 = beta(i,j,1) CJ = dcmplx(RJ(i,j,:),IJ(i,j,:)) call cderivs_x(ex(3),R,CJ,CJx) call cget_half_x(ex(3),CJx,HCJx) #ifdef OLD do k = 1,ex(3)-1 ! note our CJx(ex(3)) = (CJ(ex(3))-CJ(ex(3)-1))/dR ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR rhs = beta_rhs(R(k)+dR/2.d0,CJx(k+1),KKx(i,j,k+1)) beta(i,j,k+1) = beta(i,j,k) + rhs*dR enddo #else do k=1,ex(3)-1 RK4 = 0 rhs = beta_rhs(R(k),CJx(k),KKx(i,j,k)) call rungekutta4_scalar(dR,betah0,betah,rhs,RK4) RK4 = 1 betah1 = beta_rhs(R(k)+dR/2.d0,HCJx(k),HKKx(i,j,k)) call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4) call rswap(betah,betah1) RK4 = 2 betah1 = beta_rhs(R(k)+dR/2.d0,HCJx(k),HKKx(i,j,k)) call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4) call rswap(betah,betah1) RK4 = 3 betah1 = beta_rhs(R(k+1),CJx(k+1),KKx(i,j,k+1)) call rungekutta4_scalar(dR,betah0,betah1,rhs,RK4) call rswap(betah0,betah1) beta(i,j,k+1) = betah0 enddo ! above k takes ex(3)-1 then do not need this closing step #if 1 ! closing step k = ex(3)-1 ! note our CJx(ex(3)) = (CJ(ex(3))-CJ(ex(3)-1))/dR ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR rhs = beta_rhs(R(k)+dR/2.d0,CJx(k+1),KKx(i,j,k+1)) beta(i,j,k+1) = beta(i,j,k) + rhs*dR #endif #endif enddo enddo gont = 0 return end function NullEvol_beta !------------------------------------------------------------------------------ ! this R is indeed x function NullEvol_Q(ex,crho,sigma,R,RJ,IJ,Rk,Ik,Rnu,Inu,RB,IB,RQ,IQ,KK,Hkk,KKx,HKKx, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RQ,IQ real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KK,Hkk,KKx,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rk,Ik,Rnu,Inu,RB,IB real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: xx,dR double complex :: CQ0,CQ,CQ1,RHS double complex,dimension(ex(3)) :: CJx,HCJx,DCJx,HDCJx,Ck,Ckx,HCkx,Cnu,Cnux,HCnux,CB,CBx,HCBx double complex,dimension(ex(3)) :: HCJ,HCk,HCnu,HCB,HDCJ double complex, dimension(ex(1),ex(2),ex(3)) :: CJ,DCJ integer :: i,j,k,RK4 double complex :: Q_rhs !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ) & +sum(RK)+sum(IK)+sum(Rnu)+sum(Inu)+sum(RB)+sum(IB) & +sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Q: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Q: find NaN in IJ" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_Q: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_Q: find NaN in IQ" if(sum(RK).ne.sum(RK))write(*,*)"NullEvol_Q: find NaN in RK" if(sum(IK).ne.sum(IK))write(*,*)"NullEvol_Q: find NaN in IK" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Q: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Q: find NaN in Inu" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Q: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Q: find NaN in IB" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Q: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Q: find NaN in HKK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Q: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Q: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CQ0 = dcmplx(RQ(i,j,1),IQ(i,j,1)) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cget_half_x(ex(3),CJx,HCJx) call cget_half_x(ex(3),CJ,HCJ) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) call cget_half_x(ex(3),DCJx,HDCJx) call cget_half_x(ex(3),DCJ,HDCJ) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) call cderivs_x(ex(3),R,Ck,Ckx) call cget_half_x(ex(3),Ckx,HCkx) call cget_half_x(ex(3),Ck,HCk) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) call cderivs_x(ex(3),R,Cnu,Cnux) call cget_half_x(ex(3),Cnux,HCnux) call cget_half_x(ex(3),Cnu,HCnu) CB = dcmplx(RB(i,j,:),IB(i,j,:)) call cderivs_x(ex(3),R,CB,CBx) call cget_half_x(ex(3),CBx,HCBx) call cget_half_x(ex(3),CB,HCB) #ifdef OLD do k = 1,ex(3)-1 xx = R(k)+dR/2.d0 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR RHS = Q_rhs(xx,HCJ(k),CJx(k+1),DCJx(k+1),HKK(i,j,k),HCk(k),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1),CBx(k+1),HCnu(k),HDCJ(k),HCB(k),0) RHS = RHS+CQ0*(1.d0/dR-1.d0/xx/(1.d0-xx)) CQ0 = RHS/(1.d0/dR+1.d0/xx/(1.d0-xx)) RQ(i,j,k+1) = dreal(CQ0) IQ(i,j,k+1) = dimag(CQ0) enddo #else do k=1,ex(3)-2 RK4 = 0 RHS = Q_rhs(R(k),CJ(i,j,k),CJx(k),DCJx(k),KK(i,j,k),Ck(k),Ckx(k),Cnux(k),KKx(i,j,k),CBx(k),Cnu(k),DCJ(i,j,k),CB(k),CQ0) call rungekutta4_cplxscalar(dR,CQ0,CQ,RHS,RK4) RK4 = 1 CQ1 = Q_rhs(R(k)+dR/2.d0,HCJ(k),HCJx(k),HDCJx(k),HKK(i,j,k),HCk(k),HCkx(k),HCnux(k),HKKx(i,j,k), & HCBx(k),HCnu(k),HDCJ(k),HCB(k),CQ) call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4) call cswap(CQ,CQ1) RK4 = 2 CQ1 = Q_rhs(R(k)+dR/2.d0,HCJ(k),HCJx(k),HDCJx(k),HKK(i,j,k),HCk(k),HCkx(k),HCnux(k),HKKx(i,j,k), & HCBx(k),HCnu(k),HDCJ(k),HCB(k),CQ) call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4) call cswap(CQ,CQ1) RK4 = 3 CQ1 = Q_rhs(R(k+1),CJ(i,j,k+1),CJx(k+1),DCJx(k+1),KK(i,j,k+1),Ck(k+1),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1), & CBx(k+1),Cnu(k+1),DCJ(i,j,k+1),CB(k+1),CQ) call rungekutta4_cplxscalar(dR,CQ0,CQ1,RHS,RK4) call cswap(CQ0,CQ1) RQ(i,j,k+1) = dreal(CQ0) IQ(i,j,k+1) = dimag(CQ0) enddo #if 0 k = ex(3) CQ0 = -2*CB(k) RQ(i,j,k+1) = dreal(CQ0) IQ(i,j,k+1) = dimag(CQ0) #else ! closing step k = ex(3)-1 CQ0 = dcmplx(RQ(i,j,k),IQ(i,j,k)) xx = R(k)+dR/2.d0 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR RHS = Q_rhs(xx,HCJ(k),CJx(k+1),DCJx(k+1),HKK(i,j,k), & HCk(k),Ckx(k+1),Cnux(k+1),KKx(i,j,k+1),CBx(k+1),HCnu(k),HDCJ(k),HCB(k),dcmplx(0.d0,0.d0)) RHS = RHS+CQ0*(1.d0/dR-1.d0/xx/(1.d0-xx)) CQ0 = RHS/(1.d0/dR+1.d0/xx/(1.d0-xx)) RQ(i,j,k+1) = dreal(CQ0) IQ(i,j,k+1) = dimag(CQ0) #endif #endif enddo enddo gont = 0 return end function NullEvol_Q !-------------------------------------------------------------------- ! this R is indeed x function NullEvol_U(ex,crho,sigma,R,RJ,IJ,RQ,IQ,KK,HKK,beta,RU,IU, & Rmin) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RQ,IQ,beta,KK,HKK real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RU,IU real*8,intent(in) :: Rmin ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR double complex :: CU0,CU,CU1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: CJ,CQ,HCJ,HCQ real*8,dimension(ex(3)) :: Hbeta double complex :: U_rhs !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ)+sum(beta)+sum(RU)+sum(IU)+sum(KK)+sum(HKK) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_U: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_U: find NaN in IJ" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_U: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_U: find NaN in IQ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_U: find NaN in beta" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_U: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_U: find NaN in IU" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_U: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_U: find NaN in HKK" gont = 1 return endif dR = R(2) - R(1) do j=1,ex(2) do i=1,ex(1) CU0 = dcmplx(RU(i,j,1),IU(i,j,1)) CJ = dcmplx(RJ(i,j,:),IJ(i,j,:)) CQ = dcmplx(RQ(i,j,:),IQ(i,j,:)) call cget_half_x(ex(3),CJ,HCJ) call cget_half_x(ex(3),CQ,HCQ) call rget_half_x(ex(3),beta(i,j,:),Hbeta) #ifdef OLD do k = 1,ex(3)-1 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR RHS = U_rhs(R(k)+dR/2,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k)) CU0 = CU0+RHS*dR RU(i,j,k+1) = dreal(CU0) IU(i,j,k+1) = dimag(CU0) enddo #else do k=1,ex(3)-2 RK4 = 0 RHS = U_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),CQ(k),CJ(k)) call rungekutta4_cplxscalar(dR,CU0,CU,RHS,RK4) RK4 = 1 CU1 = U_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k)) call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4) call cswap(CU,CU1) RK4 = 2 CU1 = U_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k)) call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4) call cswap(CU,CU1) RK4 = 3 CU1 = U_rhs(R(k+1),Rmin,beta(i,j,k+1),KK(i,j,k+1),CQ(k+1),CJ(k+1)) call rungekutta4_cplxscalar(dR,CU0,CU1,RHS,RK4) call cswap(CU0,CU1) RU(i,j,k+1) = dreal(CU0) IU(i,j,k+1) = dimag(CU0) enddo ! above k takes ex(3)-1 then do not need closing step #if 1 ! closing step k = ex(3)-1 CU0 = dcmplx(RU(i,j,k),IU(i,j,k)) ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR RHS = U_rhs(R(k)+dR/2,Rmin,Hbeta(k),HKK(i,j,k),HCQ(k),HCJ(k)) CU0 = CU0+RHS*dR RU(i,j,k+1) = dreal(CU0) IU(i,j,k+1) = dimag(CU0) #endif #endif enddo enddo gont = 0 return end function NullEvol_U !---------------------------------------------------------------------------------------- ! this R is indeed x function NullEvol_W(ex,crho,sigma,R,RJ,IJ,RB,IB,Rnu,Inu,Rk,Ik, & RU,IU,RQ,IQ,W,beta,KK,HKK,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: W real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,RB,IB real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RU,IU,RQ,IQ,beta,KK,HKK real*8,intent(in ) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR real*8, dimension(ex(3)) :: Hbeta double complex, dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU double complex, dimension(ex(1),ex(2),ex(3)) :: CB,DCB,bDCB,CJ,DCJ,Cnu,bDCnu,Ck,bDCk double complex, dimension(ex(3)) :: HCB,HDCB,HbDCB,HCJ,HDCJ,HCnu,HbDCnu,HCk,HbDCk double complex, dimension(ex(3)) :: HbDCU,bDCUx,HbDCUx,CQ,HCQ real*8 :: Wh0,Wh1,Wh,rhs integer :: i,j,k,RK4 real*8 :: xx,W_rhs !!! sanity check dR = sum(RJ)+sum(IJ)+sum(beta)+sum(RB)+sum(IB)+sum(Rnu)+sum(Inu) & +sum(Rk)+sum(Ik)+sum(W)+sum(RU)+sum(IU)+sum(RQ)+sum(IQ)& +sum(KK)+sum(HKK) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_W: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_W: find NaN in IJ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_W: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_W: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_W: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_W: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_W: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_W: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_W: find NaN in Ik" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_W: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_W: find NaN in IU" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_W: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_W: find NaN in IQ" if(sum(W).ne.sum(W))write(*,*)"NullEvol_W: find NaN in W" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_W: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_W: find NaN in HKK" gont = 1 return endif dR = R(2) - R(1) CB = dcmplx(RB,IB) CU = dcmplx(RU,IU) Ck = dcmplx(Rk,Ik) Cnu = dcmplx(Rnu,Inu) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,Ck(:,:,k),bDCk(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,Cnu(:,:,k),bDCnu(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) Wh0 = W(i,j,1) call rget_half_x(ex(3),beta(i,j,:),Hbeta) call cderivs_x(ex(3),R,bDCU,bDCUx) call cget_half_x(ex(3),bDCUx,HbDCUx) call cget_half_x(ex(3),bDCU,HbDCU) call cget_half_x(ex(3),DCJ,HDCJ) call cget_half_x(ex(3),DCB,HDCB) call cget_half_x(ex(3),bDCB,HbDCB) call cget_half_x(ex(3),CB,HCB) call cget_half_x(ex(3),CJ,HCJ) call cget_half_x(ex(3),Cnu,HCnu) call cget_half_x(ex(3),Ck,HCk) CQ = dcmplx(RQ(i,j,:),IQ(i,j,:)) call cget_half_x(ex(3),CQ,HCQ) call cget_half_x(ex(3),bDCk,HbDCk) call cget_half_x(ex(3),bDCnu,HbDCnu) #ifdef OLD do k = 1,ex(3)-1 xx = R(k)+dR/2 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR rhs = W_rhs(xx,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),0, & HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),bDCUx(k+1),HDCJ(k)) rhs = rhs+Wh0*(1.d0/dR-1.d0/xx/(1.d0-xx)) W(i,j,k+1) = rhs/(1.d0/dR+1.d0/xx/(1.d0-xx)) enddo #else do k=1,ex(3)-2 RK4 = 0 rhs = W_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),DCB(i,j,k),CB(i,j,k),CJ(i,j,k),Cnu(i,j,k),Ck(i,j,k),Wh0, & CQ(k),bDCk(i,j,k),bDCnu(i,j,k),bDCB(i,j,k),bDCU(i,j,k),bDCUx(k),DCJ(i,j,k)) call rungekutta4_scalar(dR,Wh0,Wh,rhs,RK4) RK4 = 1 Wh1 = W_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),Wh, & HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),HbDCUx(k),HDCJ(k)) call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4) call rswap(Wh,Wh1) RK4 = 2 Wh1 = W_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),Wh, & HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),HbDCUx(k),HDCJ(k)) call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4) call rswap(Wh,Wh1) RK4 = 3 Wh1 = W_rhs(R(k+1),Rmin,beta(i,j,k+1),KK(i,j,k+1),DCB(i,j,k+1),CB(i,j,k+1),CJ(i,j,k+1),Cnu(i,j,k+1),Ck(i,j,k+1),Wh, & CQ(k+1),bDCk(i,j,k+1),bDCnu(i,j,k+1),bDCB(i,j,k+1),bDCU(i,j,k+1),bDCUx(k+1),DCJ(i,j,k+1)) call rungekutta4_scalar(dR,Wh0,Wh1,rhs,RK4) call rswap(Wh0,Wh1) W(i,j,k+1) = Wh0 enddo #if 0 k = ex(3) W(i,j,k) = dreal(bDCU(i,j,k)) #else ! closing step k = ex(3)-1 Wh0 = W(i,j,k) xx = R(k)+dR/2 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR rhs = W_rhs(xx,Rmin,Hbeta(k),HKK(i,j,k),HDCB(k),HCB(k),HCJ(k),HCnu(k),HCk(k),0.d0, & HCQ(k),HbDCk(k),HbDCnu(k),HbDCB(k),HbDCU(k),bDCUx(k+1),HDCJ(k)) rhs = rhs+Wh0*(1.d0/dR-1.d0/xx/(1.d0-xx)) W(i,j,k+1) = rhs/(1.d0/dR+1.d0/xx/(1.d0-xx)) #endif #endif enddo enddo gont = 0 return end function NullEvol_W !----------------------------------------------------------------------------------------------- ! given exact Theta_x ! this R is indeed x function NullEvol_Theta_givenx(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI,T,sst) result(gont) implicit none integer,intent(in ):: ex(1:3),sst real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin,T real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2 real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2 real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI ! gont = 0: success; gont = 1: something wrong integer::gont real*8,dimension(ex(3))::HR real*8,dimension(ex(1),ex(2),ex(3)) :: RThetax,IThetax,HRThetax,HIThetax double complex,dimension(ex(3)) :: fRHS,HfRHS real*8 :: xx,dR integer :: i,j,k,RK4 double complex :: CTheta0,CTheta,CTheta1,RHS integer,parameter :: ks=1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + & sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) HR = R+dR/2 call get_exact_null_theta_x(ex,crho,sigma,R,RThetax,IThetax,sst,Rmin,T, & quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, & gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI) call get_exact_null_theta_x(ex,crho,sigma,HR,HRThetax,HIThetax,sst,Rmin,T, & quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, & gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI) do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,ks),ITheta(i,j,ks)) fRHS = dcmplx(RThetax(i,j,:),IThetax(i,j,:)) HfRHS = dcmplx(HRThetax(i,j,:),HIThetax(i,j,:)) ! call cget_half_x(ex(3),fRHS,HfRHS) do k=ks,ex(3)-1 RK4 = 0 RHS = fRHS(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4) RK4 = 1 CTheta1 = HfRHS(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 2 CTheta1 = HfRHS(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 3 CTheta1 = fRHS(k+1) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta0,CTheta1) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo #if 0 ! closing step k = ex(3)-1 RHS = fRHS(k) CTheta0 = dcmplx(RTheta(i,j,k),ITheta(i,j,k))+RHS*dR RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) #endif enddo enddo gont = 0 return end function NullEvol_Theta_givenx !----------------------------------------------------------------------------------------------- #if 1 ! real evolve ! for eth_x, eth first, _x later ! this R is indeed x function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB real*8,dimension(ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx double complex :: Theta_rhs,Theta_rhs_o real*8 :: xx,dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + & sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) call cget_half_x(ex(3),CB(i,j,:),HCB) call cget_half_x(ex(3),DCB(i,j,:),HDCB) call cget_half_x(ex(3),bDCB(i,j,:),HbDCB) call cget_half_x(ex(3),Cnu,HCnu) call cderivs_x(ex(3),R,Cnu,Cnux) call cget_half_x(ex(3),Cnux,HCnux) call cget_half_x(ex(3),Ck,HCk) call rget_half_x(ex(3),beta(i,j,:),Hbeta) call rderivs_x(ex(3),R,beta(i,j,:),betax) call rget_half_x(ex(3),betax,Hbetax) call rderivs_x(ex(3),R,W,Wx) call rget_half_x(ex(3),Wx,HWx) call rget_half_x(ex(3),W(i,j,:),HW) call cget_half_x(ex(3),CU(i,j,:),HCU) call cderivs_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_x(ex(3),R,CU(i,j,:),CUx) call cget_half_x(ex(3),DCUx,HDCUx) call cget_half_x(ex(3),CUx,HCUx) call cget_half_x(ex(3),DCU(i,j,:),HDCU) call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx) call cget_half_x(ex(3),bDCUx,HbDCUx) call cget_half_x(ex(3),bDCU(i,j,:),HbDCU) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx) call cget_half_x(ex(3),CJx,HCJx) call cget_half_x(ex(3),CJxx,HCJxx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) call cget_half_x(ex(3),DCJx,HDCJx) ! old type code: PRD 54, 6153, Eq.(32) etc. #if 0 ! start up part k = 1 RHS = Theta_rhs_o(R(k)+dR/2.d0,Rmin,Hbeta(k),betax(k),HKK(i,j,k),KKx(i,j,k),HCU(k),CUx(k),DCUx(k),HbDCU(k),bDCUx(k), & HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k),HCJ(k),HDCJ(k), & CJx(k),CJxx(k),DCJx(k),HbDCB(k),HCnu(k),Cnux(k),HCk(k),CTheta0) CTheta1 = RHS-(1-2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR)*CTheta0 CTheta0 = CTheta1/(1+2.d0*(R(k)+dR/2.d0)*(1.d0-R(k)-dR/2.d0)/dR) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) do k=1,ex(3)-2 RHS = Theta_rhs_o(R(k+1),Rmin,beta(i,j,k+1),betax(k+1),KK(i,j,k+1),KKx(i,j,k+1),CU(i,j,k+1),CUx(k+1),DCUx(k+1),bDCU(i,j,k+1),bDCUx(k+1), & DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(k+1),CJ(i,j,k+1),DCJ(i,j,k+1), & CJx(k+1),CJxx(k+1),DCJx(k+1),bDCB(i,j,k+1),Cnu(k+1),Cnux(k+1),Ck(k+1),CTheta0) CTheta1 = RHS-(1-R(k+1)*(1.d0-R(k+1))/dR)*(dcmplx(RTheta(i,j,k),ITheta(i,j,k))) CTheta0 = CTheta1/(1+R(k+1)*(1.d0-R(k+1))/dR) RTheta(i,j,k+2) = dreal(CTheta0) ITheta(i,j,k+2) = dimag(CTheta0) enddo #endif #ifdef OLD do k = 1,ex(3)-1 xx = R(k)+dR/2 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR ! note our fxx(ex(3)) = (f(ex(3))-2.d0*f(ex(3)-1)+f(ex(3)-2))/dR RHS = Theta_rhs(xx,Rmin,Hbeta(k),betax(k+1),HKK(i,j,k),KKx(i,j,k+1),HCU(k),CUx(k+1),DCUx(k+1),HbDCU(k),bDCUx(k+1), & HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k+1),HCJ(k),HDCJ(k), & CJx(k+1),CJxx(k+1),DCJx(k+1),HbDCB(k),HCnu(k),Cnux(k+1),HCk(k),0) RHS = RHS+CTheta0*(1.d0/dR-0.5d0/xx/(1.d0-xx)) CTheta0 = RHS/(1.d0/dR+0.5d0/xx/(1.d0-xx)) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo #else do k=1,ex(3)-2 RK4 = 0 RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), & DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), & CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),CTheta0) call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4) RK4 = 1 CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),Hbetax(k),HKK(i,j,k),HKKx(i,j,k), & HCU(k),HCUx(k),HDCUx(k),HbDCU(k),HbDCUx(k), & HDCU(k),HCB(k),HDCB(k),HW(k),HWx(k),HCJ(k),HDCJ(k), & HCJx(k),HCJxx(k),HDCJx(k),HbDCB(k),HCnu(k),HCnux(k),HCk(k),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 2 CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(k),Hbetax(k),HKK(i,j,k),HKKx(i,j,k), & HCU(k),HCUx(k),HDCUx(k),HbDCU(k),HbDCUx(k), & HDCU(k),HCB(k),HDCB(k),HW(k),HWx(k),HCJ(k),HDCJ(k), & HCJx(k),HCJxx(k),HDCJx(k),HbDCB(k),HCnu(k),HCnux(k),HCk(k),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 3 CTheta1 = Theta_rhs(R(k+1),Rmin,beta(i,j,k+1),betax(k+1),KK(i,j,k+1),KKx(i,j,k+1), & CU(i,j,k+1),CUx(k+1),DCUx(k+1),bDCU(i,j,k+1),bDCUx(k+1), & DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(k+1),CJ(i,j,k+1),DCJ(i,j,k+1), & CJx(k+1),CJxx(k+1),DCJx(k+1),bDCB(i,j,k+1),Cnu(k+1),Cnux(k+1),Ck(k+1),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta0,CTheta1) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo #if 0 k = ex(3) CTheta0 = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 & +CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2 RTheta(i,j,k) = dreal(CTheta0) ITheta(i,j,k) = dimag(CTheta0) #else ! closing step k = ex(3)-1 CTheta0 = dcmplx(RTheta(i,j,k),ITheta(i,j,k)) xx = R(k)+dR/2 ! note our HCJ(ex(3)-1) = (CJ(ex(3))+CJ(ex(3)-1))/2 ! note our KKx(ex(3)) = (KK(ex(3))-KK(ex(3)-1))/dR ! note our fxx(ex(3)) = (f(ex(3))-2.d0*f(ex(3)-1)+f(ex(3)-2))/dR RHS = Theta_rhs(xx,Rmin,Hbeta(k),betax(k+1),HKK(i,j,k),KKx(i,j,k+1),HCU(k),CUx(k+1),DCUx(k+1),HbDCU(k),bDCUx(k+1), & HDCU(k),HCB(k),HDCB(k),HW(k),Wx(k+1),HCJ(k),HDCJ(k), & CJx(k+1),CJxx(k+1),DCJx(k+1),HbDCB(k),HCnu(k),Cnux(k+1),HCk(k),dcmplx(0.d0,0.d0)) RHS = RHS+CTheta0*(1.d0/dR-0.5d0/xx/(1.d0-xx)) CTheta0 = RHS/(1.d0/dR+0.5d0/xx/(1.d0-xx)) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) #endif #endif enddo enddo gont = 0 return end function NullEvol_Theta !-------------------------------------------------------------------- ! check with fake_Theta_rhs #elif 0 ! this R is indeed x function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 integer,parameter :: ks=1 double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB real*8,dimension(ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx double complex :: Theta_rhs,Theta_rhs_o real*8 :: xx,dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + & sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,ks),ITheta(i,j,ks)) Cnu = dcmplx(RTheta(i,j,:),ITheta(i,j,:)) call fake_Theta_rhs(ex(3),R,Ck,Cnu) call cget_half_x(ex(3),Ck,HCk) do k=ks,ex(3)-1 RK4 = 0 RHS = Ck(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4) RK4 = 1 CTheta1 = HCk(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 2 CTheta1 = HCk(k) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 3 CTheta1 = Ck(k+1) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta0,CTheta1) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo enddo enddo gont = 0 return end function NullEvol_Theta #else ! for eth_x, _x first, eth second ! this R is indeed x function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,HKK,KKx,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(1),ex(2),ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx double complex,dimension(ex(1),ex(2),ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB real*8,dimension(ex(1),ex(2),ex(3)) :: Hbeta,betax,Hbetax,HW,Wx,HWx double complex :: Theta_rhs real*8 :: xx,dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + & sum(KK)+sum(HKK)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK" if(sum(HKK).ne.sum(HKK))write(*,*)"NullEvol_Theta: find NaN in HKK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_Theta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) Cnu = dcmplx(Rnu,Inu) Ck = dcmplx(Rk,Ik) do j=1,ex(2) do i=1,ex(1) call cderivs_x(ex(3),R,Cnu(i,j,:),Cnux(i,j,:)) call rderivs_x(ex(3),R,beta(i,j,:),betax(i,j,:)) call rderivs_x(ex(3),R,W(i,j,:),Wx(i,j,:)) call cderivs_x(ex(3),R,CU(i,j,:),CUx(i,j,:)) call cderivs_x(ex(3),R,CJ(i,j,:),CJx(i,j,:)) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx(i,j,:)) enddo enddo do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CUx(:,:,k),DCUx(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CUx(:,:,k),bDCUx(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJx(:,:,k),DCJx(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) call cget_half_x(ex(3),CB(i,j,:),HCB(i,j,:)) call cget_half_x(ex(3),DCB(i,j,:),HDCB(i,j,:)) call cget_half_x(ex(3),bDCB(i,j,:),HbDCB(i,j,:)) call cget_half_x(ex(3),Cnu(i,j,:),HCnu(i,j,:)) call cget_half_x(ex(3),Cnux(i,j,:),HCnux(i,j,:)) call cget_half_x(ex(3),Ck(i,j,:),HCk(i,j,:)) call rget_half_x(ex(3),beta(i,j,:),Hbeta(i,j,:)) call rget_half_x(ex(3),betax(i,j,:),Hbetax(i,j,:)) call rget_half_x(ex(3),Wx(i,j,:),HWx(i,j,:)) call rget_half_x(ex(3),W(i,j,:),HW(i,j,:)) call cget_half_x(ex(3),CU(i,j,:),HCU(i,j,:)) call cget_half_x(ex(3),DCUx(i,j,:),HDCUx(i,j,:)) call cget_half_x(ex(3),CUx(i,j,:),HCUx(i,j,:)) call cget_half_x(ex(3),DCU(i,j,:),HDCU(i,j,:)) call cget_half_x(ex(3),bDCUx(i,j,:),HbDCUx(i,j,:)) call cget_half_x(ex(3),bDCU(i,j,:),HbDCU(i,j,:)) call cget_half_x(ex(3),CJx(i,j,:),HCJx(i,j,:)) call cget_half_x(ex(3),CJxx(i,j,:),HCJxx(i,j,:)) call cget_half_x(ex(3),DCJx(i,j,:),HDCJx(i,j,:)) enddo enddo do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) do k=1,ex(3)-2 RK4 = 0 RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(i,j,k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(i,j,k),DCUx(i,j,k),bDCU(i,j,k),bDCUx(i,j,k), & DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(i,j,k),CJ(i,j,k),DCJ(i,j,k), & CJx(i,j,k),CJxx(i,j,k),DCJx(i,j,k),bDCB(i,j,k),Cnu(i,j,k),Cnux(i,j,k),Ck(i,j,k),CTheta0) call rungekutta4_cplxscalar(dR,CTheta0,CTheta,RHS,RK4) RK4 = 1 CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(i,j,k),Hbetax(i,j,k),HKK(i,j,k),HKKx(i,j,k), & HCU(i,j,k),HCUx(i,j,k),HDCUx(i,j,k),HbDCU(i,j,k),HbDCUx(i,j,k), & HDCU(i,j,k),HCB(i,j,k),HDCB(i,j,k),HW(i,j,k),HWx(i,j,k),HCJ(i,j,k),HDCJ(i,j,k), & HCJx(i,j,k),HCJxx(i,j,k),HDCJx(i,j,k),HbDCB(i,j,k),HCnu(i,j,k),HCnux(i,j,k),HCk(i,j,k),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 2 CTheta1 = Theta_rhs(R(k)+dR/2.d0,Rmin,Hbeta(i,j,k),Hbetax(i,j,k),HKK(i,j,k),HKKx(i,j,k), & HCU(i,j,k),HCUx(i,j,k),HDCUx(i,j,k),HbDCU(i,j,k),HbDCUx(i,j,k), & HDCU(i,j,k),HCB(i,j,k),HDCB(i,j,k),HW(i,j,k),HWx(i,j,k),HCJ(i,j,k),HDCJ(i,j,k), & HCJx(i,j,k),HCJxx(i,j,k),HDCJx(i,j,k),HbDCB(i,j,k),HCnu(i,j,k),HCnux(i,j,k),HCk(i,j,k),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta,CTheta1) RK4 = 3 CTheta1 = Theta_rhs(R(k+1),Rmin,beta(i,j,k+1),betax(i,j,k+1),KK(i,j,k+1), & KKx(i,j,k+1),CU(i,j,k+1),CUx(i,j,k+1),DCUx(i,j,k+1),bDCU(i,j,k+1),bDCUx(i,j,k+1), & DCU(i,j,k+1),CB(i,j,k+1),DCB(i,j,k+1),W(i,j,k+1),Wx(i,j,k+1),CJ(i,j,k+1),DCJ(i,j,k+1), & CJx(i,j,k+1),CJxx(i,j,k+1),DCJx(i,j,k+1),bDCB(i,j,k+1),Cnu(i,j,k+1),Cnux(i,j,k+1),Ck(i,j,k+1),CTheta) call rungekutta4_cplxscalar(dR,CTheta0,CTheta1,RHS,RK4) call cswap(CTheta0,CTheta1) RTheta(i,j,k+1) = dreal(CTheta0) ITheta(i,j,k+1) = dimag(CTheta0) enddo k = ex(3) CTheta0 = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(i,j,k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 & +CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2 RTheta(i,j,k) = dreal(CTheta0) ITheta(i,j,k) = dimag(CTheta0) enddo enddo gont = 0 return end function NullEvol_Theta #endif #elif (RKorAM == 1) !------------------------------------------------------------------------------ ! this R is indeed x function NullEvol_beta(ex,crho,sigma,R,RJ,IJ,beta,KKx,HKKx) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KKx,HKKx ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR,beta_rhs double complex, dimension(ex(3)):: CJ,CJx real*8, dimension(ex(3)) :: rhs integer :: i,j,k real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1 real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(beta)+sum(KKx)+sum(HKKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_beta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_beta: find NaN in IJ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_beta: find NaN in beta" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_beta: find NaN in KKx" if(sum(HKKx).ne.sum(HKKx))write(*,*)"NullEvol_beta: find NaN in HKKx" gont = 1 return endif dR = R(2) - R(1) do j=1,ex(2) do i=1,ex(1) CJ = dcmplx(RJ(i,j,:),IJ(i,j,:)) #if 0 call cderivs_sw_x(ex(3),R,CJ,CJx) #else call cderivs_x(ex(3),R,CJ,CJx) #endif do k=1,ex(3) rhs(k) = beta_rhs(R(k),CJx(k),KKx(i,j,k)) enddo k = 1 beta(i,j,k+1) = beta(i,j,k) + (rhs(k+1)+rhs(k))*dR/2 k = 2 beta(i,j,k+1) = beta(i,j,k) + (F5o12*rhs(k+1) + F2o3*rhs(k) - F1o12*rhs(k-1))*dR do k=3,ex(3)-1 beta(i,j,k+1) = beta(i,j,k) + (F3o8*rhs(k+1) + F19o24*rhs(k) - F5o24*rhs(k-1) + F1o24*rhs(k-2))*dR enddo enddo enddo gont = 0 return end function NullEvol_beta !------------------------------------------------------------------------------ ! this R is indeed x function NullEvol_Q(ex,crho,sigma,R,RJ,IJ,Rk,Ik,Rnu,Inu,RB,IB,RQ,IQ,KK,Hkk,KKx,HKKx, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RQ,IQ real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,KK,KKx,HKK,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rk,Ik,Rnu,Inu,RB,IB real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: xx,dR double complex,dimension(ex(3)) :: CQ,RHS real*8, dimension(ex(3)) :: gunc double complex,dimension(ex(3)) :: CJx,DCJx,Ck,Ckx,Cnu,Cnux,CB,CBx double complex, dimension(ex(1),ex(2),ex(3)) :: CJ,DCJ integer :: i,j,k double complex :: ZEO,Q_rhs real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1 real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ) & +sum(RK)+sum(IK)+sum(Rnu)+sum(Inu)+sum(RB)+sum(IB) & +sum(KK)+sum(KKx) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Q: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Q: find NaN in IJ" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_Q: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_Q: find NaN in IQ" if(sum(RK).ne.sum(RK))write(*,*)"NullEvol_Q: find NaN in RK" if(sum(IK).ne.sum(IK))write(*,*)"NullEvol_Q: find NaN in IK" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Q: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Q: find NaN in Inu" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Q: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Q: find NaN in IB" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Q: find NaN in KK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Q: find NaN in KKx" gont = 1 return endif dR = R(2) - R(1) ZEO = dcmplx(0.d0,0.d0) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CQ(1) = dcmplx(RQ(i,j,1),IQ(i,j,1)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) CB = dcmplx(RB(i,j,:),IB(i,j,:)) #if 0 call cderivs_sw_x(ex(3),R,CJ(i,j,:),CJx) call cderivs_sw_x(ex(3),R,DCJ(i,j,:),DCJx) call cderivs_sw_x(ex(3),R,Ck,Ckx) call cderivs_sw_x(ex(3),R,Cnu,Cnux) call cderivs_sw_x(ex(3),R,CB,CBx) #else call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) call cderivs_x(ex(3),R,Ck,Ckx) call cderivs_x(ex(3),R,Cnu,Cnux) call cderivs_x(ex(3),R,CB,CBx) #endif do k = 1,ex(3) RHS(k) = Q_rhs(R(k),CJ(i,j,k),CJx(k),DCJx(k),KK(i,j,k),Ck(k),Ckx(k),Cnux(k),KKx(i,j,k),CBx(k),Cnu(k),DCJ(i,j,k),CB(k),ZEO) gunc(k) = -2/R(k)/(1-R(k)) enddo k = 1 CQ(k+1) = CQ(k) + (RHS(k+1)+RHS(k)+CQ(k)*gunc(k))*dR/2 CQ(k+1) = CQ(k+1)/(1-0.5*dR*gunc(k+1)) k = 2 CQ(k+1) = CQ(k) + (F5o12*RHS(k+1) + F2o3*(RHS(k)+CQ(k)*gunc(k)) - F1o12*(RHS(k-1)+CQ(k-1)*gunc(k-1)))*dR CQ(k+1) = CQ(k+1)/(1-F5o12*dR*gunc(k+1)) do k=3,ex(3)-2 CQ(k+1) = CQ(k) + (F3o8*RHS(k+1) + F19o24*(RHS(k)+CQ(k)*gunc(k)) - F5o24*(RHS(k-1)+CQ(k-1)*gunc(k-1)) & + F1o24*(RHS(k-2)+CQ(k-2)*gunc(k-2)))*dR CQ(k+1) = CQ(k+1)/(1-F3o8*dR*gunc(k+1)) enddo k = ex(3) CQ(k) = -2*CB(k) RQ(i,j,:) = dreal(CQ) IQ(i,j,:) = dimag(CQ) enddo enddo gont = 0 return end function NullEvol_Q !-------------------------------------------------------------------- ! this R is indeed x function NullEvol_U(ex,crho,sigma,R,RJ,IJ,RQ,IQ,KK,HKK,beta,RU,IU, & Rmin) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RQ,IQ,beta,KK,HKK real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RU,IU real*8,intent(in) :: Rmin ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR double complex,dimension(ex(3)) :: CU0,RHS integer :: i,j,k double complex :: U_rhs double complex,dimension(ex(3)) :: CJ,CQ real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1 real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RQ)+sum(IQ)+sum(beta)+sum(RU)+sum(IU)+sum(KK) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_U: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_U: find NaN in IJ" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_U: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_U: find NaN in IQ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_U: find NaN in beta" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_U: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_U: find NaN in IU" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_U: find NaN in KK" gont = 1 return endif dR = R(2) - R(1) do j=1,ex(2) do i=1,ex(1) CU0(1) = dcmplx(RU(i,j,1),IU(i,j,1)) CJ = dcmplx(RJ(i,j,:),IJ(i,j,:)) CQ = dcmplx(RQ(i,j,:),IQ(i,j,:)) do k = 1,ex(3) RHS(k) = U_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),CQ(k),CJ(k)) enddo k = 1 CU0(k+1) = CU0(k) + (RHS(k+1)+RHS(k))*dR/2 k = 2 CU0(k+1) = CU0(k) + (F5o12*RHS(k+1) + F2o3*RHS(k) - F1o12*RHS(k-1))*dR do k=3,ex(3)-1 CU0(k+1) = CU0(k) + (F3o8*RHS(k+1) + F19o24*RHS(k) - F5o24*RHS(k-1) & + F1o24*RHS(k-2))*dR enddo RU(i,j,:) = dreal(CU0) IU(i,j,:) = dimag(CU0) enddo enddo gont = 0 return end function NullEvol_U !---------------------------------------------------------------------------------------- ! this R is indeed x function NullEvol_W(ex,crho,sigma,R,RJ,IJ,RB,IB,Rnu,Inu,Rk,Ik, & RU,IU,RQ,IQ,W,beta,KK,HKK,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: W real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RJ,IJ,RB,IB real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: RU,IU,RQ,IQ,beta,KK,HKK real*8,intent(in ) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont real*8 :: dR double complex, dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU double complex, dimension(ex(1),ex(2),ex(3)) :: CB,DCB,bDCB,CJ,DCJ,Cnu,bDCnu,Ck,bDCk double complex, dimension(ex(3)) :: bDCUx,CQ integer :: i,j,k real*8, dimension(ex(3)) :: rhs,gunc real*8 :: zeo,W_rhs real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1 real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(beta)+sum(RB)+sum(IB)+sum(Rnu)+sum(Inu) & +sum(Rk)+sum(Ik)+sum(W)+sum(RU)+sum(IU)+sum(RQ)+sum(IQ)& +sum(KK) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_W: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_W: find NaN in IJ" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_W: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_W: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_W: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_W: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_W: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_W: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_W: find NaN in Ik" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_W: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_W: find NaN in IU" if(sum(RQ).ne.sum(RQ))write(*,*)"NullEvol_W: find NaN in RQ" if(sum(IQ).ne.sum(IQ))write(*,*)"NullEvol_W: find NaN in IQ" if(sum(W).ne.sum(W))write(*,*)"NullEvol_W: find NaN in W" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_W: find NaN in KK" gont = 1 return endif dR = R(2) - R(1) zeo = 0.d0 CB = dcmplx(RB,IB) CU = dcmplx(RU,IU) Ck = dcmplx(Rk,Ik) Cnu = dcmplx(Rnu,Inu) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,Ck(:,:,k),bDCk(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,Cnu(:,:,k),bDCnu(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) #if 0 call cderivs_sw_x(ex(3),R,bDCU,bDCUx) #else call cderivs_x(ex(3),R,bDCU,bDCUx) #endif CQ = dcmplx(RQ(i,j,:),IQ(i,j,:)) do k = 1,ex(3) rhs(k) = W_rhs(R(k),Rmin,beta(i,j,k),KK(i,j,k),DCB(i,j,k),CB(i,j,k),CJ(i,j,k),Cnu(i,j,k),Ck(i,j,k),zeo, & CQ(k),bDCk(i,j,k),bDCnu(i,j,k),bDCB(i,j,k),bDCU(i,j,k),bDCUx(k),DCJ(i,j,k)) gunc(k) = -2/R(k)/(1-R(k)) enddo k = 1 W(i,j,k+1) = W(i,j,k) + (rhs(k+1)+rhs(k)+W(i,j,k)*gunc(k))*dR/2 W(i,j,k+1) = W(i,j,k+1)/(1-0.5*dR*gunc(k+1)) k = 2 W(i,j,k+1) = W(i,j,k) + (F5o12*rhs(k+1) + F2o3*(rhs(k)+W(i,j,k)*gunc(k)) - F1o12*(rhs(k-1)+W(i,j,k-1)*gunc(k-1)))*dR W(i,j,k+1) = W(i,j,k+1)/(1-F5o12*dR*gunc(k+1)) do k=3,ex(3)-2 W(i,j,k+1) = W(i,j,k) + (F3o8*rhs(k+1) + F19o24*(rhs(k)+W(i,j,k)*gunc(k)) - F5o24*(rhs(k-1)+W(i,j,k-1)*gunc(k-1)) & + F1o24*(rhs(k-2)+W(i,j,k-2)*gunc(k-2)))*dR W(i,j,k+1) = W(i,j,k+1)/(1-F3o8*dR*gunc(k+1)) enddo k = ex(3) W(i,j,k) = dreal(bDCU(i,j,k)) enddo enddo gont = 0 return end function NullEvol_W !-------------------------------------------------------------------- ! this R is indeed x function NullEvol_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,KK,HKK,KKx,HKKx,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W,KK,KKx,HKK,HKKx real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex,dimension(ex(3)) :: CTheta0,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: Cnu,Ck,CUx,DCUx,bDCUx double complex,dimension(ex(3)) :: Cnux,CJx,CJxx,DCJx real*8,dimension(ex(3)) :: betax,Wx,gunc double complex :: Theta_rhs,ZEO real*8 :: dR real*8,parameter :: F5o12=2.d0/1.2d1,F2o3=2.d0/3.d0,F1o12=1.d0/1.2d1 real*8,parameter :: F3o8=3.d0/8.d0,F19o24=1.9d1/2.4d1,F5o24=5.d0/2.4d1,F1o24=1.d0/2.4d1 !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) + & sum(KK)+sum(KKx)+sum(W) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" if(sum(KK).ne.sum(KK))write(*,*)"NullEvol_Theta: find NaN in KK" if(sum(KKx).ne.sum(KKx))write(*,*)"NullEvol_Theta: find NaN in KKx" if(sum(W).ne.sum(W))write(*,*)"NullEvol_Theta: find NaN in W" gont = 1 return endif dR = R(2) - R(1) ZEO = dcmplx(0.d0,0.d0) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CTheta0(1) = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) #if 0 call cderivs_sw_x(ex(3),R,Cnu,Cnux) call rderivs_sw_x(ex(3),R,beta(i,j,:),betax) call rderivs_sw_x(ex(3),R,W,Wx) call cderivs_sw_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_sw_x(ex(3),R,CU(i,j,:),CUx) call cderivs_sw_x(ex(3),R,bDCU(i,j,:),bDCUx) call cderivs_sw_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_sw_x(ex(3),R,CJ(i,j,:),CJxx) call cderivs_sw_x(ex(3),R,DCJ(i,j,:),DCJx) #else call cderivs_x(ex(3),R,Cnu,Cnux) call rderivs_x(ex(3),R,beta(i,j,:),betax) call rderivs_x(ex(3),R,W,Wx) call cderivs_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_x(ex(3),R,CU(i,j,:),CUx) call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) #endif do k = 1,ex(3) rhs(k) = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(i,j,k),KKx(i,j,k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), & DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), & CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),ZEO) gunc(k) = -1/R(k)/(1-R(k)) enddo k = 1 CTheta0(k+1) = CTheta0(k) + (RHS(k+1)+RHS(k)+CTheta0(k)*gunc(k))*dR/2 CTheta0(k+1) = CTheta0(k+1)/(1-0.5*dR*gunc(k+1)) k = 2 CTheta0(k+1) = CTheta0(k) + (F5o12*RHS(k+1) + F2o3*(RHS(k)+CTheta0(k)*gunc(k)) - F1o12*(RHS(k-1)+CTheta0(k-1)*gunc(k-1)))*dR CTheta0(k+1) = CTheta0(k+1)/(1-F5o12*dR*gunc(k+1)) do k=3,ex(3)-2 CTheta0(k+1) = CTheta0(k) + (F3o8*RHS(k+1) + F19o24*(RHS(k)+CTheta0(k)*gunc(k)) - F5o24*(RHS(k-1)+CTheta0(k-1)*gunc(k-1)) & + F1o24*(RHS(k-2)+CTheta0(k-2)*gunc(k-2)))*dR CTheta0(k+1) = CTheta0(k+1)/(1-F3o8*dR*gunc(k+1)) enddo k = ex(3) CTheta0(k) = -KK(i,j,k)*DCU(i,j,k)-(CU(i,j,k)*Cnu(k)+dconjg(CU(i,j,k))*DCJ(i,j,k))/2 & +CJ(i,j,k)*(bDCU(i,j,k)-dconjg(bDCU(i,j,k)))/2 - W(i,j,k)*CJ(i,j,k)/2 RTheta(i,j,:) = dreal(CTheta0) ITheta(i,j,:) = dimag(CTheta0) enddo enddo gont = 0 return end function NullEvol_Theta #else #error "not recognized RKorAM" #endif !===================================================================================================================================== ! basic tool routines subroutine rswap(r1,r2) implicit none !~~~~~~% Input parameters: real*8,intent(inout) :: r1,r2 real*8 :: r r = r1 r1= r2 r2= r return end subroutine rswap !---- subroutine cswap(r1,r2) implicit none !~~~~~~% Input parameters: double complex,intent(inout) :: r1,r2 double complex :: r r = r1 r1= r2 r2= r return end subroutine cswap ! center type finite difference !==================================================================================== !---- subroutine rderivs_x(lx,X,f,fx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(lx) :: f real*8,intent(out),dimension(lx) :: fx real*8 :: dX dX = X(2)-X(1) #ifdef OLD fx(1:lx-1) = (f(2:lx)-f(1:lx-1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #else #if (ghost_width == 2) fx(2:lx-1) = (f(3:lx)-f(1:lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #elif (ghost_width == 3) fx(3:lx-2) = (f(1:lx-4)-8.d0*f(2:lx-3)+8.d0*f(4:lx-1)-f(5:lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX ! fx(1) =-(2.5d1*f(1)-4.8d1*f(2)+3.6d1*f(3)-1.6d1*f(4)+3.d0*f(5))/1.2d1/dX ! fx(2) =-(3.d0*f(1)+1.d1*f(2)-1.8d1*f(3)+6.d0*f(4)-f(5))/1.2d1/dX #elif (ghost_width == 4) fx(4:lx-3) = (-f(1:lx-6)+9.d0*f(2:lx-5)-4.5d1*f(3:lx-4)+4.5d1*f(5:lx-2)-9.d0*f(6:lx-1)+f(7:lx))/6.d1/dX fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #elif (ghost_width == 5) fx(5:lx-4) = (3.d0*f(1:lx-8)-3.2d1*f(2:lx-7)+1.68d2*f(3:lx-6)-6.72d2*f(4:lx-5)+ & 6.72d2*f(6:lx-3)-1.68d2*f(7:lx-2)+3.2d1*f(8:lx-1)-3.d0*f(9:lx))/8.4d2/dX fx(4) = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX fx(lx-3) = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #endif #endif return end subroutine rderivs_x !---- subroutine rderivs_x_point(lx,X,f,fx,k) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx,k real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(lx) :: f real*8,intent(out) :: fx real*8 :: dX dX = X(2)-X(1) #ifdef OLD if(k .eq. lx)then fx = (f(lx)-f(lx-1))/dX else fx = (f(k+1)-f(k))/dX endif #else #if (ghost_width == 2) if(k .gt. 1 .and. k .lt. lx) then fx = (f(k+1)-f(k-1))/2.d0/dX elseif(k.eq.1) then fx = (f(2)-f(1))/dX elseif(k.eq.lx) then fx = (f(lx)-f(lx-1))/dX endif #elif (ghost_width == 3) if(k .gt. 2 .and. k .lt. lx-1) then fx = (f(k-2)-8.d0*f(k-1)+8.d0*f(k+1)-f(k+2))/1.2d1/dX elseif(k.eq.1) then fx = (f(2)-f(1))/dX elseif(k.eq.lx) then fx = (f(lx)-f(lx-1))/dX elseif(k.eq.2) then fx = (f(3)-f(1))/2.d0/dX elseif(k.eq.lx-1) then fx = (f(lx)-f(lx-2))/2.d0/dX endif #elif (ghost_width == 4) if(k .gt. 3 .and. k .lt. lx-2) then fx = (-f(k-3)+9.d0*f(k-2)-4.5d1*f(k-1)+4.5d1*f(k+1)-9.d0*f(k+2)+f(k+3))/6.d1/dX elseif(k.eq.1) then fx = (f(2)-f(1))/dX elseif(k.eq.lx) then fx = (f(lx)-f(lx-1))/dX elseif(k.eq.2) then fx = (f(3)-f(1))/2.d0/dX elseif(k.eq.lx-1) then fx = (f(lx)-f(lx-2))/2.d0/dX elseif(k.eq.3) then fx = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX elseif(k.eq.lx-2) then fx = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX endif #elif (ghost_width == 5) if(k .gt. 4 .and. k .lt. lx-3) then fx = (3.d0*f(k-4)-3.2d1*f(k-3)+1.68d2*f(k-2)-6.72d2*f(k-1)+ & 6.72d2*f(k+1)-1.68d2*f(k+2)+3.2d1*f(k+3)-3.d0*f(k+4))/8.4d2/dX elseif(k.eq.1) then fx = (f(2)-f(1))/dX elseif(k.eq.lx) then fx = (f(lx)-f(lx-1))/dX elseif(k.eq.2) then fx = (f(3)-f(1))/2.d0/dX elseif(k.eq.lx-1) then fx = (f(lx)-f(lx-2))/2.d0/dX elseif(k.eq.3) then fx = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX elseif(k.eq.lx-2) then fx = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX elseif(k.eq.4) then fx = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX elseif(k.eq.lx-3) then fx = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX endif #endif #endif return end subroutine rderivs_x_point !---- subroutine cderivs_x(lx,X,f,fx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X double complex,intent(in),dimension(lx) :: f double complex,intent(out),dimension(lx) :: fx real*8 :: dX dX = X(2)-X(1) #ifdef OLD fx(1:lx-1) = (f(2:lx)-f(1:lx-1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #else #if (ghost_width == 2) fx(2:lx-1) = (f(3:lx)-f(1:lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #elif (ghost_width == 3) fx(3:lx-2) = (f(1:lx-4)-8.d0*f(2:lx-3)+8.d0*f(4:lx-1)-f(5:lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX ! fx(1) =-(2.5d1*f(1)-4.8d1*f(2)+3.6d1*f(3)-1.6d1*f(4)+3.d0*f(5))/1.2d1/dX ! fx(2) =-(3.d0*f(1)+1.d1*f(2)-1.8d1*f(3)+6.d0*f(4)-f(5))/1.2d1/dX #elif (ghost_width == 4) fx(4:lx-3) = (-f(1:lx-6)+9.d0*f(2:lx-5)-4.5d1*f(3:lx-4)+4.5d1*f(5:lx-2)-9.d0*f(6:lx-1)+f(7:lx))/6.d1/dX fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #elif (ghost_width == 5) fx(5:lx-4) = (3.d0*f(1:lx-8)-3.2d1*f(2:lx-7)+1.68d2*f(3:lx-6)-6.72d2*f(4:lx-5)+ & 6.72d2*f(6:lx-3)-1.68d2*f(7:lx-2)+3.2d1*f(8:lx-1)-3.d0*f(9:lx))/8.4d2/dX fx(4) = (-f(1)+9.d0*f(2)-4.5d1*f(3)+4.5d1*f(5)-9.d0*f(6)+f(7))/6.d1/dX fx(lx-3) = (-f(lx-6)+9.d0*f(lx-5)-4.5d1*f(lx-4)+4.5d1*f(lx-2)-9.d0*f(lx-1)+f(lx))/6.d1/dX fx(3) = (f(1)-8.d0*f(2)+8.d0*f(4)-f(5))/1.2d1/dX fx(lx-2) = (f(lx-4)-8.d0*f(lx-3)+8.d0*f(lx-1)-f(lx))/1.2d1/dX fx(2) = (f(3)-f(1))/2.d0/dX fx(lx-1) = (f(lx)-f(lx-2))/2.d0/dX fx(1) = (f(2)-f(1))/dX fx(lx) = (f(lx)-f(lx-1))/dX #endif #endif return end subroutine cderivs_x !---- subroutine cdderivs_x(lx,X,f,fxx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X double complex,intent(in),dimension(lx) :: f double complex,intent(out),dimension(lx) :: fxx real*8 :: dX dX = X(2)-X(1) dX = dX*dX #ifdef OLD fxx(1:lx-2) = (f(3:lx)-2.0*f(2:lx-1)+f(1:lx-2))/dX fxx(lx-1) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX fxx(lx ) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX #else #if (ghost_width == 2) fxx(2:lx-1) = (f(3:lx)-2.d0*f(2:lx-1)+f(1:lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 3) fxx(3:lx-2) = (-f(1:lx-4)+1.6d1*f(2:lx-3)-3.d1*f(3:lx-2)+1.6d1*f(4:lx-1)-f(5:lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 4) fxx(4:lx-3) = (2.d0*f(1:lx-6)-2.7d1*f(2:lx-5)+2.7d2*f(3:lx-4)-4.9d2*f(4:lx-3) & +2.7d2*f(5:lx-2)-2.7d1*f(6:lx-1)+2.d0*f(7:lx))/1.8d2/dX fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 5) fxx(5:lx-4) = (-9.d0*f(1:lx-8)+1.28d2*f(2:lx-7)-1.008d3*f(3:lx-6)+8.064d3*f(4:lx-5)-1.435d4*f(5:lx-4) & +8.064d3*f(6:lx-3)-1.008d3*f(7:lx-2)+1.28d2*f(8:lx-1)-9.d0*f(9:lx))/5.04d3/dX fxx(4) = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) & +2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX fxx(lx-3) = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) & +2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #endif #endif return end subroutine cdderivs_x !---- subroutine rdderivs_x(lx,X,f,fxx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(lx) :: f real*8,intent(out),dimension(lx) :: fxx real*8 :: dX dX = X(2)-X(1) dX = dX*dX #ifdef OLD fxx(1:lx-2) = (f(3:lx)-2.0*f(2:lx-1)+f(1:lx-2))/dX fxx(lx-1) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX fxx(lx ) = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX #else #if (ghost_width == 2) fxx(2:lx-1) = (f(3:lx)-2.d0*f(2:lx-1)+f(1:lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 3) fxx(3:lx-2) = (-f(1:lx-4)+1.6d1*f(2:lx-3)-3.d1*f(3:lx-2)+1.6d1*f(4:lx-1)-f(5:lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 4) fxx(4:lx-3) = (2.d0*f(1:lx-6)-2.7d1*f(2:lx-5)+2.7d2*f(3:lx-4)-4.9d2*f(4:lx-3) & +2.7d2*f(5:lx-2)-2.7d1*f(6:lx-1)+2.d0*f(7:lx))/1.8d2/dX fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #elif (ghost_width == 5) fxx(5:lx-4) = (-9.d0*f(1:lx-8)+1.28d2*f(2:lx-7)-1.008d3*f(3:lx-6)+8.064d3*f(4:lx-5)-1.435d4*f(5:lx-4) & +8.064d3*f(6:lx-3)-1.008d3*f(7:lx-2)+1.28d2*f(8:lx-1)-9.d0*f(9:lx))/5.04d3/dX fxx(4) = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) & +2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX fxx(lx-3) = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) & +2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX fxx(3) = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX fxx(lx-2) = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX fxx(2) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx-1) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX fxx(1) = (f(3)-2.d0*f(2)+f(1))/dX fxx(lx) = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX #endif #endif return end subroutine rdderivs_x !---- subroutine rdderivs_x_point(lx,X,f,fxx,k) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx,k real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(lx) :: f real*8,intent(out) :: fxx real*8 :: dX dX = X(2)-X(1) dX = dX*dX #ifdef OLD if(k.lt.lx-1) then fxx = (f(k+2)-2.0*f(k+1)+f(k))/dX elseif(k.eq.lx-1) then fxx = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX elseif(k.eq.lx) then fxx = (f(lx)-2.0*f(lx-1)+f(lx-2))/dX endif #else #if (ghost_width == 2) if(k.gt.1 .and. k.lt.lx) then fxx = (f(k+1)-2.d0*f(k)+f(k-1))/dX elseif(k.eq.1) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX endif #elif (ghost_width == 3) if(k.gt.2 .and. k.lt.lx-1) then fxx = (-f(k-2)+1.6d1*f(k-1)-3.d1*f(k)+1.6d1*f(k+1)-f(k+2))/1.2d1/dX elseif(k.eq.1) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX elseif(k.eq.2) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx-1) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX endif #elif (ghost_width == 4) if(k.gt.3 .and. k.lt.lx-2)then fxx = (2.d0*f(k-3)-2.7d1*f(k-2)+2.7d2*f(k-1)-4.9d2*f(k) & +2.7d2*f(k+1)-2.7d1*f(k+2)+2.d0*f(k+3))/1.8d2/dX elseif(k.eq.1) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX elseif(k.eq.2) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx-1) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX elseif(k.eq.3) then fxx = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX elseif(k.eq.lx-2) then fxx = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX endif #elif (ghost_width == 5) if(k.gt.4 .and. k.lt.lx-3) then fxx = (-9.d0*f(k-4)+1.28d2*f(k-3)-1.008d3*f(k-2)+8.064d3*f(k-1)-1.435d4*f(k) & +8.064d3*f(k+1)-1.008d3*f(k+2)+1.28d2*f(k+3)-9.d0*f(k+4))/5.04d3/dX elseif(k.eq.1) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX elseif(k.eq.2) then fxx = (f(3)-2.d0*f(2)+f(1))/dX elseif(k.eq.lx-1) then fxx = (f(lx)-2.d0*f(lx-1)+f(lx-2))/dX elseif(k.eq.3) then fxx = (-f(1)+1.6d1*f(2)-3.d1*f(3)+1.6d1*f(4)-f(5))/1.2d1/dX elseif(k.eq.lx-2) then fxx = (-f(lx-4)+1.6d1*f(lx-3)-3.d1*f(lx-2)+1.6d1*f(lx-1)-f(lx))/1.2d1/dX elseif(k.eq.4) then fxx = (2.d0*f(1)-2.7d1*f(2)+2.7d2*f(3)-4.9d2*f(4) & +2.7d2*f(5)-2.7d1*f(6)+2.d0*f(7))/1.8d2/dX elseif(k.eq.lx-3) then fxx = (2.d0*f(lx-6)-2.7d1*f(lx-5)+2.7d2*f(lx-4)-4.9d2*f(lx-3) & +2.7d2*f(lx-2)-2.7d1*f(lx-1)+2.d0*f(lx))/1.8d2/dX endif #endif #endif return end subroutine rdderivs_x_point !---- subroutine rdderivs_xy_point(lx,ly,X,Y,f,fxy,i,j) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx,ly,i,j real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(ly) :: Y real*8,intent(in),dimension(lx,ly) :: f real*8,intent(out) :: fxy real*8 :: dX,dY dX = X(2)-X(1) dY = Y(2)-Y(1) !! we only consider inner points #if (ghost_width == 2) if(i>1 .and. j>1.and.i<lx .and.j<ly)then fxy = (f(i-1,j-1)-f(i+1,j-1)-f(i-1,j+1)+f(i+1,j+1))/(4*dX*dY) else fxy = 0.d0 endif #elif (ghost_width == 3) if(i>2 .and. j>2.and.i<lx-1 .and.j<ly-1)then fxy = ( (f(i-2,j-2)-8*f(i-1,j-2)+8*f(i+1,j-2)-f(i+2,j-2)) & -8 *(f(i-2,j-1)-8*f(i-1,j-1)+8*f(i+1,j-1)-f(i+2,j-1)) & +8 *(f(i-2,j+1)-8*f(i-1,j+1)+8*f(i+1,j+1)-f(i+2,j+1)) & - (f(i-2,j+2)-8*f(i-1,j+2)+8*f(i+1,j+2)-f(i+2,j+2)))/(144*dX*dY) else fxy = 0.d0 endif #elif (ghost_width == 4) if(i>3 .and. j>3.and.i<lx-2 .and.j<ly-2)then fxy = (- (-f(i-3,j-3)+9*f(i-2,j-3)-45*f(i-1,j-3)+45*f(i+1,j-3)-9*f(i+2,j-3)+f(i+3,j-3)) & +9 *(-f(i-3,j-2)+9*f(i-2,j-2)-45*f(i-1,j-2)+45*f(i+1,j-2)-9*f(i+2,j-2)+f(i+3,j-2)) & -45*(-f(i-3,j-1)+9*f(i-2,j-1)-45*f(i-1,j-1)+45*f(i+1,j-1)-9*f(i+2,j-1)+f(i+3,j-1)) & +45*(-f(i-3,j+1)+9*f(i-2,j+1)-45*f(i-1,j+1)+45*f(i+1,j+1)-9*f(i+2,j+1)+f(i+3,j+1)) & -9 *(-f(i-3,j+2)+9*f(i-2,j+2)-45*f(i-1,j+2)+45*f(i+1,j+2)-9*f(i+2,j+2)+f(i+3,j+2)) & + (-f(i-3,j+3)+9*f(i-2,j+3)-45*f(i-1,j+3)+45*f(i+1,j+3)-9*f(i+2,j+3)+f(i+3,j+3)))/(3.6d3*dX*dY) else fxy = 0.d0 endif #elif (ghost_width == 5) if(i>4 .and. j>4.and.i<lx-3 .and.j<ly-3)then fxy = ( 3 *( 3*f(i-4,j-4)-32*f(i-3,j-4)+168*f(i-2,j-4)-672*f(i-1,j-4) & -3*f(i+4,j-4)+32*f(i+3,j-4)-168*f(i+2,j-4)+672*f(i+1,j-4)) & -32 *( 3*f(i-4,j-3)-32*f(i-3,j-3)+168*f(i-2,j-3)-672*f(i-1,j-3) & -3*f(i+4,j-3)+32*f(i+3,j-3)-168*f(i+2,j-3)+672*f(i+1,j-3)) & +168*( 3*f(i-4,j-2)-32*f(i-3,j-2)+168*f(i-2,j-2)-672*f(i-1,j-2) & -3*f(i+4,j-2)+32*f(i+3,j-2)-168*f(i+2,j-2)+672*f(i+1,j-2)) & -672*( 3*f(i-4,j-1)-32*f(i-3,j-1)+168*f(i-2,j-1)-672*f(i-1,j-1) & -3*f(i+4,j-1)+32*f(i+3,j-1)-168*f(i+2,j-1)+672*f(i+1,j-1)) & +672*( 3*f(i-4,j+1)-32*f(i-3,j+1)+168*f(i-2,j+1)-672*f(i-1,j+1) & -3*f(i+4,j+1)+32*f(i+3,j+1)-168*f(i+2,j+1)+672*f(i+1,j+1)) & -168*( 3*f(i-4,j+2)-32*f(i-3,j+2)+168*f(i-2,j+2)-672*f(i-1,j+2) & -3*f(i+4,j+2)+32*f(i+3,j+2)-168*f(i+2,j+2)+672*f(i+1,j+2)) & +32 *( 3*f(i-4,j+3)-32*f(i-3,j+3)+168*f(i-2,j+3)-672*f(i-1,j+3) & -3*f(i+4,j+3)+32*f(i+3,j+3)-168*f(i+2,j+3)+672*f(i+1,j+3)) & -3 *( 3*f(i-4,j+4)-32*f(i-3,j+4)+168*f(i-2,j+4)-672*f(i-1,j+4) & -3*f(i+4,j+4)+32*f(i+3,j+4)-168*f(i+2,j+4)+672*f(i+1,j+4)) )/(7.056d5*dX*dY) else fxy = 0.d0 endif #endif return end subroutine rdderivs_xy_point ! side-winded type finite difference (arXiv:1208.3891) !==================================================================================== subroutine rderivs_sw_x(lx,X,f,fx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X real*8,intent(in),dimension(lx) :: f real*8,intent(out),dimension(lx) :: fx real*8 :: dX integer :: i dX = X(2)-X(1) #if (ghost_width == 2) !#error do i=1,lx-3 fx(i) = (-3.d0*f(i)+4.d0*f(i+1)-f(i+2))/2.d0/dX ! 向前差分,2价精度 enddo do i=lx-2,lx-1 fx(i) = (f(i+1)-f(i-1))/2.d0/dX ! 中心差分,2价精度 enddo i=lx fx(i) = (3.d0*f(i)-4.d0*f(i-1)+f(i-2))/2.d0/dX ! 向后差分,2价精度 #elif (ghost_width == 3) do i=1,lx-5 fx(i) = (-2.5d1/1.2d1*f(i)+4*f(i+1)-3*f(i+2)+4.d0/3.d0*f(i+3)-0.25d0*f(i+4))/dX ! 向前差分,4价精度 enddo do i=lx-4,lx-2 fx(i) = (f(i-2)-8.d0*f(i-1)+8.d0*f(i+1)-f(i+2))/1.2d1/dX ! 中心差分,4价精度 enddo do i=lx-1,lx fx(i) = (2.5d1/1.2d1*f(i)-4*f(i-1)+3*f(i-2)-4.d0/3.d0*f(i-3)+0.25d0*f(i-4))/dX ! 向后差分,4价精度 enddo #elif (ghost_width == 4) !#error do i=1,lx-7 fx(i) = (-147.d0*f(i)+360.d0*f(i+1)-450.d0*f(i+2)+400.d0*f(i+3)-225.d0*f(i+4)+72.d0*f(i+5)-10.d0*f(i+6)) /60.d0/dX ! 向前差分,6价精度 enddo do i=lx-6,lx-3 fx(i) = (f(i-3)-9.d0*f(i-2)+45.d0*f(i-1)-45.d0*f(i+1)+9.0*f(i+2)-f(i+3)) /60.d0/dX ! 中心差分,6价精度 enddo do i=lx-2,lx fx(i) = (147.d0*f(i)-360.d0*f(i-1)+450.d0*f(i-2)-400.d0*f(i-3)+225.d0*f(i-4)-72.d0*f(i-5)+10.d0*f(i-6)) /60.d0/dX ! 向后差分,6价精度 enddo #elif (ghost_width == 5) !#error do i=1,lx-9 fx(i) = (-(761.d0/280.d0)*f(i)+8.d0*f(i+1)-14.d0*f(i+2)+(56.d0/3.d0)*f(i+3)-(35.d0/2.d0)*f(i+4) & +(56.d0/5.d0)*f(i+5)-(14.d0/3.d0)*f(i+6)+(8.d0/7.d0)*f(i+7)-(1.d0/8.d0)*f(i+8)) /dX ! 向前差分,8价精度 enddo do i=lx-8,lx-4 fx(i) = (3.d0*f(i-4)-32.d0*f(i-3)+168.d0*f(i-2)-672.d0*f(i-1)+672.d0*f(i+1)-168.d0*f(i+2)+32.d0*f(i+3)-3.d0*f(i+4)) /840.d0/dX ! 中心差分,8价精度 enddo do i=lx-3,lx fx(i) = ((761.d0/280.d0)*f(i)-8.d0*f(i-1)+14.d0*f(i-2)-(56.d0/3.d0)*f(i-3)+(35.d0/2.d0)*f(i-4) & -(56.d0/5.d0)*f(i-5)+(14.d0/3.d0)*f(i-6)-(8.d0/7.d0)*f(i-7)+(1.d0/8.d0)*f(i-8)) /dX ! 向后差分,8价精度 enddo #endif return end subroutine rderivs_sw_x !---- subroutine cderivs_sw_x(lx,X,f,fx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X double complex,intent(in),dimension(lx) :: f double complex,intent(out),dimension(lx) :: fx real*8 :: dX integer :: i dX = X(2)-X(1) #if (ghost_width == 2) !#error do i=1,lx-3 fx(i) = (-3.d0*f(i)+4.d0*f(i+1)-f(i+2))/2.d0/dX ! 向前差分,2价精度 enddo do i=lx-2,lx-1 ! 宽度为 2 fx(i) = (f(i+1)-f(i-1))/2.d0/dX ! 中心差分,2价精度 enddo i=lx fx(i) = (3.d0*f(i)-4.d0*f(i-1)+f(i-2))/2.d0/dX ! 向后差分,2价精度 #elif (ghost_width == 3) do i=1,lx-5 fx(i) = (-2.5d1/1.2d1*f(i)+4*f(i+1)-3*f(i+2)+4.d0/3.d0*f(i+3)-0.25d0*f(i+4))/dX ! 向前差分,4价精度 enddo do i=lx-4,lx-2 ! 宽度为 3 fx(i) = (f(i-2)-8.d0*f(i-1)+8.d0*f(i+1)-f(i+2))/1.2d1/dX ! 中心差分,4价精度 enddo do i=lx-1,lx fx(i) = (2.5d1/1.2d1*f(i)-4*f(i-1)+3*f(i-2)-4.d0/3.d0*f(i-3)+0.25d0*f(i-4))/dX ! 向后差分,4价精度 enddo #elif (ghost_width == 4) !#error do i=1,lx-7 fx(i) = (-147.d0*f(i)+360.d0*f(i+1)-450.d0*f(i+2)+400.d0*f(i+3)-225.d0*f(i+4)+72.d0*f(i+5)-10.d0*f(i+6)) /60.d0/dX ! 向前差分,6价精度 enddo do i=lx-6,lx-3 ! 宽度为 4 fx(i) = (f(i-3)-9.d0*f(i-2)+45.d0*f(i-1)-45.d0*f(i+1)+9.0*f(i+2)-f(i+3)) /60.d0/dX ! 中心差分,6价精度 enddo do i=lx-2,lx fx(i) = (147.d0*f(i)-360.d0*f(i-1)+450.d0*f(i-2)-400.d0*f(i-3)+225.d0*f(i-4)-72.d0*f(i-5)+10.d0*f(i-6)) /60.d0/dX ! 向后差分,6价精度 enddo #elif (ghost_width == 5) !#error do i=1,lx-9 fx(i) = (-(761.d0/280.d0)*f(i)+8.d0*f(i+1)-14.d0*f(i+2)+(56.d0/3.d0)*f(i+3)-(35.d0/2.d0)*f(i+4) & +(56.d0/5.d0)*f(i+5)-(14.d0/3.d0)*f(i+6)+(8.d0/7.d0)*f(i+7)-(1.d0/8.d0)*f(i+8)) /dX ! 向前差分,8价精度 enddo do i=lx-8,lx-4 ! 宽度为 5 fx(i) = (3.d0*f(i-4)-32.d0*f(i-3)+168.d0*f(i-2)-672.d0*f(i-1)+672.d0*f(i+1)-168.d0*f(i+2)+32.d0*f(i+3)-3.d0*f(i+4)) /840.d0/dX ! 中心差分,8价精度 enddo do i=lx-3,lx fx(i) = ((761.d0/280.d0)*f(i)-8.d0*f(i-1)+14.d0*f(i-2)-(56.d0/3.d0)*f(i-3)+(35.d0/2.d0)*f(i-4) & -(56.d0/5.d0)*f(i-5)+(14.d0/3.d0)*f(i-6)-(8.d0/7.d0)*f(i-7)+(1.d0/8.d0)*f(i-8)) /dX ! 向后差分,8价精度 enddo #endif return end subroutine cderivs_sw_x !---- subroutine cdderivs_sw_x(lx,X,f,fxx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: X double complex,intent(in),dimension(lx) :: f double complex,intent(out),dimension(lx) :: fxx real*8 :: dX integer :: i dX = X(2)-X(1) dX = dX*dX #if (ghost_width == 2) !#error do i=1,lx-3 fxx(i) = (f(i)-2.d0*f(i+1)+f(i+2))/dX ! 向前差分,2价精度 enddo do i=lx-2,lx-1 fxx(i) = (f(i-1)-2.d0*f(i)+f(i+1))/dX ! 中心差分,2价精度 enddo i=lx fxx(i) = (f(i)-2.d0*f(i-1)+f(i-2))/dX ! 向后差分,2价精度 #elif (ghost_width == 3) do i=1,lx-5 fxx(i) = ((15.d0/4)*f(i)-(77.d0/6)*f(i+1)+(107.d0/6)*f(i+2)-13.d0*f(i+3)+(61.d0/12.d0)*f(i+4)-(5.d0/6.d0)*f(i+5))/dX ! 向前差分,4价精度 enddo do i=lx-4,lx-2 ! 宽度为 3 fxx(i) = (-f(i-2)+16.d0*f(i-1)-30.d0*f(i)+16.d0*f(i+1)-f(i+2))/12.d0/dX ! 中心差分,4价精度 enddo do i=lx-1,lx fxx(i) = ((15.d0/4)*f(i)-(77.d0/6)*f(i-1)+(107.d0/6)*f(i-2)-13.d0*f(i-3)+(61.d0/12.d0)*f(i-4)-(5.d0/6.d0)*f(i-5))/dX ! 向后差分,4价精度 enddo #elif (ghost_width == 4) do i=1,lx-7 fxx(i) = (812.d0*f(i)-3132.d0*f(i+1)+5265.d0*f(i+2)-5080.d0*f(i+3)+2970.d0*f(i+4)-972.d0*f(i+5)+137.d0*f(i+6)) /180.d0/dX ! 向前差分,6价精度 enddo do i=lx-6,lx-3 ! 宽度为 4 fxx(i) = (2.d0*f(i-3)-27.d0*f(i-2)+270.d0*f(i-1)-490.d0*f(i)+270.d0*f(i+1)-27.d0*f(i+2)+2.d0*f(i+3)) /180.d0/dX ! 中心差分,6价精度 enddo do i=lx-2,lx fxx(i) = (812.d0*f(i)-3132.d0*f(i-1)+5265.d0*f(i-2)-5080.d0*f(i-3)+2970.d0*f(i-4)-972.d0*f(i-5)+137.d0*f(i-6)) /180.d0/dX ! 向后差分,6价精度 enddo #elif (ghost_width == 5) do i=1,lx-9 fxx(i) = ( (29531.d0/5040.d0)*f(i) - (962.d0/35.d0)*f(i+1) + (621.d0/10.d0)*f(i+2) & - (4006.d0/45.d0)*f(i+3) + (691.d0/8.d0)*f(i+4) - (282.d0/5.d0)*f(i+5) & + (2143.d0/90.d0)*f(i+6) - (206.d0/35.d0)*f(i+7) + (363.d0/560.d0)*f(i+8) ) /dX ! 向前差分,8阶精度 enddo do i=lx-8,lx-4 ! 宽度为 5 fxx(i) = ( - 9.d0*f(i-4) + 128.d0*f(i-3) - 1008.d0*f(i-2) + 8064.d0*f(i-1) & - 14350.d0*f(i) + 8064.d0*f(i+1) - 1008.d0*f(i+2) + 128.d0*f(i+3) - 9.d0*f(i+4) ) / 5040.d0/dX ! 中心差分,8阶精度 enddo do i=lx-3,lx fxx(i) = ( (29531.d0/5040.d0)*f(i) - (962.d0/35.d0)*f(i-1) + (621.d0/10.d0)*f(i-2) & - (4006.d0/45.d0)*f(i-3) + (691.d0/8.d0)*f(i-4) - (282.d0/5.d0)*f(i-5) & + (2143.d0/90.d0)*f(i-6) - (206.d0/35.d0)*f(i-7) + (363.d0/560.d0)*f(i-8) ) /dX ! 向后差分,8阶精度 enddo #endif return end subroutine cdderivs_sw_x !-------------------- subroutine rget_half_x(lx,f,hf) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx real*8,intent(in),dimension(lx) :: f real*8,intent(out),dimension(lx) :: hf #ifdef OLD hf(1:lx-1) = (f(2:lx)+f(1:lx-1))/2.d0 hf(lx) = hf(lx-1) #else hf(lx) = 0.d0 #if (ghost_width == 2) hf(1:lx-1) = (f(1:lx-1)+f(2:lx))/2.d0 #elif (ghost_width == 3) hf(2:lx-2) = (-f(1:lx-3)+9.d0*f(2:lx-2)+9.d0*f(3:lx-1)-f(4:lx))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #elif (ghost_width == 4) hf(3:lx-3) = 3.d0/2.56d2*(f(1:lx-5)+f(6:lx))-2.5d1/2.56d2*(f(2:lx-4)+f(5:lx-1))+7.5d1/1.28d2*(f(3:lx-3)+f(4:lx-2)) hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #elif (ghost_width == 5) hf(4:lx-4) = (-5.d0*(f(1:lx-7)+f(8:lx))+4.9d1*(f(2:lx-6)+f(7:lx-1))-2.45d2*(f(3:lx-5)+f(6:lx-2))+1.225d3*(f(4:lx-4)+f(5:lx-3))) & /2.048d3 hf(3) = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4)) hf(lx-3) = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3)) hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #endif #endif return end subroutine rget_half_x !-------------------- subroutine rget_half_x_point(lx,f,hf,k) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx,k real*8,intent(in),dimension(lx) :: f real*8,intent(out) :: hf #ifdef OLD if(k.eq.lx .or. k.eq.lx-1)then hf = (f(lx)+f(lx-1))/2.d0 else hf = (f(k+1)+f(k))/2.d0 endif #else if(k .eq. lx) hf = 0.d0 #if (ghost_width == 2) if(k .lt. lx) hf = (f(k+1)+f(k))/2.d0 #elif (ghost_width == 3) if(k .eq. 1) hf = (f(2)+f(1))/2.d0 if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0 if(k .gt. 1 .and. k .lt. lx-1) hf = (-f(k-1)+9.d0*f(k)+9.d0*f(k+1)-f(k+2))/1.6d1 #elif (ghost_width == 4) if(k .eq. 1) hf = (f(2)+f(1))/2.d0 if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0 if(k .eq. 2) hf = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 if(k .eq. lx-2) hf =(-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 if(k .gt. 2 .and. k .lt. lx-2) hf = 3.d0/2.56d2*(f(k-2)+f(k+3))-2.5d1/2.56d2*(f(k-1)+f(k+2))+7.5d1/1.28d2*(f(k)+f(k+1)) #elif (ghost_width == 5) if(k .eq. 1) hf = (f(2)+f(1))/2.d0 if(k .eq. lx-1) hf = (f(lx)+f(lx-1))/2.d0 if(k .eq. 2) hf = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 if(k .eq. lx-2) hf =(-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 if(k .eq. 3) hf = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4)) if(k .eq. lx-3) hf = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3)) if(k .gt. 3 .and. k .lt. lx-3) then hf = (-5.d0*(f(k-3)+f(k+4))+4.9d1*(f(k-2)+f(k+3))-2.45d2*(f(k-1)+f(k+2))+1.225d3*(f(k)+f(k+1)))/2.048d3 endif #endif #endif return end subroutine rget_half_x_point !-------------------- subroutine cget_half_x(lx,f,hf) implicit none !~~~~~~% Input parameters: integer,intent(in) :: lx double complex,intent(in),dimension(lx) :: f double complex,intent(out),dimension(lx) :: hf #ifdef OLD hf(1:lx-1) = (f(2:lx)+f(1:lx-1))/2.d0 hf(lx) = hf(lx-1) #else hf(lx) = 0.d0 #if (ghost_width == 2) hf(1:lx-1) = (f(1:lx-1)+f(2:lx))/2.d0 #elif (ghost_width == 3) hf(2:lx-2) = (-f(1:lx-3)+9.d0*f(2:lx-2)+9.d0*f(3:lx-1)-f(4:lx))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #elif (ghost_width == 4) hf(3:lx-3) = 3.d0/2.56d2*(f(1:lx-5)+f(6:lx))-2.5d1/2.56d2*(f(2:lx-4)+f(5:lx-1))+7.5d1/1.28d2*(f(3:lx-3)+f(4:lx-2)) hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #elif (ghost_width == 5) hf(4:lx-4) = (-5.d0*(f(1:lx-7)+f(8:lx))+4.9d1*(f(2:lx-6)+f(7:lx-1))-2.45d2*(f(3:lx-5)+f(6:lx-2))+1.225d3*(f(4:lx-4)+f(5:lx-3))) & /2.048d3 hf(3) = 3.d0/2.56d2*(f(1)+f(6))-2.5d1/2.56d2*(f(2)+f(5))+7.5d1/1.28d2*(f(3)+f(4)) hf(lx-3) = 3.d0/2.56d2*(f(lx)+f(lx-5))-2.5d1/2.56d2*(f(lx-1)+f(lx-4))+7.5d1/1.28d2*(f(lx-2)+f(lx-3)) hf(2) = (-f(1)+9.d0*f(2)+9.d0*f(3)-f(4))/1.6d1 hf(lx-2) = (-f(lx-1)+9.d0*f(lx-2)+9.d0*f(lx-3)-f(lx-4))/1.6d1 hf(1) = (f(2)+f(1))/2.d0 hf(lx-1) = (f(lx)+f(lx-1))/2.d0 #endif #endif return end subroutine cget_half_x !--------- subroutine derivs_eth(ex,X,Y,f,eth_f,spin,e, & quR1,quR2,quI1,quI2,gR,gI) implicit none !~~~~~~% Input parameters: ! spin: spin weight of f; e: eth (1) or eth bar (-1) integer,intent(in) :: spin,e integer,intent(in),dimension(2) :: ex real*8,intent(in),dimension(ex(1)) :: X real*8,intent(in),dimension(ex(2)) :: Y double complex,intent(in),dimension(ex(1),ex(2)) :: f double complex,intent(out),dimension(ex(1),ex(2)) :: eth_f real*8,intent(in),dimension(ex(1),ex(2)) :: quR1,quR2,quI1,quI2,gR,gI double complex,dimension(ex(1),ex(2)) :: fx,fy double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama integer :: i,j !sanity check if(e.ne.1 .and.e.ne.-1)then write(*,*) "derivs_eth: bad input e = ",e return endif qu1 = dcmplx(quR1,quI1) qu2 = dcmplx(quR2,quI2) gama = dcmplx(gR,gI) do j=1,ex(2) call cderivs_x(ex(1),X,f(:,j),fx(:,j)) enddo do i=1,ex(1) call cderivs_x(ex(2),Y,f(i,:),fy(i,:)) enddo if(e.eq.1)then eth_f = qu1*fx+qu2*fy+spin*gama*f else eth_f = dconjg(qu1)*fx+dconjg(qu2)*fy-spin*dconjg(gama)*f endif return end subroutine derivs_eth !--------- ! dqu: q^B @_B q^A ! bdqu: bar{q}^B @_B q^A ! dg: q^B @_B Gamma ! bdg: q^B @_B bar{Gamma} subroutine dderivs_eth(ex,X,Y,f,eth_f,spin,e1,e2, & quR1,quR2,quI1,quI2,gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI) implicit none !~~~~~~% Input parameters: ! spin: spin weight of f; e: eth (1) or eth bar (-1) integer,intent(in) :: spin,e1,e2 integer,intent(in),dimension(2) :: ex real*8,intent(in),dimension(ex(1)) :: X real*8,intent(in),dimension(ex(2)) :: Y double complex,intent(in),dimension(ex(1),ex(2)) :: f double complex,intent(out),dimension(ex(1),ex(2)) :: eth_f real*8,intent(in),dimension(ex(1),ex(2)) :: quR1,quR2,quI1,quI2,gR,gI real*8,intent(in),dimension(ex(1),ex(2)) :: dquR1,dquR2,dquI1,dquI2 real*8,intent(in),dimension(ex(1),ex(2)) :: bdquR1,bdquR2,bdquI1,bdquI2 real*8,intent(in),dimension(ex(1),ex(2)) :: dgR,dgI,bdgR,bdgI double complex,dimension(ex(1),ex(2)) :: fx,fy,fxx,fxy,fyy double complex,dimension(ex(1),ex(2)) :: qu1,qu2,dqu1,dqu2,bdqu1,bdqu2,gama,dgama,bdgama integer :: i,j !sanity check if((e1.ne.1 .and. e1.ne.-1).or.(e2.ne.1 .and. e2.ne.-1))then write(*,*) "dderivs_eth: bad input e1 = ",e1,"e2 = ",e2 return endif qu1 = dcmplx(quR1,quI1) qu2 = dcmplx(quR2,quI2) gama = dcmplx(gR,gI) dqu1 = dcmplx(dquR1,dquI1) dqu2 = dcmplx(dquR2,dquI2) bdqu1 = dcmplx(bdquR1,bdquI1) bdqu2 = dcmplx(bdquR2,bdquI2) dgama = dcmplx(dgR,dgI) bdgama = dcmplx(bdgR,bdgI) do j=1,ex(2) call cderivs_x(ex(1),X,f(:,j),fx(:,j)) call cdderivs_x(ex(1),X,f(:,j),fxx(:,j)) enddo do i=1,ex(1) call cderivs_x(ex(2),Y,f(i,:),fy(i,:)) call cdderivs_x(ex(2),Y,f(i,:),fyy(i,:)) enddo do j=1,ex(2) call cderivs_x(ex(1),X,fy(:,j),fxy(:,j)) enddo if(e1.eq.1.and.e2.eq.1)then ! eth_eth eth_f = qu1*qu1*fxx+2.d0*qu1*qu2*fxy+qu2*qu2*fyy & +(dqu1+(2*spin+1)*gama*qu1)*fx+(dqu2+(2*spin+1)*gama*qu2)*fy & +spin*(dgama+(spin+1)*gama*gama)*f elseif(e1.eq.1.and.e2.eq.-1)then ! eth_ethb eth_f = qu1*dconjg(qu1)*fxx+qu1*dconjg(qu2)*fxy+qu2*dconjg(qu1)*fxy+qu2*dconjg(qu2)*fyy & +(dconjg(bdqu1)-spin*dconjg(gama)*qu1+(spin-1)*gama*dconjg(qu1))*fx & +(dconjg(bdqu2)-spin*dconjg(gama)*qu2+(spin-1)*gama*dconjg(qu2))*fy & -spin*(bdgama+(spin-1)*gama*dconjg(gama))*f elseif(e1.eq.-1.and.e2.eq.1)then ! ethb_eth eth_f = qu1*dconjg(qu1)*fxx+qu1*dconjg(qu2)*fxy+qu2*dconjg(qu1)*fxy+qu2*dconjg(qu2)*fyy & +(bdqu1+spin*gama*dconjg(qu1)-(spin+1)*dconjg(gama)*qu1)*fx & +(bdqu2+spin*gama*dconjg(qu2)-(spin+1)*dconjg(gama)*qu2)*fy & +spin*(dconjg(bdgama)-(spin+1)*gama*dconjg(gama))*f else ! ethb_ethb eth_f = dconjg(qu1*qu1)*fxx+2.d0*dconjg(qu1*qu2)*fxy+dconjg(qu2*qu2)*fyy & +(dconjg(dqu1)-(2*spin-1)*dconjg(gama*qu1))*fx+(dconjg(dqu2)-(2*spin-1)*dconjg(gama*qu2))*fy & -spin*(dconjg(dgama)-(spin-1)*dconjg(gama*gama))*f endif return end subroutine dderivs_eth !------------------------------------------------------------------------------------------ ! CQG 24S327, Eq.(19)--(22) !------------------------------------------------------------------------------------------ subroutine setup_dyad(ex,crho,sigma,R, & quR1,quR2,quI1,quI2, & qlR1,qlR2,qlI1,qlI2, & gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI, & gx,gy,gz,sst,Rmin) implicit none !~~~~~~% Input parameters: integer,intent(in ):: ex(1:3),sst real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8,intent(in) :: Rmin real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2 real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2 real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: gR,gI real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2 real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2 real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: gx,gy,gz real*8 :: thetac,thetas,sr,ss,cr,cs,srss,crcs,tcts,tcts2 real*8 :: sr2,ss2,cr2,cs2,tc2,ts2 integer :: i,j,k real*8 :: ggr,tgrho,tgsigma real*8,dimension(ex(1),ex(2),ex(3)) :: tp do j=1,ex(2) do i=1,ex(1) sr = dsin(crho(i)) ss = dsin(sigma(j)) cr = dcos(crho(i)) cs = dcos(sigma(j)) srss = sr*ss crcs = cr*cs sr2 = sr*sr ss2 = ss*ss cr2 = cr*cr cs2 = cs*cs thetac = dsqrt((1.d0-srss)/2.d0) thetas = dsqrt((1.d0+srss)/2.d0) tc2 = thetac*thetac ts2 = thetas*thetas tcts = thetac*thetas tcts2 = tcts*tcts qlR1(i,j,:) = thetac*cs/4.d0/tcts2 qlI1(i,j,:) = thetas*cs/4.d0/tcts2 qlR2(i,j,:) = thetac*cr/4.d0/tcts2 qlI2(i,j,:) =-thetas*cr/4.d0/tcts2 quR1(i,j,:) = 2.d0*tcts*thetas/cs quI1(i,j,:) = 2.d0*tcts*thetac/cs quR2(i,j,:) = 2.d0*tcts*thetas/cr quI2(i,j,:) =-2.d0*tcts*thetac/cr gR(i,j,:) = (crcs*crcs*(sr+ss)+(cr*cr-cs*cs)*(ss-sr))/4.d0/thetac/crcs gI(i,j,:) = (crcs*crcs*(sr-ss)+(cs*cs-cr*cr)*(ss+sr))/4.d0/thetas/crcs dquR1(i,j,:) = 3.d0*cr2*ss*(cs2+cr2*ss2)/2.d0/cr/cs2 dquI1(i,j,:) =-thetac*thetas*sr*(cs2+3.d0*cr2*ss2)/cr/cs2 dquR2(i,j,:) = 3.d0*cs2*sr*(cr2+cs2*sr2)/2.d0/cs/cr2 dquI2(i,j,:) = thetac*thetas*ss*(cr2+3.d0*cs2*sr2)/cs/cr2 bdquR1(i,j,:) = sr*cs2*(cs2+cr2*ss2)/2.d0/cr/cs2 bdquI1(i,j,:) = thetac*thetas*ss*(8.d0*tc2*ts2-3.d0*cs2*sr2-cr2)/cr/cs2 bdquR2(i,j,:) = ss*cr2*(cr2+cs2*sr2)/2.d0/cs/cr2 bdquI2(i,j,:) =-thetac*thetas*sr*(8.d0*tc2*ts2-3.d0*cr2*ss2-cs2)/cs/cr2 ggr = tc2*ts2/cr2/cs2 dgR(i,j,:) =-3.d0*ggr*sr*ss*(cs2+cr2) dgI(i,j,:) = 6.d0*ggr*thetac*thetas*(cs2-cr2) bdgR(i,j,:) = ggr*(2.d0*cr2*cs2+cr2+cs2) bdgI(i,j,:) = 0.d0 do k=1,ex(3) ggr = R(k)*Rmin/(1.d0-R(k)) tgrho = dtan(crho(i)) tgsigma = dtan(sigma(j)) select case (sst) case (0) gz(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gx(i,j,k) = gz(i,j,k)*tgrho gy(i,j,k) = gz(i,j,k)*tgsigma case (1) gz(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gx(i,j,k) = gz(i,j,k)*tgrho gy(i,j,k) = gz(i,j,k)*tgsigma case (2) gx(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gy(i,j,k) = gx(i,j,k)*tgrho gz(i,j,k) = gx(i,j,k)*tgsigma case (3) gx(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gy(i,j,k) = gx(i,j,k)*tgrho gz(i,j,k) = gx(i,j,k)*tgsigma case (4) gy(i,j,k) = ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gx(i,j,k) = gy(i,j,k)*tgrho gz(i,j,k) = gy(i,j,k)*tgsigma case (5) gy(i,j,k) = -ggr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) gx(i,j,k) = gy(i,j,k)*tgrho gz(i,j,k) = gy(i,j,k)*tgsigma case default write(*,*) "setup_dyad: not recognized sst = ",sst return end select enddo enddo enddo ! note the difference between right hand coordinate and left hand coordinate if(sst==1 .or. sst==3 .or. sst==4)then quI1 = -quI1 quI2 = -quI2 qlI1 = -qlI1 qlI2 = -qlI2 gI = -gI dquI1 = -dquI1 dquI2 = -dquI2 bdquI1 = -bdquI1 bdquI2 = -bdquI2 dgI = -dgI bdgI = -bdgI endif return end subroutine setup_dyad !--------------------------------------------------------------------------------- ! interface to c !--------------------------------------------------------------------------------- subroutine eth_derivs(ex,X,Y,Rfi,Ifi,Reth_f,Ieth_f,spin,e, & quR1,quR2,quI1,quI2,gR,gI) implicit none !~~~~~~% Input parameters: ! spin: spin weight of f; e: eth (1) or eth bar (-1) integer,intent(in) :: spin,e integer,intent(in),dimension(3) :: ex real*8,intent(in),dimension(ex(1)) :: X real*8,intent(in),dimension(ex(2)) :: Y real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: Rfi,Ifi real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: Reth_f,Ieth_f real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2,gR,gI double complex,dimension(ex(1),ex(2)) :: f,eth_f,fx,fy double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama integer :: k do k=1,ex(3) f = dcmplx(Rfi(:,:,k),Ifi(:,:,k)) call derivs_eth(ex(1:2),X,Y,f,eth_f,spin,e, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) Reth_f(:,:,k) = dreal(eth_f) Ieth_f(:,:,k) = dimag(eth_f) enddo return end subroutine eth_derivs !--------------------------------------------------------------------------------- subroutine eth_dderivs(ex,X,Y,Rfi,Ifi,Reth_f,Ieth_f,spin,e1,e2, & quR1,quR2,quI1,quI2,gR,gI, & dquR1,dquR2,dquI1,dquI2, & bdquR1,bdquR2,bdquI1,bdquI2, & dgR,dgI,bdgR,bdgI) implicit none !~~~~~~% Input parameters: ! spin: spin weight of f; e: eth (1) or eth bar (-1) integer,intent(in) :: spin,e1,e2 integer,intent(in),dimension(3) :: ex real*8,intent(in),dimension(ex(1)) :: X real*8,intent(in),dimension(ex(2)) :: Y real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: Rfi,Ifi real*8,intent(out),dimension(ex(1),ex(2),ex(3)) :: Reth_f,Ieth_f real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2,gR,gI real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dquR1,dquR2,dquI1,dquI2 real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: bdquR1,bdquR2,bdquI1,bdquI2 real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: dgR,dgI,bdgR,bdgI double complex,dimension(ex(1),ex(2)) :: f,eth_f,fx,fy double complex,dimension(ex(1),ex(2)) :: qu1,qu2,gama integer :: k do k=1,ex(3) f = dcmplx(Rfi(:,:,k),Ifi(:,:,k)) call dderivs_eth(ex(1:2),X,Y,f,eth_f,spin,e1,e2, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k), & dquR1(:,:,k),dquR2(:,:,k),dquI1(:,:,k),dquI2(:,:,k), & bdquR1(:,:,k),bdquR2(:,:,k),bdquI1(:,:,k),bdquI2(:,:,k), & dgR(:,:,k),dgI(:,:,k),bdgR(:,:,k),bdgI(:,:,k)) Reth_f(:,:,k) = dreal(eth_f) Ieth_f(:,:,k) = dimag(eth_f) enddo return end subroutine eth_dderivs !--------------------------------------------------------------------------------- ! fill symmetric boundary buffer points !--------------------------------------------------------------------------------- subroutine fill_symmetric_boundarybuffer(ex,crho,sigma,R,drho,dsigma, & quR1,quR2,quI1,quI2,qlR1,qlR2,qlI1,qlI2, & Rv,Iv,Symmetry,sst,spin) implicit none !~~~~~~% Input parameters: integer,intent(in),dimension(3) :: ex integer,intent(in) :: Symmetry,sst,spin real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8,intent(in) :: drho,dsigma real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: quR1,quR2,quI1,quI2 real*8,intent(in),dimension(ex(1),ex(2),ex(3)) :: qlR1,qlR2,qlI1,qlI2 real*8,intent(inout),dimension(ex(1),ex(2),ex(3)) :: Rv,Iv integer :: i,j,k,t select case (Symmetry) case (0) return case (1) if((sst==2.or.sst==4).and.dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then do k=1,ex(3) do j=1,ghost_width do i=1,ex(1) #ifdef Vertex #ifdef Cell #error Both Cell and Vertex are defined #endif t = 2*ghost_width+2-j #endif #ifdef Cell #ifdef Vertex #error Both Cell and Vertex are defined #endif t = 2*ghost_width+1-j #endif call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), & quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), & Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin) enddo enddo enddo endif if((sst==3.or.sst==5).and.dabs(sigma(ex(2))-ghost_width*dsigma) < dsigma/2.d0)then do k=1,ex(3) do j=ex(2)-ghost_width+1,ex(2) do i=1,ex(1) t = ex(2)-j+1 #ifdef Vertex #ifdef Cell #error Both Cell and Vertex are defined #endif t = ex(2)-2*ghost_width-1+t #endif #ifdef Cell #ifdef Vertex #error Both Cell and Vertex are defined #endif t = ex(2)-2*ghost_width+t #endif call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), & quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), & Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin) enddo enddo enddo endif case (2) if(dabs(crho(1)+ghost_width*drho) < drho/2.d0)then if(dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then do k=1,ex(3) do j=1,ghost_width do i=ghost_width+1,ex(1) #ifdef Vertex #ifdef Cell #error Both Cell and Vertex are defined #endif t = 2*ghost_width+2-j #endif #ifdef Cell #ifdef Vertex #error Both Cell and Vertex are defined #endif t = 2*ghost_width+1-j #endif call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), & quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), & Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin) enddo enddo enddo endif do k=1,ex(3) do j=1,ex(2) do i=1,ghost_width #ifdef Vertex #ifdef Cell #error Both Cell and Vertex are defined #endif t = 2*ghost_width+2-i #endif #ifdef Cell #ifdef Vertex #error Both Cell and Vertex are defined #endif t = 2*ghost_width+1-i #endif call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), & quR1(t,j,k),quR2(t,j,k),quI1(t,j,k),quI2(t,j,k),qlR1(t,j,k),qlR2(t,j,k),qlI1(t,j,k),qlI2(t,j,k), & Rv(i,j,k),Iv(i,j,k),Rv(t,j,k),Iv(t,j,k),1,spin) enddo enddo enddo else if(dabs(sigma(1)+ghost_width*dsigma) < dsigma/2.d0)then do k=1,ex(3) do j=1,ghost_width do i=1,ex(1) #ifdef Vertex #ifdef Cell #error Both Cell and Vertex are defined #endif t = 2*ghost_width+2-j #endif #ifdef Cell #ifdef Vertex #error Both Cell and Vertex are defined #endif t = 2*ghost_width+1-j #endif call symmetrymap(quR1(i,j,k),quR2(i,j,k),quI1(i,j,k),quI2(i,j,k),qlR1(i,j,k),qlR2(i,j,k),qlI1(i,j,k),qlI2(i,j,k), & quR1(i,t,k),quR2(i,t,k),quI1(i,t,k),quI2(i,t,k),qlR1(i,t,k),qlR2(i,t,k),qlI1(i,t,k),qlI2(i,t,k), & Rv(i,j,k),Iv(i,j,k),Rv(i,t,k),Iv(i,t,k),2,spin) enddo enddo enddo endif endif end select return end subroutine fill_symmetric_boundarybuffer !------------------------------------------------------------------------------------------- ! note Theta has spin weight 2, but the its determinant does not satisfy our ! assumption, so in principle this routine can not be applied to Theta ! but since we never calculate eth derivative on Theta, it does not matter !------------------------------------------------------------------------------------------- subroutine symmetrymap(oquR1,oquR2,oquI1,oquI2,oqlR1,oqlR2,oqlI1,oqlI2, & iquR1,iquR2,iquI1,iquI2,iqlR1,iqlR2,iqlI1,iqlI2, & oRv,oIv,iRv,iIv,ind,spin) implicit none !~~~~~~% Input parameters: integer,intent(in) :: ind,spin real*8,intent(in) :: iquR1,iquR2,iquI1,iquI2,iqlR1,iqlR2,iqlI1,iqlI2 real*8,intent(in) :: iRv,iIv real*8,intent(in) :: oquR1,oquR2,oquI1,oquI2,oqlR1,oqlR2,oqlI1,oqlI2 real*8,intent(out) :: oRv,oIv double complex :: iqu1,iqu2,iql1,iql2,oqu1,oqu2,oql1,oql2 real*8 :: h11,h12,h22,KK double complex :: r iqu1 = dcmplx(iquR1,iquI1) iqu2 = dcmplx(iquR2,iquI2) iql1 = dcmplx(iqlR1,iqlI1) iql2 = dcmplx(iqlR2,iqlI2) oqu1 = dcmplx(oquR1,oquI1) oqu2 = dcmplx(oquR2,oquI2) oql1 = dcmplx(oqlR1,oqlI1) oql2 = dcmplx(oqlR2,oqlI2) ! spin weight select case (spin) case (-2) r = dcmplx(iRv,iIv) #if 1 KK=dsqrt(iRv*iRv+iIv*iIv+1.d0) h11 = dreal(r*iql1*iql1+KK*dconjg(iql1)*iql1) h12 =-dreal(r*iql1*iql2+KK*dconjg(iql1)*iql2) h22 = dreal(r*iql2*iql2+KK*dconjg(iql2)*iql2) #else h11 = dreal(r*iql1*iql1) h12 =-dreal(r*iql1*iql2) h22 = dreal(r*iql2*iql2) #endif r = h11*dconjg(oqu1)*dconjg(oqu1)+2.d0*h12*dconjg(oqu1)*dconjg(oqu2)+h22*dconjg(oqu2)*dconjg(oqu2) oRv = dreal(r) oIv = dimag(r) case (-1) r = dcmplx(iRv,iIv) h11 = dreal(r*iql1) h22 = dreal(r*iql2) if(ind == 1) h11 = -h11 if(ind == 2) h22 = -h22 r = h11*dconjg(oqu1)+h22*dconjg(oqu2) oRv = dreal(r) oIv = dimag(r) case (0) oRv = iRv oIv = iIv case(1) r = dcmplx(iRv,iIv) h11 = dreal(r*dconjg(iql1)) h22 = dreal(r*dconjg(iql2)) if(ind == 1) h11 = -h11 if(ind == 2) h22 = -h22 r = h11*oqu1+h22*oqu2 oRv = dreal(r) oIv = dimag(r) case (2) r = dcmplx(iRv,iIv) #if 1 KK=dsqrt(iRv*iRv+iIv*iIv+1.d0) h11 = dreal(r*dconjg(iql1)*dconjg(iql1)+KK*dconjg(iql1)*iql1) h12 =-dreal(r*dconjg(iql1)*dconjg(iql2)+KK*dconjg(iql1)*iql2) h22 = dreal(r*dconjg(iql2)*dconjg(iql2)+KK*dconjg(iql2)*iql2) #else h11 = dreal(r*dconjg(iql1)*dconjg(iql1)) h12 =-dreal(r*dconjg(iql1)*dconjg(iql2)) h22 = dreal(r*dconjg(iql2)*dconjg(iql2)) #endif r = h11*oqu1*oqu1+2.d0*h12*oqu1*oqu2+h22*oqu2*oqu2 oRv = dreal(r) oIv = dimag(r) case default write(*,*)"symmetrymap does not yet support spin weight ",spin end select return end subroutine symmetrymap !------------------------- subroutine calculate_K(ex,X,Y,Z,RJ,IJ,KK,HKK,KKx,HKKx) implicit none !~~~~~~% Input parameters: integer,intent(in) :: ex(3) real*8,intent(in) :: X(ex(1)),Y(ex(2)),Z(ex(3)) real*8,dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ real*8,dimension(ex(1),ex(2),ex(3)),intent(out) :: KK,HKK,KKx,HKKx integer :: i,j KK = dsqrt(1.d0+RJ*RJ+IJ*IJ) do j=1,ex(2) do i=1,ex(1) call rget_half_x(ex(3),KK(i,j,:),HKK(i,j,:)) call rderivs_x(ex(3),Z,KK(i,j,:),KKx(i,j,:)) call rget_half_x(ex(3),KKx(i,j,:),HKKx(i,j,:)) enddo enddo return end subroutine calculate_K !-------------------------------------------------------------------- ! this R is indeed x function Eq_Theta(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI) result(gont) implicit none integer,intent(in ):: ex(1:3) real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: beta,W real*8, dimension(ex(1),ex(2),ex(3)),intent(in) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: Cnu,Ck,CUx,DCUx,bDCUx double complex,dimension(ex(3)) :: Cnux,CJx,CJxx,DCJx double complex,dimension(ex(3)) :: fCTheta,CThetax real*8,dimension(ex(3)) :: KK,KKx,betax,Wx double complex :: Theta_rhs,Theta_rhs_o real*8 :: dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=1,ex(2) do i=1,ex(1) CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) fCTheta = dcmplx(RTheta(i,j,:),ITheta(i,j,:)) call cderivs_x(ex(3),R,fCTheta,CThetax) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) call cderivs_x(ex(3),R,Cnu,Cnux) call rderivs_x(ex(3),R,beta(i,j,:),betax) KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:)) call rderivs_x(ex(3),R,KK,KKx) call rderivs_x(ex(3),R,W,Wx) call cderivs_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_x(ex(3),R,CU(i,j,:),CUx) call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) RTheta(i,j,1) = 0.d0 ITheta(i,j,1) = 0.d0 do k=1,ex(3)-1 #if 0 call check_daxiao(beta(i,j,k)) call check_daxiao(CU(i,j,k)) call check_daxiao(bDCU(i,j,k)) call check_daxiao(DCU(i,j,k)) call check_daxiao(CB(i,j,k)) call check_daxiao(DCB(i,j,k)) call check_daxiao(W(i,j,k)) call check_daxiao(CJ(i,j,k)) call check_daxiao(DCJ(i,j,k)) call check_daxiao(bDCB(i,j,k)) call check_daxiao(Cnu(k)) call check_daxiao(Ck(k)) call check_daxiao(fCTheta(k)) call check_daxiao(betax(k)) call check_daxiao(KKx(k)) call check_daxiao(CUx(k)) call check_daxiao(DCUx(k)) call check_daxiao(bDCUx(k)) call check_daxiao(Wx(k)) call check_daxiao(CJx(k)) call check_daxiao(CJxx(k)) call check_daxiao(DCJx(k)) call check_daxiao(Cnux(k)) #endif RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(k),KKx(k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), & DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), & CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),fCTheta(k)) RHS = RHS - CThetax(k) RTheta(i,j,k+1) = dreal(RHS) ITheta(i,j,k+1) = dimag(RHS) enddo write(*,*)RTheta(i,j,:) stop enddo enddo gont = 0 return end function Eq_Theta subroutine check_daxiao(f) implicit none double complex,intent(in) :: f real*8,parameter :: eps=1.d-5 if(cdabs(f)>eps) write(*,*) f return end subroutine check_daxiao subroutine check_factor(T,crho,sigma,R,sst,Rmin) implicit none integer,intent(in) :: sst real*8,intent(in) :: T,crho,sigma,R,Rmin real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf double complex :: Yslm,II,Jr,swtf,ff hgr = R*Rmin/(1.d0-R) tgrho = dtan(crho) tgsigma = dtan(sigma) tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0) ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/2.d0) select case (sst) case (0) z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (1) z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (2) x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (3) x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (4) y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case (5) y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case default write(*,*) "get_null_boundary: not recognized sst = ",sst return end select gt = dacos(z/hgr) gp = datan2(y,x) swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma) if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf) select case (sst) case (0,1) swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2) case (2,3) swtf = II*swtf*dsin(gt) case (4,5) swtf = -II*swtf*dsin(gt) end select write(*,*) dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,0,gt,gp)*swtf**2 return end subroutine check_factor subroutine getdxs(T,crho,sigma,R,betax,KKx,CUx,DCUx,bDCUx,Wx,CJx,CJxx,DCJx,Cnux,CThetax,sst,Rmin) implicit none integer,intent(in) :: sst real*8,intent(in) :: T,crho,sigma,R,Rmin real*8,intent(out) :: betax,KKx,Wx double complex,intent(out) :: CUx,DCUx,bDCUx,CJx,CJxx,DCJx,Cnux,CThetax real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf double complex :: Yslm,II,Jr,swtf,ff double complex :: beta0,C1,C2 integer :: nu,m call initial_null_paramter(beta0,C1,C2,nu,m) II = dcmplx(0.d0,1.d0) hgr = R*Rmin/(1.d0-R) tgrho = dtan(crho) tgsigma = dtan(sigma) tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0) ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/2.d0) select case (sst) case (0) z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (1) z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (2) x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (3) x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (4) y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case (5) y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case default write(*,*) "get_null_boundary: not recognized sst = ",sst return end select gt = dacos(z/hgr) gp = datan2(y,x) swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma) if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf) select case (sst) case (0,1) swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2) case (2,3) swtf = II*swtf*dsin(gt) case (4,5) swtf = -II*swtf*dsin(gt) end select betax = 0.d0 Jr = -(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0/hgr**2& +2.d0*nu*nu*C2/hgr**3-3.d0*II*nu*C2/hgr**4-2.d0*C2/hgr**5 Wx = dreal(Yslm(0,2,m,gt,gp))*dreal(Jr*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin ! Wx = dreal(Jr*cdexp(II*nu*T))*(Rmin+hgr)**2/Rmin KKx = 0.d0 Jr = -2.d0*beta0/hgr/hgr-C1/hgr**3-II*nu*C2/hgr**4-C2/hgr**5 rf = dreal(Jr*cdexp(II*nu*T)) CUx = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf*(Rmin+hgr)**2/Rmin ! CUx = rf*(Rmin+hgr)**2/Rmin DCUx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin ! DCUx = rf*(Rmin+hgr)**2/Rmin bDCUx =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*rf*(Rmin+hgr)**2/Rmin ! bDCUx = rf*(Rmin+hgr)**2/Rmin Jr = -C1/4.d0/hgr**2+C2/4.d0/hgr**4 rf = dreal(Jr*cdexp(II*nu*T)) CJx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin ! CJx = rf*(Rmin+hgr)**2/Rmin Cnux =-dble((2-1)*(2+2))*dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf*(Rmin+hgr)**2/Rmin ! Cnux = rf*(Rmin+hgr)**2/Rmin DCJx = 0.d0 rf = dreal(Jr*II*nu*cdexp(II*nu*T)) CThetax = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**2/Rmin ! CThetax = rf*(Rmin+hgr)**2/Rmin Jr = C1/2.d0/hgr**3-C2/hgr**5 rf = dreal(Jr*cdexp(II*nu*T)) CJxx = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf*(Rmin+hgr)**4/Rmin**2+2.d0*(Rmin+hgr)/Rmin*CJx ! CJxx = rf*(Rmin+hgr)**4/Rmin**2+2.d0*(Rmin+hgr)/Rmin*CJx #if 0 DCUx = DCUx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 CJx = CJx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 CJxx = CJxx*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 CThetax = CThetax*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 #endif return end subroutine getdxs subroutine getndxs(T,crho,sigma,R,beta,KK,CU,bDCU,DCU,CB,DCB,W,CJ,DCJ,bDCB,Cnu,Ck,CTheta,sst,Rmin) implicit none integer,intent(in) :: sst real*8,intent(in) :: T,crho,sigma,R,Rmin real*8,intent(out) :: beta,KK,W double complex,intent(out) :: CU,bDCU,DCU,CB,DCB,CJ,DCJ,bDCB,Cnu,Ck,CTheta real*8 ::x,y,z,hgr,gt,gp,tgrho,tgsigma,tc,ts,rf double complex :: Yslm,II,Jr,swtf,ff double complex :: beta0,C1,C2 integer :: nu,m call initial_null_paramter(beta0,C1,C2,nu,m) II = dcmplx(0.d0,1.d0) hgr = R*Rmin/(1.d0-R) tgrho = dtan(crho) tgsigma = dtan(sigma) tc = dsqrt((1.d0-dsin(crho)*dsin(sigma))/2.d0) ts = dsqrt((1.d0+dsin(crho)*dsin(sigma))/2.d0) select case (sst) case (0) z = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (1) z = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = z*tgrho y = z*tgsigma case (2) x = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (3) x = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) y = x*tgrho z = x*tgsigma case (4) y = hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case (5) y = -hgr/dsqrt(1+tgrho*tgrho+tgsigma*tgsigma) x = y*tgrho z = y*tgsigma case default write(*,*) "get_null_boundary: not recognized sst = ",sst return end select gt = dacos(z/hgr) gp = datan2(y,x) swtf = 2.d0*tc*ts*(ts+II*tc)/dcos(sigma) if(sst==1 .or. sst==3 .or. sst==4) swtf = dconjg(swtf) select case (sst) case (0,1) swtf = swtf/(dcos(gp)+II*dcos(gt)*dsin(gp))*(dcos(gt)**2+dsin(gt)**2*dcos(gp)**2) case (2,3) swtf = II*swtf*dsin(gt) case (4,5) swtf = -II*swtf*dsin(gt) end select beta = dreal(Yslm(0,2,m,gt,gp))*dreal(beta0*cdexp(II*nu*T)) ! beta = dreal(beta0*cdexp(II*nu*T)) CB = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*dreal(beta0*cdexp(II*nu*T)) ! CB = dreal(beta0*cdexp(II*nu*T)) DCB = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*dreal(beta0*cdexp(II*nu*T)) ! DCB = dreal(beta0*cdexp(II*nu*T)) bDCB =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*dreal(beta0*cdexp(II*nu*T)) ! bDCB = dreal(beta0*cdexp(II*nu*T)) Jr = (2.4d1*II*nu*beta0-3.d0*nu*nu*C1+nu**4*C2)/6.d0+(3.d0*II*nu*C1-6.d0*beta0-II*nu**3*C2)/3.d0/hgr& -nu*nu*C2/hgr/hgr+II*nu*C2/hgr**3+C2/2.d0/hgr**4 W = dreal(Yslm(0,2,m,gt,gp))*dreal(Jr*cdexp(II*nu*T)) ! W = dreal(Jr*cdexp(II*nu*T)) Jr = (2.4d1*beta0+3.d0*II*nu*C1-II*nu**3*C2)/3.6d1+C1/4.d0/hgr-C2/1.2d1/hgr**3 rf = dreal(Jr*cdexp(II*nu*T)) CJ = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf ! CJ = rf DCJ = 0.d0 Cnu =-dsqrt(dble((2+2)*(2-2+1)*(2-1)*2*(2+1)*(2+2)))*Yslm(1,2,m,gt,gp)*swtf*rf ! Cnu = rf KK = dsqrt(1.d0+cdabs(CJ)**2) Ck = 0.d0 rf = dreal(Jr*II*nu*cdexp(II*nu*T)) CTheta = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf ! CTheta = rf Jr = (-2.4d1*II*nu*beta0+3.d0*nu*nu*C1-nu**4*C2)/36.d0+2.d0*beta0/hgr& +C1/2.d0/hgr/hgr+II*nu*C2/3.d0/hgr**3+C2/4.d0/hgr**4 rf = dreal(Jr*cdexp(II*nu*T)) CU = dsqrt(dble(2*(2+1)))*Yslm(1,2,m,gt,gp)*swtf*rf ! CU = rf DCU = dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2*rf ! DCU = rf bDCU =-dble(2*(2+1))*Yslm(0,2,m,gt,gp)*rf ! bDCU = rf #if 0 DCU = DCU*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 DCB = DCB*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 CTheta = CTheta*dsqrt(dble((2-1)*2*(2+1)*(2+2)))*Yslm(2,2,m,gt,gp)*swtf**2 #endif return end subroutine getndxs !-------------------------------------------------------------------- ! this R is indeed x function Eq_Theta_2(ex,crho,sigma,R,RJ,IJ,RU,IU,beta,RB,IB, & Rnu,Inu,Rk,Ik,RTheta,ITheta,W,Rmin, & qlR1,qlR2,qlI1,qlI2,quR1,quR2,quI1,quI2,gR,gI, & T,sst) result(gont) implicit none integer,intent(in ):: ex(1:3),sst real*8,intent(in),dimension(ex(1))::crho real*8,intent(in),dimension(ex(2))::sigma real*8,intent(in),dimension(ex(3))::R real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: beta,W real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RJ,IJ,RU,IU,RB,IB,Rnu,Inu,Rk,Ik real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: RTheta,ITheta real*8,intent(in) :: Rmin,T real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: qlR1,qlR2,qlI1,qlI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: quR1,quR2,quI1,quI2 real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gR,gI ! gont = 0: success; gont = 1: something wrong integer::gont double complex,dimension(ex(1),ex(2),ex(3)) :: CU,DCU,bDCU,CB,DCB,bDCB,CJ,DCJ double complex :: CTheta0,CTheta,CTheta1,RHS integer :: i,j,k,RK4 double complex,dimension(ex(3)) :: Cnu,Ck,HCnu,HCk,HCU,HDCU,CUx,HCUx,DCUx,HDCUx,HbDCU,bDCUx,HbDCUx double complex,dimension(ex(3)) :: Cnux,HCnux,HCJ,HDCJ,CJx,HCJx,CJxx,HCJxx,DCJx,HDCJx,HCB,HDCB,HbDCB double complex,dimension(ex(3)) :: fCTheta,CThetax real*8,dimension(ex(3)) :: KK,KKx,HKK,HKKx,Hbeta,betax,Hbetax,HW,Wx,HWx double complex :: Theta_rhs,Theta_rhs_o real*8 :: dR !!! sanity check dR = sum(RJ)+sum(IJ)+sum(RU)+sum(IU)+sum(beta)+sum(RB)+sum(IB) + & sum(Rnu)+sum(Inu)+sum(Rk)+sum(Ik)+sum(RTheta)+sum(ITheta) if(dR.ne.dR) then if(sum(RJ).ne.sum(RJ))write(*,*)"NullEvol_Theta: find NaN in RJ" if(sum(IJ).ne.sum(IJ))write(*,*)"NullEvol_Theta: find NaN in IJ" if(sum(RU).ne.sum(RU))write(*,*)"NullEvol_Theta: find NaN in RU" if(sum(IU).ne.sum(IU))write(*,*)"NullEvol_Theta: find NaN in IU" if(sum(beta).ne.sum(beta))write(*,*)"NullEvol_Theta: find NaN in beta" if(sum(RB).ne.sum(RB))write(*,*)"NullEvol_Theta: find NaN in RB" if(sum(IB).ne.sum(IB))write(*,*)"NullEvol_Theta: find NaN in IB" if(sum(Rnu).ne.sum(Rnu))write(*,*)"NullEvol_Theta: find NaN in Rnu" if(sum(Inu).ne.sum(Inu))write(*,*)"NullEvol_Theta: find NaN in Inu" if(sum(Rk).ne.sum(Rk))write(*,*)"NullEvol_Theta: find NaN in Rk" if(sum(Ik).ne.sum(Ik))write(*,*)"NullEvol_Theta: find NaN in Ik" if(sum(RTheta).ne.sum(RTheta))write(*,*)"NullEvol_Theta: find NaN in RTheta" if(sum(ITheta).ne.sum(ITheta))write(*,*)"NullEvol_Theta: find NaN in ITheta" gont = 1 return endif dR = R(2) - R(1) CU = dcmplx(RU,IU) CB = dcmplx(RB,IB) CJ = dcmplx(RJ,IJ) do k=1,ex(3) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),DCU(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CU(:,:,k),bDCU(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),DCB(:,:,k),1,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CB(:,:,k),bDCB(:,:,k),1,-1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) call derivs_eth(ex(1:2),crho,sigma,CJ(:,:,k),DCJ(:,:,k),2,1, & quR1(:,:,k),quR2(:,:,k),quI1(:,:,k),quI2(:,:,k),gR(:,:,k),gI(:,:,k)) enddo do j=ghost_width+1,ex(2)-ghost_width do i=ghost_width+1,ex(1)-ghost_width CTheta0 = dcmplx(RTheta(i,j,1),ITheta(i,j,1)) fCTheta = dcmplx(RTheta(i,j,:),ITheta(i,j,:)) call cderivs_x(ex(3),R,fCTheta,CThetax) Cnu = dcmplx(Rnu(i,j,:),Inu(i,j,:)) Ck = dcmplx(Rk(i,j,:),Ik(i,j,:)) call cget_half_x(ex(3),CB(i,j,:),HCB) call cget_half_x(ex(3),DCB(i,j,:),HDCB) call cget_half_x(ex(3),bDCB(i,j,:),HbDCB) call cget_half_x(ex(3),Cnu,HCnu) call cderivs_x(ex(3),R,Cnu,Cnux) call cget_half_x(ex(3),Cnux,HCnux) call cget_half_x(ex(3),Ck,HCk) call rget_half_x(ex(3),beta(i,j,:),Hbeta) call rderivs_x(ex(3),R,beta(i,j,:),betax) call rget_half_x(ex(3),betax,Hbetax) KK = dsqrt(1.d0+RJ(i,j,:)*RJ(i,j,:)+IJ(i,j,:)*IJ(i,j,:)) call rget_half_x(ex(3),KK,HKK) call rderivs_x(ex(3),R,KK,KKx) call rget_half_x(ex(3),KKx,HKKx) call rderivs_x(ex(3),R,W,Wx) call rget_half_x(ex(3),Wx,HWx) call rget_half_x(ex(3),W(i,j,:),HW) call cget_half_x(ex(3),CU(i,j,:),HCU) call cderivs_x(ex(3),R,DCU(i,j,:),DCUx) call cderivs_x(ex(3),R,CU(i,j,:),CUx) call cget_half_x(ex(3),DCUx,HDCUx) call cget_half_x(ex(3),CUx,HCUx) call cget_half_x(ex(3),DCU(i,j,:),HDCU) call cderivs_x(ex(3),R,bDCU(i,j,:),bDCUx) call cget_half_x(ex(3),bDCUx,HbDCUx) call cget_half_x(ex(3),bDCU(i,j,:),HbDCU) call cderivs_x(ex(3),R,CJ(i,j,:),CJx) call cdderivs_x(ex(3),R,CJ(i,j,:),CJxx) call cget_half_x(ex(3),CJx,HCJx) call cget_half_x(ex(3),CJxx,HCJxx) call cderivs_x(ex(3),R,DCJ(i,j,:),DCJx) call cget_half_x(ex(3),DCJx,HDCJx) RTheta(i,j,1) = 0.d0 ITheta(i,j,1) = 0.d0 do k=1,ex(3)-1 ! call getndxs(T,crho(i),sigma(j),R(k),beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k), & ! CB(i,j,k),DCB(i,j,k),W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k),sst,Rmin) ! call getdxs(T,crho(i),sigma(j),R(k),betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k), & ! Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k),sst,Rmin) RHS = Theta_rhs(R(k),Rmin,beta(i,j,k),betax(k),KK(k),KKx(k),CU(i,j,k),CUx(k),DCUx(k),bDCU(i,j,k),bDCUx(k), & DCU(i,j,k),CB(i,j,k),DCB(i,j,k),W(i,j,k),Wx(k),CJ(i,j,k),DCJ(i,j,k), & CJx(k),CJxx(k),DCJx(k),bDCB(i,j,k),Cnu(k),Cnux(k),Ck(k),fCTheta(k)) RHS = RHS - CThetax(k) #if 0 if(cdabs(RHS)>1.d-9)then #if 0 write(*,*)beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k),CB(i,j,k),DCB(i,j,k) write(*,*)W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k) call getndxs(T,crho(i),sigma(j),R(k),beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k), & CB(i,j,k),DCB(i,j,k),W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k),sst,Rmin) write(*,*)"VS" write(*,*)beta(i,j,k),KK(k),CU(i,j,k),bDCU(i,j,k),DCU(i,j,k),CB(i,j,k),DCB(i,j,k) write(*,*)W(i,j,k),CJ(i,j,k),DCJ(i,j,k),bDCB(i,j,k),Cnu(k),Ck(k),fCTheta(k) #endif write(*,*)betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k) write(*,*)Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k) call getdxs(T,crho(i),sigma(j),R(k),betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k), & Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k),sst,Rmin) write(*,*)"VS" write(*,*)betax(k),KKx(k),CUx(k),DCUx(k),bDCUx(k) write(*,*)Wx(k),CJx(k),CJxx(k),DCJx(k),Cnux(k),CThetax(k) ! write(*,*)RHS ! call check_factor(T,crho(i),sigma(j),R(k),sst,Rmin) stop endif #endif RTheta(i,j,k+1) = dreal(RHS) ITheta(i,j,k+1) = dimag(RHS) enddo enddo enddo gont = 0 return end function Eq_Theta_2