Ev2.f

Go to the documentation of this file.

C************************************************************
C
C       NVAR=1 - SQUARED FORM OF SURFACES IS USED:
C       A1*X**2+A2*X*Y+A3*Y**2+A4*X+A5*Y+A6=0
C
C       NVAR=2 - BILINEAR FORM OF SURFACES IS USED
C
C       NVAR = 11 - TEST FOR NVAR=1
C--------------------------------------------------------
        SUBROUTINE fun(X,F)
C
        IMPLICIT REAL*8(a-h,o-z)
C
        common/cpcoor/ nvar,ifailf,rcp,zcp,uu(5),df
        common/yscr/ y
C
        sqrt(xxx) = dsqrt(xxx)
C
        IF( nvar .EQ. 1 ) THEN
C
            a1 = 0.5d0*uu(3)
            a2 =       uu(4)
            a3 = 0.5d0*uu(5)
            a4 =       uu(1)
            a5 =       uu(2)
            a6 = df
C
            alf = 1.d0
            ifailf = 0
C
            b = (a2*x+a5) / (2.d0*a3)
            det = b**2 - (a1*x**2+a4*x+a6)/a3
            IF( det .LT. 0 ) THEN
C                WRITE(6,*) ' PROGRAM TERMINATED: DET= ',DET
C                STOP
                 det = - det
                 alf = 1.d+3
                 ifailf = 1
            ENDIF
            det = sqrt( det ) * alf
            y  =  -b + det
            y2 =  -b - det
C
            f  = (x-rcp)**2 + (y -zcp)**2
            f2 = (x-rcp)**2 + (y2-zcp)**2
C
            IF (f2 .LT. f ) THEN
              f = f2
              y = y2
            ENDIF
C
        ENDIF
C
C
        IF( nvar .EQ. 11 ) THEN
C..TEST
          a1 = 1.d0
          a2 = 0.d0
          a3 = 1.d0
          a4 = 0.d0
          a5 = 0.d0
          a6 = -1.d0
C......
            b = (a2*x+a5) / (2.d0*a3)
            det = sqrt( b**2 - (a1*x**2+a4*x+a6)/a3 )
            y  =  -b + det
            y2 =  -b - det
C
            f  = (x-rcp)**2 + (y -zcp)**2
            f2 = (x-rcp)**2 + (y2-zcp)**2
C
            IF (f2 .LT. f )  f = f2
C
        ENDIF
C
        RETURN
        END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C***************************************************************
C Definition of (IC,JC)-cell, wich contain given (RC,ZC)-point
C---------------------------------------------------------------
C
          SUBROUTINE ijmesh( NR,R,NZ,Z, RC,ZC, IC,JC)
C
         IMPLICIT REAL*8(a-h,o-z)
C
         dimension r(1),z(1)
C
          nout  = 17
          nter  = 6
C
         IF( (rc.LT.r(1)) .OR. (r(nr).LT.rc) ) THEN
           ! WRITE(NTER,*) 'PROGRAM WAS INTERRUPTED IN SUBR. "IJMESH" '
           ! WRITE(NTER,*) 'POINT "RC" IS OUT OF "BOX" : RC =', RC
            !WRITE(NOUT,*) 'PROGRAM WAS INTERRUPTED IN SUBR. "IJMESH" '
            !WRITE(NOUT,*) 'POINT "RC" IS OUT OF "BOX" : RC =', RC
            stop
         END IF
         IF( (zc.LT.z(1)) .OR. (z(nz).LT.zc) ) THEN
           ! WRITE(NTER,*) 'PROGRAM WAS INTERRUPTED IN SUBR. "IJMESH" '
           ! WRITE(NTER,*) 'POINT "ZC" IS OUT OF "BOX" : ZC =', ZC
            !WRITE(NOUT,*) 'PROGRAM WAS INTERRUPTED IN SUBR. "IJMESH" '
            !WRITE(NOUT,*) 'POINT "ZC" IS OUT OF "BOX" : ZC =', ZC
            stop
         END IF
