c ****************************************************************** c * * c * PROGRAM TO SOLVE THE TIME-DEPENDENT GINZBURG-LANDAU * c * EQUATIONS ON RECTANGULAR DOMAINS WITH RECTANGULAR HOLES * c * * c ****************************************************************** PROGRAM TDGL c **************** Initialization of variables ********************* PARAMETER (IMAX=100,JMAX=100) INTEGER Nx ! number of cells in x INTEGER Ny ! number of cells in x REAL Ax ! mesh spacing in x REAL Ay ! mesh spacing in y REAL dt ! time step REAL k ! kappa ( = lambda / xi ) REAL Eta ! eta REAL T ! temperature ( scaled with Tc ) REAL He ! applied external magnetic field REAL Ho ! initial applied external magnetic field REAL Hi ! magnetic field in the interior hole REAL dHe ! applied external magnetic field increment REAL Ruido ! Standard deviation of random noise REAL Eo ! Parameter for random noise INTEGER TotSteps! number of time steps INTEGER Step ! step counter INTEGER OutSteps! number of steps between outputs INTEGER OutStp2 ! number of steps between run-time information INTEGER HeSteps ! number of steps between applied field increments INTEGER idum ! auxiliary INTEGER i,j ! auxiliary CHARACTER*36 Label ! auxiliary INTEGER iHe ! auxiliary for naming restart files COMPLEX F(IMAX+1,JMAX+1) ! order parameter COMPLEX Ux(IMAX,JMAX+1) ! link variable in x COMPLEX Uy(IMAX+1,JMAX) ! link variable in y COMPLEX bloop(IMAX,JMAX) ! exp(-i a_x a_y B_z) c The following three arrays are auxiliary to avoid copies COMPLEX G(IMAX+1,JMAX+1) ! order parameter COMPLEX Vx(IMAX,JMAX+1) ! link variable in x COMPLEX Vy(IMAX+1,JMAX) ! link variable in y c LOGICAL Bulk(IMAX,JMAX) ! bulk cell indicator INTEGER iswnod (IMAX+1,JMAX+1) ! vertex indicator INTEGER iswlx (IMAX,JMAX+1) ! x-link indicator INTEGER iswly (IMAX+1,JMAX) ! y-link indicator REAL rAx2 ! auxiliar REAL rAy2 ! auxiliar REAL rAxAy ! auxiliar REAL k2 ! auxiliar REAL rEta ! auxiliar REAL T1 ! auxiliar REAL rT1 ! auxiliar REAL k2rT12 ! auxiliar REAL AxAy ! auxiliar REAL AxAyHe ! auxiliar REAL AxAyHi ! auxiliar REAL AxAydHi ! auxiliar COMPLEX EiAxAyHe! auxiliar, approx. exp(-i*ax*ay*He) COMPLEX EiAxAyHi! auxiliar, approx. exp(-i*ax*ay*Hi) REAL rotA ! auxiliar COMPLEX EiMAxAydHi ! var. aux., approx. exp(-i*M*ax*ay*dHi) COMPLEX EiMAxAyHi ! var. aux., approx. exp(-i*M*ax*ay*dHi) REAL FFF ! auxiliar REAL FreeEnergy ! ! REAL Bz ! mean field REAL Magnetization ! ! COMPLEX GaussDev! complex numbers generator (Gaussian distrib.) INTEGER d ! auxiliar REAL Flux ! total fluxoid REAL Flux1 ! internal fluxoid REAL PATHx1 ! REAL PATHx2 ! REAL PATHy1 ! REAL PATHy2 ! circulations COMPLEX PTHx1 ! COMPLEX PTHx2 ! COMPLEX PTHy1 ! COMPLEX PTHy2 ! circulations REAL MinF ! REAL MinF1 ! INTEGER N ! Hole N boundary INTEGER S ! Hole S boundary INTEGER E ! Hole E boundary INTEGER W ! Hole W boundary INTEGER M ! auxiliar REAL rM ! auxiliar c ****************************************************************** iHe = 9000000 c ********************** Input data ************************** OPEN (8,FILE='tdgl.ini',FORM='FORMATTED') READ (8,*) Label,Nx WRITE (6,*) Label,'Nx: ',Nx READ (8,*) Label,Ny WRITE (6,*) Label,'Ny: ',Ny READ (8,*) Label,Ax WRITE (6,*) Label,'Ax: ',Ax READ (8,*) Label,Ay WRITE (6,*) Label,'Ay: ',Ay READ (8,*) Label,dt WRITE (6,*) Label,'dt: ',dt READ (8,*) Label,k WRITE (6,*) Label,'k: ',k READ (8,*) Label,Eta WRITE (6,*) Label,'Eta: ',Eta READ (8,*) Label,T WRITE (6,*) Label,'T: ',T READ (8,*) Label,He WRITE (6,*) Label,'He: ',He READ (8,*) Label,dHe WRITE (6,*) Label,'dHe: ',dHe READ (8,*) Label,Eo WRITE (6,*) Label,'Eo: ',Eo READ (8,*) Label,TotSteps WRITE (6,*) Label,'TotSteps: ',TotSteps READ (8,*) Label,OutSteps WRITE (6,*) Label,'OutSteps: ',OutSteps READ (8,*) Label,OutStp2 WRITE (6,*) Label,'OutStp2: ',OutStp2 READ (8,*) Label,HeSteps WRITE (6,*) Label,'HeSteps: ',HeSteps READ (8,*) Label,d WRITE (6,*) Label,'d: ',d READ (8,*) Label,N WRITE (6,*) Label,'N: ',N READ (8,*) Label,S WRITE (6,*) Label,'S: ',S READ (8,*) Label,E WRITE (6,*) Label,'E: ',E READ (8,*) Label,W WRITE (6,*) Label,'W: ',W CLOSE (8) c ****************************************************************** c **************** Initialization ********************************** M=(N-S)*(E-W)+2*(N-S)+2*(E-W) rM=1.0 / M Ho=He Hi=0.0 k2=k*k rEta=1/Eta T1=1-T rT1=1/T1 k2rT12=(k*k)/(T1*T1) rAx2=1/(Ax*Ax) rAy2=1/(Ay*Ay) rAxAy=1/(Ax*Ay) AxAy=Ax*Ay AxAyHe=Ax*Ay*He EiAxAyHe=COS(AxAyHe)-(0,1)*SIN(AxAyHe) Ruido=dt*SQRT(0.5235987755983*Eo*dt*T) idum=0 c ****************************************************************** c ********************* Prepare output ***************************** OPEN (13,FILE='s.dat',STATUS='NEW',FORM='FORMATTED') WRITE (13,*) '# Step',' He',' Hi', * ' Bz',' Magnetization', * ' FreeEnergy',' Flux', * ' Flux1' * ,' MinF',' MinF1' CLOSE (13) OPEN (15,FILE='c.dat',STATUS='NEW',FORM='FORMATTED') WRITE (15,*) ' He',' Bz',' Flux',' Hi' CLOSE (15) c ****************************************************************** c *************************** Restart ****************************** OPEN (11,FILE='tdgl.ree',STATUS='OLD',FORM='FORMATTED',ERR=10) WRITE (6,*) 'Reading the restart file: tdgl.ree... ' READ (11,*) Hi READ (11,*) ((F(i,j),i=1,Nx+1),j=1,Ny+1) READ (11,*) ((Ux(i,j),i=1,Nx),j=1,Ny+1) READ (11,*) ((Uy(i,j),i=1,Nx+1),j=1,Ny) CLOSE (11, STATUS='KEEP') GO TO 20 10 WRITE (6,*) 'Initializing to a perfect Meissner state... ' DO i=1,Nx+1 DO j=1,Ny+1 F(i,j)=1 END DO END DO DO i=1,Nx DO j=1,Ny+1 Ux(i,j)=1 END DO END DO DO i=1,Nx+1 DO j=1,Ny Uy(i,j)=1 END DO END DO 20 DO i=1,Nx DO j=1,Ny Bulk(i,j)=.TRUE. END DO END DO DO i=1,Nx+1 DO j=1,Ny+1 iswnod(i,j)=0 END DO END DO DO i=1,Nx DO j=1,Ny+1 iswlx(i,j)=0 END DO END DO DO i=1,Nx+1 DO j=1,Ny iswly(i,j)=0 END DO END DO DO i=W,E-1 DO j=S,N-1 Bulk(i,j)=.FALSE. END DO END DO DO i=1,Nx DO j=1,Ny if (Bulk(i,j)) then iswnod(i,j)=iswnod(i,j)+1 iswnod(i+1,j)=iswnod(i+1,j)+1 iswnod(i,j+1)=iswnod(i,j+1)+1 iswnod(i+1,j+1)=iswnod(i+1,j+1)+1 iswlx(i,j)=iswlx(i,j)+1 iswlx(i,j+1)=iswlx(i,j+1)+1 iswly(i,j)=iswly(i,j)+1 iswly(i+1,j)=iswly(i+1,j)+1 end if END DO END DO DO i=W+1,E-1 DO j=S+1,N-1 F(i,j)=0 END DO END DO AxAyHi=AxAy*Hi EiMAxAyHi=COS(M*AxAyHi)-(0,1)*SIN(M*AxAyHi) c ****************************************************************** WRITE (6,*) 'Starting Calculations... ' c *********************** Main loop ******************************** DO Step=2,TotSteps,2 c ********* Time integration step of TDGL equations **************** DO i=2,Nx DO j=2,Ny IF(iswnod(i,j).ge.3) THEN G(i,j)=F(i,j)+dt*rEta * *( rAx2*(Ux(i,j)*F(i+1,j) * -2*F(i,j) * +CONJG(Ux(i-1,j))*F(i-1,j)) * +rAy2*(Uy(i,j)*F(i,j+1) * -2*F(i,j) * +CONJG(Uy(i,j-1))*F(i,j-1)) * -T1*(F(i,j)*CONJG(F(i,j))-1 )*F(i,j)) * +GaussDev(idum)*Ruido END IF END DO END DO DO i=1,Nx DO j=1,Ny if (Bulk(i,j)) * bloop(i,j) = Ux(i,j)*Uy(i+1,j)* * CONJG(Ux(i,j+1))*CONJG(Uy(i,j)) END DO END DO DO i=1,Nx DO j=2,Ny IF(iswlx(i,j).eq.2) THEN Vx(i,j)=Ux(i,j)+dt* * ((0,-1)*T1*Ux(i,j)*AIMAG(CONJG(F(i,j))*Ux(i,j)*F(i+1,j)) * -k2*rAy2*Ux(i,j)* * (bloop(i,j)*CONJG(bloop(i,j-1))-1)) c Vx(i,j)=Vx(i,j)/CABS(Vx(i,j)) END IF END DO END DO DO i=2,Nx DO j=1,Ny IF(iswly(i,j).eq.2) THEN Vy(i,j)=Uy(i,j)+dt* * ((0,-1)*T1*Uy(i,j)*AIMAG(CONJG(F(i,j))*Uy(i,j)*F(i,j+1)) * -k2*rAx2*Uy(i,j)* * (bloop(i-1,j)*CONJG(bloop(i,j))-1)) c Vy(i,j)=Vy(i,j)/CABS(Vy(i,j)) END IF END DO END DO DO j=2,Ny G(1,j)=G(2,j)*Vx(1,j) G(Nx+1,j)=G(Nx,j)*CONJG(Vx(Nx,j)) END DO DO i=2,Nx G(i,1)=G(i,2)*Vy(i,1) G(i,Ny+1)=G(i,Ny)*CONJG(Vy(i,Ny)) END DO c cosmetic G(1,1)=(G(1,2)+G(2,1))*0.50 G(1,Ny+1)=(G(1,Ny)+G(2,Ny+1))*0.50 G(Nx+1,1)=(G(Nx+1,2)+G(Nx,1))*0.50 G(Nx+1,Ny+1)=(G(Nx,Ny+1)+G(Nx+1,Ny))*0.50 DO i=2,Nx-1 Vx(i,1)=EiAxAyHe* * (CONJG(Vy(i+1,1))*Vx(i,2)*Vy(i,1)) c Vx(i,1)=Vx(i,1)/CABS(Vx(i,1)) Vx(i,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Vy(i,Ny))*Vx(i,Ny)*Vy(i+1,Ny)) c Vx(i,Ny+1)=Vx(i,Ny+1)/CABS(Vx(i,Ny+1)) END DO Vx(1,1)=EiAxAyHe* * (CONJG(Vy(2,1))*Vx(1,2)*Uy(1,1)) c Vx(1,1)=Vx(1,1)/CABS(Vx(1,1)) Vx(1,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Uy(1,Ny))*Vx(1,Ny)*Vy(2,Ny)) c Vx(1,Ny+1)=Vx(1,Ny+1)/CABS(Vx(1,Ny+1)) Vx(Nx,1)=EiAxAyHe* * (CONJG(Uy(Nx+1,1))*Vx(Nx,2)*Vy(Nx,1)) c Vx(Nx,1)=Vx(Nx,1)/CABS(Vx(Nx,1)) Vx(Nx,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Vy(Nx,Ny))*Vx(Nx,Ny)*Uy(Nx+1,Ny)) c Vx(Nx,Ny+1)=Vx(Nx,Ny+1)/CABS(Vx(Nx,Ny+1)) DO j=1,Ny Vy(1,j)=CONJG(EiAxAyHe)* * (Vx(1,j)*Vy(2,j)*CONJG(Vx(1,j+1))) c Vy(1,j)=Vy(1,j)/CABS(Vy(1,j)) Vy(Nx+1,j)=EiAxAyHe* * (Vx(Nx,j+1)*Vy(Nx,j)*CONJG(Vx(Nx,j))) c Vy(Nx+1,j)=Vy(Nx+1,j)/CABS(Vy(Nx+1,j)) END DO if (E.gt.W.and.N.gt.S) then DO j=S+1,N-1 G(E,j)=G(E+1,j)*Vx(E,j) G(W,j)=G(W-1,j)*CONJG(Vx(W-1,j)) END DO DO i=W+1,E-1 G(i,N)=G(i,N+1)*Vy(i,N) G(i,S)=G(i,S-1)*CONJG(Vy(i,S-1)) END DO c G(W,N)=(G(W+1,N)+G(W-1,N)+G(W,N+1)+G(W,N-1))*0.250 c G(W,S)=(G(W+1,S)+G(W-1,S)+G(W,S+1)+G(W,S-1))*0.250 c G(E,N)=(G(E+1,N)+G(E-1,N)+G(E,N+1)+G(E,N-1))*0.250 c G(E,S)=(G(E+1,S)+G(E-1,S)+G(E,S+1)+G(E,S-1))*0.250 PTHx1=Vy(E,S-1)*CONJG(Vy(W,S-1)) PTHx2=Vy(E,N)*CONJG(Vy(W,N)) PTHy1=Vx(W-1,S)*CONJG(Vx(W-1,N)) PTHy2=Vx(E,S)*CONJG(Vx(E,N)) DO i=W,E-1 PTHx1=PTHx1*Vx(i,S-1) PTHx1=PTHx1/CABS(PTHx1) PTHx2=PTHx2*CONJG(Vx(i,N+1)) PTHx2=PTHx2/CABS(PTHx2) END DO DO j=S,N-1 PTHy1=PTHy1*CONJG(Vy(W-1,j)) PTHy1=PTHy1/CABS(PTHy1) PTHy2=PTHy2*Vy(E+1,j) PTHy2=PTHy2/CABS(PTHy2) END DO EiMAxAydHi=PTHx1*PTHx2*PTHy1*PTHy2*CONJG(EiMAxAyHi) EiMAxAydHi=EiMAxAydHi/CABS(EiMAxAydHi) EiMAxAyHi=PTHx1*PTHx2*PTHy1*PTHy2 EiMAxAyHi=EiMAxAyHi/CABS(EiMAxAyHi) AxAydHi=-ATAN2(AIMAG(EiMAxAydHi),REAL(EiMAxAydHi))*rM AxAyHi=AxAyHi+AxAydHi EiAxAyHi=COS(AxAyHi)-(0,1)*SIN(AxAyHi) DO i=W,E-1 Vx(i,N)=EiAxAyHi* * (CONJG(Vy(i+1,N))*Vx(i,N+1)*Vy(i,N)) c Vx(i,N)=Vx(i,N)/CABS(Vx(i,N)) Vx(i,S)=CONJG(EiAxAyHi)* * (CONJG(Vy(i,S-1))*Vx(i,S-1)*Vy(i+1,S-1)) c Vx(i,S)=Vx(i,S)/CABS(Vx(i,S)) END DO DO j=S,N-1 Vy(E,j)=CONJG(EiAxAyHi)* * (Vx(E,j)*Vy(E+1,j)*CONJG(Vx(E,j+1))) c Vy(E,j)=Vy(E,j)/CABS(Vy(E,j)) Vy(W,j)=EiAxAyHi* * (Vx(W-1,j+1)*Vy(W-1,j)*CONJG(Vx(W-1,j))) c Vy(W,j)=Vy(W,j)/CABS(Vy(W,j)) END DO end if c ****************************************************************** DO i=2,Nx DO j=2,Ny IF (iswnod(i,j).ge.3) THEN F(i,j)=G(i,j)+dt*rEta * *( rAx2*(Vx(i,j)*G(i+1,j) * -2*G(i,j) * +CONJG(Vx(i-1,j))*G(i-1,j)) * +rAy2*(Vy(i,j)*G(i,j+1) * -2*G(i,j) * +CONJG(Vy(i,j-1))*G(i,j-1)) * -T1*(G(i,j)*CONJG(G(i,j))-1 )*G(i,j)) * +GaussDev(idum)*Ruido END IF END DO END DO DO i=1,Nx DO j=1,Ny if (Bulk(i,j)) * bloop(i,j) = Vx(i,j)*Vy(i+1,j)* * CONJG(Vx(i,j+1))*CONJG(Vy(i,j)) END DO END DO DO i=1,Nx DO j=2,Ny IF(iswlx(i,j).eq.2) THEN Ux(i,j)=Vx(i,j)+dt* * ((0,-1)*T1*Vx(i,j)*AIMAG(CONJG(G(i,j))*Vx(i,j)*G(i+1,j)) * -k2*rAy2*Vx(i,j)* * (bloop(i,j)*CONJG(bloop(i,j-1))-1)) Ux(i,j)=Ux(i,j)/CABS(Ux(i,j)) END IF END DO END DO DO i=2,Nx DO j=1,Ny IF(iswly(i,j).eq.2) THEN Uy(i,j)=Vy(i,j)+dt* * ((0,-1)*T1*Vy(i,j)*AIMAG(CONJG(G(i,j))*Vy(i,j)*G(i,j+1)) * -k2*rAx2*Vy(i,j)* * (bloop(i-1,j)*CONJG(bloop(i,j))-1)) Uy(i,j)=Uy(i,j)/CABS(Uy(i,j)) END IF END DO END DO DO j=2,Ny F(1,j)=F(2,j)*Ux(1,j) F(Nx+1,j)=F(Nx,j)*CONJG(Ux(Nx,j)) END DO DO i=2,Nx F(i,1)=F(i,2)*Uy(i,1) F(i,Ny+1)=F(i,Ny)*CONJG(Uy(i,Ny)) END DO c cosmetic F(1,1)=(F(1,2)+F(2,1))*0.50 F(1,Ny+1)=(F(1,Ny)+F(2,Ny+1))*0.50 F(Nx+1,1)=(F(Nx+1,2)+F(Nx,1))*0.50 F(Nx+1,Ny+1)=(F(Nx,Ny+1)+F(Nx+1,Ny))*0.50 DO i=2,Nx-1 Ux(i,1)=EiAxAyHe* * (CONJG(Uy(i+1,1))*Ux(i,2)*Uy(i,1)) Ux(i,1)=Ux(i,1)/CABS(Ux(i,1)) Ux(i,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Uy(i,Ny))*Ux(i,Ny)*Uy(i+1,Ny)) Ux(i,Ny+1)=Ux(i,Ny+1)/CABS(Ux(i,Ny+1)) END DO Ux(1,1)=EiAxAyHe* * (CONJG(Uy(2,1))*Ux(1,2)*Vy(1,1)) Ux(1,1)=Ux(1,1)/CABS(Ux(1,1)) Ux(1,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Vy(1,Ny))*Ux(1,Ny)*Uy(2,Ny)) Ux(1,Ny+1)=Ux(1,Ny+1)/CABS(Ux(1,Ny+1)) Ux(Nx,1)=EiAxAyHe* * (CONJG(Vy(Nx+1,1))*Ux(Nx,2)*Uy(Nx,1)) Ux(Nx,1)=Ux(Nx,1)/CABS(Ux(Nx,1)) Ux(Nx,Ny+1)=CONJG(EiAxAyHe)* * (CONJG(Uy(Nx,Ny))*Ux(Nx,Ny)*Vy(Nx+1,Ny)) Ux(Nx,Ny+1)=Ux(Nx,Ny+1)/CABS(Ux(Nx,Ny+1)) DO j=1,Ny Uy(1,j)=CONJG(EiAxAyHe)* * (Ux(1,j)*Uy(2,j)*CONJG(Ux(1,j+1))) Uy(1,j)=Uy(1,j)/CABS(Uy(1,j)) Uy(Nx+1,j)=EiAxAyHe* * (Ux(Nx,j+1)*Uy(Nx,j)*CONJG(Ux(Nx,j))) Uy(Nx+1,j)=Uy(Nx+1,j)/CABS(Uy(Nx+1,j)) END DO if (E.gt.W.and.N.gt.S) then DO j=S+1,N-1 F(E,j)=F(E+1,j)*Ux(E,j) F(W,j)=F(W-1,j)*CONJG(Ux(W-1,j)) END DO DO i=W+1,E-1 F(i,N)=F(i,N+1)*Uy(i,N) F(i,S)=F(i,S-1)*CONJG(Uy(i,S-1)) END DO c F(W,N)=(F(W+1,N)+F(W-1,N)+F(W,N+1)+F(W,N-1))*0.250 c F(W,S)=(F(W+1,S)+F(W-1,S)+F(W,S+1)+F(W,S-1))*0.250 c F(E,N)=(F(E+1,N)+F(E-1,N)+F(E,N+1)+F(E,N-1))*0.250 c F(E,S)=(F(E+1,S)+F(E-1,S)+F(E,S+1)+F(E,S-1))*0.250 PTHx1=Uy(E,S-1)*CONJG(Uy(W,S-1)) PTHx2=Uy(E,N)*CONJG(Uy(W,N)) PTHy1=Ux(W-1,S)*CONJG(Ux(W-1,N)) PTHy2=Ux(E,S)*CONJG(Ux(E,N)) DO i=W,E-1 PTHx1=PTHx1*Ux(i,S-1) PTHx1=PTHx1/CABS(PTHx1) PTHx2=PTHx2*CONJG(Ux(i,N+1)) PTHx2=PTHx2/CABS(PTHx2) END DO DO j=S,N-1 PTHy1=PTHy1*CONJG(Uy(W-1,j)) PTHy1=PTHy1/CABS(PTHy1) PTHy2=PTHy2*Uy(E+1,j) PTHy2=PTHy2/CABS(PTHy2) END DO EiMAxAydHi=PTHx1*PTHx2*PTHy1*PTHy2*CONJG(EiMAxAyHi) EiMAxAydHi=EiMAxAydHi/CABS(EiMAxAydHi) EiMAxAyHi=PTHx1*PTHx2*PTHy1*PTHy2 EiMAxAyHi=EiMAxAyHi/CABS(EiMAxAyHi) AxAydHi=-ATAN2(AIMAG(EiMAxAydHi),REAL(EiMAxAydHi))*rM AxAyHi=AxAyHi+AxAydHi EiAxAyHi=COS(AxAyHi)-(0,1)*SIN(AxAyHi) DO i=W,E-1 Ux(i,N)=EiAxAyHi* * (CONJG(Uy(i+1,N))*Ux(i,N+1)*Uy(i,N)) Ux(i,N)=Ux(i,N)/CABS(Ux(i,N)) Ux(i,S)=CONJG(EiAxAyHi)* * (CONJG(Uy(i,S-1))*Ux(i,S-1)*Uy(i+1,S-1)) Ux(i,S)=Ux(i,S)/CABS(Ux(i,S)) END DO DO j=S,N-1 Uy(E,j)=CONJG(EiAxAyHi)* * (Ux(E,j)*Uy(E+1,j)*CONJG(Ux(E,j+1))) Uy(E,j)=Uy(E,j)/CABS(Uy(E,j)) Uy(W,j)=EiAxAyHi* * (Ux(W-1,j+1)*Uy(W-1,j)*CONJG(Ux(W-1,j))) Uy(W,j)=Uy(W,j)/CABS(Uy(W,j)) END DO end if c ****************************************************************** c *********************** Miscellaneous outputs ******************** IF (MOD(Step,OutStp2) .LE. 1) THEN Hi=rAxAy*AxAyHi FreeEnergy=0 Bz=0 MaxJsy=0 c condensation part DO i=1,Nx+1 DO j=1,Ny+1 fff = CABS(F(i,j)) FreeEnergy = FreeEnergy * + iswnod(i,j)*0.25d0*(0.5d0*fff**4-fff**2) END DO END DO c kinetic energy (in x) DO i=1,Nx DO j=1,Ny+1 FreeEnergy = FreeEnergy * +0.5d0*rT1*iswlx(i,j) * *(rAx2*CABS(Ux(i,j)*F(i+1,j)-F(i,j))**2) END DO END DO c kinetic energy (in y) DO i=1,Nx+1 DO j=1,Ny FreeEnergy = FreeEnergy * +0.5d0*rT1*iswly(i,j) * *(rAy2*CABS(Uy(i,j)*F(i,j+1)-F(i,j))**2) END DO END DO c field energy DO j=1,Ny DO i=1,Nx IF (BULK(i,j)) THEN rotA=-rAxAy*ATAN2(AIMAG(bloop(i,j)), * REAL(bloop(i,j))) FreeEnergy=FreeEnergy * +k2rT12*rotA*(rotA-2*He) ELSE rotA=Hi FreeEnergy=FreeEnergy+k2rT12*rotA*(rotA-2*He) END IF Bz=Bz+rotA END DO END DO c Bz is the mean field Bz=Bz/Nx/Ny FreeEnergy=FreeEnergy*AxAy PATHx1=0 PATHx2=0 PATHy1=0 PATHy2=0 MinF=10 DO i=1+d,Nx-d j=1+d PATHx1=PATHx1+ATAN2(AIMAG(CONJG(F(i,j))*F(i+1,j)) * ,REAL(CONJG(F(i,j))*F(i+1,j))) IF (MinF .GT. CABS(F(i,j))) THEN MinF=CABS(F(i,j)) END IF j=Ny+1-d PATHx2=PATHx2+ATAN2(AIMAG(CONJG(F(i,j))*F(i+1,j)) * ,REAL(CONJG(F(i,j))*F(i+1,j))) IF (MinF .GT. CABS(F(i,j))) THEN MinF=CABS(F(i,j)) END IF END DO DO j=1+d,Ny-d i=1+d PATHy1=PATHy1+ATAN2(AIMAG(CONJG(F(i,j))*F(i,j+1)) * ,REAL(CONJG(F(i,j))*F(i,j+1))) IF (MinF .GT. CABS(F(i,j))) THEN MinF=CABS(F(i,j)) END IF i=Nx+1-d PATHy2=PATHy2+ATAN2(AIMAG(CONJG(F(i,j))*F(i,j+1)) * ,REAL(CONJG(F(i,j))*F(i,j+1))) IF (MinF .GT. CABS(F(i,j))) THEN MinF=CABS(F(i,j)) END IF END DO Flux=(PATHx1-PATHx2-PATHy1+PATHy2)*0.1591549430919 PATHx1=0 PATHx2=0 PATHy1=0 PATHy2=0 MinF1=10 DO i=W,E-1 j=S PATHx1=PATHx1+ATAN2(AIMAG(CONJG(F(i,j))*F(i+1,j)) * ,REAL(CONJG(F(i,j))*F(i+1,j))) IF (MinF1 .GT. CABS(F(i,j))) THEN MinF1=CABS(F(i,j)) END IF j=N PATHx2=PATHx2+ATAN2(AIMAG(CONJG(F(i,j))*F(i+1,j)) * ,REAL(CONJG(F(i,j))*F(i+1,j))) IF (MinF1 .GT. CABS(F(i,j))) THEN MinF1=CABS(F(i,j)) END IF END DO DO j=S,N-1 i=W PATHy1=PATHy1+ATAN2(AIMAG(CONJG(F(i,j))*F(i,j+1)) * ,REAL(CONJG(F(i,j))*F(i,j+1))) IF (MinF1 .GT. CABS(F(i,j))) THEN MinF1=CABS(F(i,j)) END IF i=E PATHy2=PATHy2+ATAN2(AIMAG(CONJG(F(i,j))*F(i,j+1)) * ,REAL(CONJG(F(i,j))*F(i,j+1))) IF (MinF1 .GT. CABS(F(i,j))) THEN MinF1=CABS(F(i,j)) END IF END DO Flux1=(PATHx1-PATHx2-PATHy1+PATHy2)*0.1591549430919 Magnetization = (Bz - He)*0.079577472 WRITE (6,*) 'Step:',Step,' -----------------------------' WRITE (6,*) 'Energy: ',FreeEnergy WRITE (6,*) 'He: ',He ,' Hi: ',Hi WRITE (6,*) 'Bz: ',Bz , * ' Magnetization:',Magnetization WRITE (6,*) 'MinF: ',MinF,' NumVort: ',Flux OPEN (13,ACCESS='APPEND',FILE='s.dat',FORM='FORMATTED') WRITE (13,321) Step,He,Hi,Bz,Magnetization, * FreeEnergy,Flux,Flux1,MinF,MinF1 321 format(i8,' ',e14.7,' ',e14.7,' ',e14.7,' ',e14.7, * ' ',e14.7,' ',e14.7,' ',e14.7,' ',e14.7,' ',e14.7) CLOSE(13) IF (MOD(Step,OutSteps) .LE. 1) THEN OPEN (15,ACCESS='APPEND',FILE='c.dat',FORM='FORMATTED') WRITE (15,*) He,Bz,Flux,Hi CLOSE (15) iaux = Step + 70000000 write (Label,51) iaux 51 FORMAT (i8,'.gra') open (11,file=Label) write (11,*) ' Step = ', Step,' HH = ',He, * ' Mesh of ',Nx,' x ',Ny write (11,52) ((cabs(f(i,j)),i=1,Nx+1),j=1,Ny+1) 52 format (5(f8.4)) close (11) END IF END IF c ****************************************************************** c *************************** Salida ******************************* IF (MOD(Step,OutSteps) .LE. 1) THEN iHe = iHe + 1 WRITE (Label,50) iHe 50 FORMAT (I7,'.ree') WRITE (6,*) 'Creating a restart file: ',Label OPEN (10,FILE=Label,FORM='FORMATTED') WRITE (10,*) Hi , He WRITE (10,*) ((F(i,j),i=1,Nx+1),j=1,Ny+1) WRITE (10,*) ((Ux(i,j),i=1,Nx),j=1,Ny+1) WRITE (10,*) ((Uy(i,j),i=1,Nx+1),j=1,Ny) CLOSE (10, STATUS='KEEP') END IF c ****************************************************************** IF (MOD(Step,HeSteps) .LE. 1) THEN He=Ho+dHe*Step AxAyHe=Ax*Ay*He EiAxAyHe=(1-0.5*AxAyHe**2+AxAyHe**4/24.)+ * (0,1)*(-AxAyHe+AxAyHe**3/6.-AxAyHe**5/120.) END IF c ****************************************************************** END DO c ****************************************************************** c ****************************************************************** CLOSE (9, STATUS='KEEP') CLOSE (12, STATUS='KEEP') CLOSE (14, STATUS='KEEP') CLOSE (16, STATUS='KEEP') CLOSE (17, STATUS='KEEP') c **************************** Bye ********************************* WRITE (6,*) 'END' END c ****************************************************************** c ****************************************************************** c ********************* Gaussian noise ***************************** COMPLEX FUNCTION GaussDev (idum) REAL MyRAND 1 v1=2.*MyRAND(idum)-1. v2=2.*MyrAND(idum)-1. r=v1**2+v2**2 IF(r .ge. 1 .or. r .eq. 0)GO TO 1 fac=sqrt(-2.*log(r)/r) GAUSSDEV=CMPLX(v1*fac,v2*fac) RETURN END REAL FUNCTION MyRAND(idum) DIMENSION r(97) PARAMETER (m1=259200,ia1=7141,ic1=54773,rm1=1./m1) PARAMETER (m2=134456,ia2=8121,ic2=28411,rm2=1./m2) PARAMETER (m3=243000,ia3=4561,ic3=51349) DATA iff /0/ IF (idum .lt. 0 .or. iff .eq. 0) THEN iff=1 ix1=MOD(ic1-idum,m1) ix1=MOD(ia1*ix1+ic1,m1) ix2=MOD(ix1,m2) ix1=MOD(ia1*ix1+ic1,m1) ix3=MOD(ix1,m3) DO j=1,97 ix1=MOD(ia1*ix1+ic1,m1) ix2=MOD(ia2*ix2+ic2,m2) r(j)=(FLOAT(ix1)+FLOAT(ix2)*rm2)*rm1 END DO idum=1 END IF ix1=MOD(ia1*ix1+ic1,m1) ix2=MOD(ia2*ix2+ic2,m2) ix3=MOD(ia3*ix3+ic3,m3) j=1+(97*ix3)/m3 IF(j .GT. 97 .OR. j .LT. 1)PAUSE MyRAND=r(j) r(j)=(FLOAT(ix1)+FLOAT(ix2)*rm2)*rm1 RETURN END