#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.i2 .and. j>2.and.i3 .and. j>3.and.i4 .and. j>4.and.ieps) 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