C
         DO 1 i=1,nr
            IF( r(i) .GT. rc ) THEN
                ic = i - 1
                go to 2
            ENDIF
 1       CONTINUE
 2       CONTINUE
C
         DO 3 j=1,nz
            IF( z(j) .GT. zc ) THEN
                jc = j - 1
                go to 4
            ENDIF
 3       CONTINUE
 4       CONTINUE
         RETURN
         END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C***********************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C      SUM OF ANGLES FOR POINT (UC,VC) AND
C          CLOSED  CURVE U,V(N)
C      ASSUMING    U,V(N) .NE. U,V(1) !!
C----------------------------------------------------------
C
         REAL*8 FUNCTION sumang (N,U,V,UC,VC)
C
         include'double.inc'
C
         dimension u(1:n), v(1:n)
C
         sumang=0.
         reg=1.d-10
C
           DO 1 i =1,n
           du1=u(i)-uc
           dv1=v(i)-vc
         IF(i.NE.n) THEN
           du2=u(i+1)-uc
           dv2=v(i+1)-vc
         ELSE
           du2=u(1)-uc
           dv2=v(1)-vc
         ENDIF
         ang=dacos((du1*du2+dv1*dv2)/(dsqrt((du1*du1+dv1*dv1)*
     &                                (du2*du2+dv2*dv2))+reg))
         IF(du1*dv2-du2*dv1.GT.0.) THEN
          sumang=sumang+ang
         ELSE
          sumang=sumang-ang
         ENDIF
 1        CONTINUE
C
         RETURN
         END
C**********************************************************
C     MIN. DISTANCE FOR POINT (UC,VC) AND
C              CURVE U,V(N)
C----------------------------------------------------------
         REAL*8 FUNCTION distan(N,U,V,UC,VC)
C
         include'double.inc'
C
         dimension u(1:n), v(1:n)
C
         sqrt(xxx) = dsqrt(xxx)
C
         distan = 100.d0
C
         DO 1 i =1,n
            ttt = (u(i)-uc)**2 + (v(i)-vc)**2
            ttt = sqrt(ttt)
            IF(ttt.LT.distan) distan = ttt
 1       CONTINUE
C
         RETURN
         END
C
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C   POLOIDAL mag. field companents Bp_r, Bp_z  and
C   POLOIDAL mag. field projection on direction FIPRO
C   COMPUTE FOR GIVEN "PF_PROBE" POINTS:
C-----------------------------------------------------------------------
C
        SUBROUTINE bp_prob( NOUT,   NTER,   KEYPRI,
     *                      npro,   rpro,   zpro,   fipro,
     *                      brpro,  bzpro,  bpcom )
C
        include 'double.inc'
C
        dimension rpro(*),  zpro(*),  fipro(*),
     *            brpro(*), bzpro(*), bpcom(*)
C-----------------------------------------------------------------------
C
      IF( npro.NE.0 )  THEN
         do 1455 i=1,npro
            call probe( rpro(i), zpro(i), fipro(i),
     *                  brpr, bzpr, bpnum )
            brpro(i) = brpr
            bzpro(i) = bzpr
            bpcom(i) = bpnum
 1455    continue
C
C        WRITE(NOUT, * )  '   '
C        WRITE(NOUT, * ) '--------------------------------------------'
C        WRITE(NOUT, * ) 'Poloidal mag. field probe projection        '
C        WRITE(NOUT, * ) '  BPCOM(i), i=1,...,NPRO = ', NPRO
C        WRITE(NOUT, * ) '--------------------------------------------'
C        WRITE(NOUT,101) (BPCOM(I), I=1,NPRO)
C        WRITE(NOUT, * ) '   '
C
      END IF
C                 
  101  FORMAT(2x,5e14.7)
C
         RETURN
         END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C  Compute POLOIDAL mag. field companents Bp_r, Bp_z
