Silverfrost Forums

Welcome to our forums

help to fix this code

Back to Forum: General To End of Page
pashupati
Posts: 8
Back to Top
21 Aug 2012 1:45 #10661

I can only get first data point after that it shows error.

Please help

  PROGRAM TDGL
  
  
  IMPLICIT NONE

!c **************** Initialization of variables ********************* INTEGER, 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
  
  integer iaux

!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.5d0fff4-fff2) END DO END DO !c kinetic energy (in x) DO i=1,Nx DO j=1,Ny+1 FreeEnergy = FreeEnergy & +0.5d0rT1iswlx(i,j) & (rAx2CABS(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.5d0rT1iswly(i,j) & (rAy2CABS(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=-rAxAyATAN2(AIMAG(bloop(i,j)), & REAL(bloop(i,j))) FreeEnergy=FreeEnergy & +k2rT12rotA(rotA-2He) ELSE rotA=Hi FreeEnergy=FreeEnergy+k2rT12rotA*(rotA-2He) END IF Bz=Bz+rotA END DO END DO !c Bz is the mean field Bz=Bz/Nx/Ny FreeEnergy=FreeEnergyAxAy

          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=v12+v22 IF(r .ge. 1 .or. r .eq. 0)GO TO 1 fac=sqrt(-2.log(r)/r) GAUSSDEV=CMPLX(v1fac,v2fac) 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

tdgl.ini

'Number of cells along x::' 64 'Number of cells along y::' 64 'Mesh spacing along x::' 0.5 'Mesh spacing along y::' 0.5 'Time step::' 0.015 'Kappa::' 2.0 'Eta::' 1. 'Temperature::' 0.5 'Initial applied field::' 0.0 'Field increment::' 1.d-06 'Noise constant::' 1.000000E-05 'Total number of steps::' 250000 'Steps between outputs::' 50000 'Steps between stdouts::' 500 'Steps between increments in He::' 2 'Distance from path to boundary::' 1 'N boundary of hole::' 49 'S boundary of hole::' 17 'E boundary of hole::' 49 'W boundary of hole::' 17

[/code]

jjgermis
Posts: 400 Nürnberg, Germany
Back to Top
22 Aug 2012 8:52 #10665

Thanks for posting some code - this is alway helps to get down to the problem. However, since the number of code lines posted in the forum is limited, you should try to reduce the code so that one can at least reproduce it.

I always like cut, paste and compile examples 😃

Give it try - and you will get the help you are looking for.

LitusSaxonicum
Posts: 2284 Yateley, Hants, UK
Back to Top
22 Aug 2012 10:36 #10667

If the code is long, you could always upload it to one of the (legal) file sharing websites and provide a link to it. I often use DropBox (www.dropbox.com) to share files that are too long even for e mail attachments or are of types like .EXE that many e mail systems won't accept.

Eddie

jjgermis
Posts: 400 Nürnberg, Germany
Back to Top
22 Aug 2012 10:45 #10668

Hi Eddie,

thanks for the DropBox tipp. It seems like one need to install some software before one can use it. We for example are not allowed to install software that is not approved - difficult situation. But on the otherhand our IT always try to avoid the worst case scenario.

pashupati
Posts: 8
Back to Top
23 Aug 2012 2:04 #10677

Anyone willing to help me to take a look at the code? Please help

my email :dklpashupati at yahoo.com

jjgermis
Posts: 400 Nürnberg, Germany
Back to Top
23 Aug 2012 2:16 #10678

We are willing to help if you post a working example of your code. The part of the program given above is only variable declaration. It is not fair to expect from anybody to figure out what is happening with such little information.

The only 'real' code is the line where a file is opened. So for a start you must at least specify: 1.) are you able to read the file? 2.) occur the error before/after any calls to subroutines or functions 3.) what sort of files are you reading.

You need to reduce you problem. Once you have done that we can help you. Give it a try 😉

pashupati
Posts: 8
Back to Top
23 Aug 2012 2:35 #10679

Sorry about the code, I didn't realize that only a part of the code is posted 😦 . I am reading .ini files with the input parameters. After I start 'run', it reads the parameter and generate the first data point. After that it gives 'floating point overflow'. floating point co-processor fault at address 0040a7e5 in the line 917'. I am running this in window 7 32 bit laptop.

I used a gfortan in mac computer the same code works fine. I am not planning to buy Mac. 😢

jjgermis
Posts: 400 Nürnberg, Germany
Back to Top
23 Aug 2012 2:49 #10680

You have to post something like this:

open(8,file=tdgl.ini')

C
C your code here
C

close(8)

and if possible what the file looks like. Without this information you will sit on your problem - not what you want 😉

34455
jhgfd
jhgfd
12345 6788
pashupati
Posts: 8
Back to Top
23 Aug 2012 5:40 #10681 
   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)
pashupati
Posts: 8
Back to Top
23 Aug 2012 5:41 #10682 
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
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)
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)
LitusSaxonicum
Posts: 2284 Yateley, Hants, UK
Back to Top
23 Aug 2012 5:42 #10683

Jaques,

Drop Box is only one of several legal file sharing sites. You only need to instal software to post to it, as reading from it is a simple internet download. You get 2 Gb free, but if you buy a licence, you can have much more space, and then your IT section could add it to the list of accepted applications (or whatever system you chose would probably work the same). In my experience, IT organisations distrust free software, but they are happy to pay for stuff you can get free - and then they won't pay for the updates!

I found it useful for sharing EXE files.

Eddie

pashupati
Posts: 8
Back to Top
23 Aug 2012 5:48 #10684

Code is very long May I email to someone willing to help me?

email: dklpashupati at yahoo.com

Thank you !

jjgermis
Posts: 400 Nürnberg, Germany
Back to Top
24 Aug 2012 6:51 #10685

Seeing some code is a help - we are getting there 😃

However, where does the error occur. Can you isolate it: 1.) are you able to read the file? 2.) occur the error before/after any calls to subroutines or functions 3.) what sort of files are you reading.

pashupati
Posts: 8
Back to Top
24 Aug 2012 6:25 #10689

I can read the input file and it gives the first data point (valid point) and it terminate due floating point error

DanRRight
Posts: 2877 South Pole, Antarctica
Back to Top
25 Aug 2012 5:34 #10695

In the DOS prompt window (better use for that the Total Commander - to me it's #1 program on all computers and the best and the most useful app of all times) compile it like this

FTN95 yourfile.for /checkmate /link

Yourfile.for is the name of your file of course. If this is free form source add /free to this line

Then run the debugger

SDBG yourfile.exe

You will find all your errors in seconds. If not ask here again what's specifically wrong. Post here how many errors Apple compiler missed

LitusSaxonicum
Posts: 2284 Yateley, Hants, UK
Back to Top
30 Aug 2012 6:56 #10703

Since only the start of the program is given, and no example data, this is a case for clairvoyance or mind reading. While I make no claims to expertise in either, I have a feeling I know the answer.

Is it possible (probable?) that the program expects one data value per line because of the form of the READ statements, but the datafile doesn't have that layout? Or that LABEL (a CHARACTER*36 entity) isn't in quotation marks in one or more lines?

Oddly enough, my crystal ball appears to have developed a pattern akin to the old 80-column Fortran coding sheets .....

Mystic Ed

Please login to reply.