C  and POLOIDAL mag. field projection on direction FPR
C-----------------------------------------------------------------------
C
        SUBROUTINE probe( RPR, ZPR, FPR, BRPR, BZPR, BPNUM )
C
        include 'double.inc'
        include 'param.inc'
        include 'comblc.inc'
C
        sin(xxxx) = dsin(xxxx)
        cos(xxxx) = dcos(xxxx)
C-----------------------------------------------------------------------
C
          CALL ijmesh( ni,r,nj,z, rpr,zpr, ic,jc )
C
          drdz  = ( r(ic+1) - r(ic) )*( z(jc+1) - z(jc) )
C
          dpsdr = ( (u(ic+1,jc  ) - u(ic,jc  )) * (z(jc+1) - zpr  ) +
     *              (u(ic+1,jc+1) - u(ic,jc+1)) * (zpr     - z(jc))  )
     *          / drdz
          dpsdz = ( (u(ic  ,jc+1) - u(ic  ,jc)) * (r(ic+1) - rpr  ) +
     *              (u(ic+1,jc+1) - u(ic+1,jc)) * (rpr     - r(ic))  )
     *          / drdz
          brpr  = - dpsdz / rpr
          bzpr  =   dpsdr / rpr
C
          bpnum = brpr*cos(fpr) + bzpr*sin(fpr)
C
          bbb   = 1.d0/(2.d0*pi)
C
          brpr  = brpr  
          bzpr  = bzpr  
          bpnum = bpnum 
C
        RETURN
        END
C
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C    Compute full poloidal flux in given "loop" points
C-----------------------------------------------------------------------
C
        SUBROUTINE fl_loo( NOUT,   NTER,   KEYPRI,
     *                     nloo,   rloo,   zloo,
     *                     psloo )
C
        include 'double.inc'
C
        dimension rloo(1), zloo(1), psloo(1)
C-----------------------------------------------------------------------
C
      IF( nloo.NE.0 )  THEN
         do i=1,nloo
C
            call floop( rloo(i), zloo(i), psnum )
C
            if( i.eq.1 ) then
                    psloo(i) = psnum
            else
              psloo(i) = psnum - psloo(1)
            end if
C
         end do
C
C        WRITE(NOUT, * ) '   '
C        WRITE(NOUT, * ) '--------------------------------------------'
C        WRITE(NOUT, * ) 'Poloidal mag. flux loop values              '
C        WRITE(NOUT, * ) '  PSLOO(i), i=1,...,NLOO = ', NLOO
C        WRITE(NOUT, * ) '--------------------------------------------'
C        WRITE(NOUT,101) (psloo(i), i=1,nloo)
C        WRITE(NOUT, * ) '   '
C
      END IF
C                 
  101  FORMAT(2x,5e14.7)
C
         RETURN
         END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C    Compute full poloidal flux in given "loop" point
C-----------------------------------------------------------------------
        SUBROUTINE floop( RPR, ZPR, PSPR )
C
        include 'double.inc'
        include 'param.inc'
        include 'comblc.inc'
C
C-----------------------------------------------------------------------
C
          CALL ijmesh( ni,r,nj,z, rpr,zpr, ic,jc )
C
          drdz  = ( r(ic+1) - r(ic) )*( z(jc+1) - z(jc) )
C
          pspr  = ( u(ic  ,jc  ) * (r(ic+1) - rpr)*(z(jc+1) - zpr) +
     *              u(ic+1,jc  ) * ( rpr - r(ic) )*(z(jc+1) - zpr) +
     *              u(ic  ,jc+1) * (r(ic+1) - rpr)*( zpr - z(jc) ) +
     *              u(ic+1,jc+1) * ( rpr - r(ic) )*( zpr - z(jc) )   )
     *          / drdz
C
          bbb   = 1.d0/(2.d0*pi)
C
          pspr   = pspr / bbb
C
        RETURN
        END
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!