sstepon.f

Go to the documentation of this file.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
       SUBROUTINE  sstepon( KLUCH, k_auto, nstep, dt, time,
     *                     voltpf, d_pf_mat,d_tcam_mat,key_dmf )

         use durs_d_modul       

c---------------------------------------------------------------
       include 'Descript.inc'
c  ----------
c
       include 'double.inc'
c
       include 'prm.inc'
       include 'comevl.inc'
       !INCLUDE 'contro.inc'
       !parameter(nstep_p=3000)
c-----------------------------------------------------------------------
       include 'comtim.inc'

         common/comeqg/  ncequi
       common /nostep/ kstep
       common /comhel/ helinp, helout
       
!       common /comsta/ epsro,platok,
!     *                 eqdfn,betplx,tokf,psax,b0,r0,
!     *                 alf0,alf1,alf2,bet0,bet1,bet2,rax,zax,
!     *                 nurs,igdf,n_psi,n_tht,i_bsh,keyctr,i_betp,i_eqdsk
     
       common /comsta/ platok,eqdfn,i_bsh
       
       common /com234/ betpol,tokout,psiout
       common /com_cam/ camtok
         common /c_kpr/ kpr
       common /key_for_pet/ key_fixbon
       common/com_flag/kastr
c-----------------------------------------------------------------------
c
       dimension  rc(nclim),  zc(nclim),  pc(nclim), psip(nclim),
     *            vc(nclim),  hc(nclim)
c
       dimension  rc1(nclim), zc1(nclim), rc2(nclim), zc2(nclim),
     *            rc3(nclim), zc3(nclim), rc4(nclim), zc4(nclim)
c
       INTEGER    necon(nilim), ntype(nclim)
       dimension  wecon(nilim)
c
       common/comst0/ res(njlim),  volk(njlim),  volkp1(njlim)
c
       common/comst1/ pjk(njlim),pjkp1(njlim),pjkp(njlim),pjkd(njlim)
       include 'parloo.inc'
       common/comloo/ rloop(nloopp),zloop(nloopp),
     &          rprob(nprobp),zprob(nprobp),fiprob(nprobp),nloop,nprob
       common/comst2/ psk(njlim), pskp1(njlim),pskp(njlim),pskm1(njlim)

       !DIMENSION  cp_com_0(n_cp_m)
c***********************************************************************

       dimension  psiplb(nstep_p),psiexb(nstep_p),psimag(nstep_p),
     *            flu_tor(nstep_p)  

c***********************************************************************

       real*8 voltpf(*)
       real*8 d_pf_mat(*)
       real*8 d_tcam_mat(*)

c***********************************************************************


         real*4       a_print(200)
         character*30 apr
       character*40 eqdfn


          ! call put_tim(dt,time)

         kstep = nstep
         timev = time
         tstep = dt


!(((((((((((((((((((((((
!       i_diag=1
!)))))))))))))))))))))))

         pi  = 3.14159265359d0
         bbb  = 1.d0 / (2.d0*pi)

       sigm   = 1.d0  !0.50d0

       tstepr = tstep
       tstart = 0.00d0

       tstop  = 20.0d0
       kstop  = 1000

       enels  = 5.0d-8
       knels  = 500

       nfrpr1 = 5
       nfrwr1 = 0

          nout  = 17
          nter  = 6
          ninfw = 7
          ninev = 10
          ngra1 = 14
          ngra2 = 15
c


       IF( kluch .NE. 0 ) go to 1111

 1149     format(a40)
c
      ! OPEN( NOUT, FILE='EVOL0.PRT', FORM='FORMATTED' )
       !write(nout,*) '***********************************************'
      ! write(nout,*) '!               OUTPUT FILE                   !'
      ! write(nout,*) '!                   OF                        !'
      ! write(nout,*) '!  the PET code (Plasma Evolution in Tokamak) !'
      ! write(nout,*) '! ------------------------------------------- !'
      ! write(nout,*) '!  free boundary tokamak equilibrium problem  !'
      ! write(nout,*) '!                    and                      !'
      ! write(nout,*) '!             circuit equations               !'
      ! write(nout,*) '***********************************************'
c
       numwr  = 0
c--------------------------------------------------------------------
c
c======iter, sof
c
c=============================================

       kstepr = kstep
       keypri = 1

c-----------------------------------------------------------------
c --- input of "BASIC" equlibrium parameters from file "durs.dat"
c
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

         if(k_auto.eq.0) go to 2005 

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         if(kastr.eq.1) go to 1150 

        write(fname,'(a,a)') path(1:kname),'durs.dat'
        open(1,file=fname)
          !open(1,file='durs.dat')

              read(1,*) n_tht
              read(1,*) n_psi
              read(1,*) igdf
              read(1,*) epsro

              read(1,*) nurs
              read(1,*) keyctr

              read(1,*)    i_eqdsk
              read(1,1149) eqdfn
 
              read(1,*) i_betp
              read(1,*) betplx
              read(1,*) tokf
              read(1,*) psax
              read(1,*) b0,r0
                             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              read(1,*) alf0 !!
              read(1,*) alf1 !! P'=alf0*(1-(1-psi)**alf1)**alf2
              read(1,*) alf2 !!
                             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              read(1,*) bet0 !!
              read(1,*) bet1 !! FF'=bet0*(1-(1-psi)**bet1)**bet2
              read(1,*) bet2 !!
                             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              read(1,*) alw0 !!
              read(1,*) alw1 !! (ro*w**2)'=alw0*(1-(1-psi)**alw1)**alw2
              read(1,*) alw2 !!
                             !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              read(1,*) rax
              read(1,*) zax

          close(1)

         if(i_eqdsk.eq.1) then

               call tab_efit(tokf,psax,eqdfn,rax,zax,b0,r0)

               nurs   = -3999
               i_betp = 0

         endif

        write(fname,'(a,a)') path(1:kname),'durs_d.dat'
        open(1,file=fname,form='formatted')
          !open(1,file='durs_d.dat')
              write(1,*) n_tht,n_psi,igdf,nurs,keyctr,i_eqdsk,i_betp
              write(1,*) epsro,betplx,tokf,psax,b0,r0,rax,zax
              write(1,*) alf0,alf1,alf2,bet0,bet1,bet2
          close(1)

 1150   continue


c***********************************************************************
c--- input of positions of "PF_PROBE" points:
c
        CALL propnt( nout, nter, ninfw, ngra1,
     *               npro, rprob, zprob, fiprob )
c***********************************************************************
c--- input of positions of "FL_LOOP" points:
c
        CALL loopnt( nout, nter, ninfw, ngra1,
     *               nloop, rloop, zloop )
c***********************************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c***********************************************************************
c--- input parameters of pfc system and passiv conductors
c
       CALL conduc( nc, ncequi, ncpfc, nfw, nbp, nvv,
     *              rc,zc, pc, vc,hc, ntype,
     *              rc1,zc1,  rc2,zc2,  rc3,zc3,  rc4,zc4,
     *              res, volk, volkp1,
     *              necon, wecon,
     *              nout, nter, ninfw, ngra1 )



c
c***********************************************************************
c--- definition induct. and  selfinduct. matrix
c--- for "EDDY" conductors:  "PPIND" from COMMON /ppidps/
c                            *******
c
          CALL l_matr( nout,  nter, nc, ncpfc,
     *               ntype, rc,   zc, vc, hc,
     *               necon,wecon )
c
c-----------------------------------------------------------------------
c --- initial condition(currents) for circuit equations
c

      DO 17 l=1,ncequi
          IF( l.LE.nequi ) THEN
           pjk(l)  = pfceqw(l)
           pjkp(l) = pjk(l)
        ELSE
           pjk(l)  = pc(ncpfc+l-nequi)
           pjkp(l) = pjk(l)
        END IF
   17 CONTINUE



        write(fname,'(a,a)') path(1:kname),'currents.wr'
        open(1,file=fname,form='formatted')
       !open(1,file='currents.wr',form='formatted')
        write(1,*) nequi,ncequi
        write(1,*)(pjk(j),j=1,ncequi)
       close(1)

        write(fname,'(a,a)') path(1:kname),'res_mat.wr'
        open(1,file=fname,form='formatted')
       !open(1,file='res_mat.wr')
        write(1,*) (res(j),j=1,ncequi)
       close(1)

c
c.....  toroidal currents of passive conductor structures 
c
!        CALL PASCUR( NOUT, NTER, NEQUI, NFW, NBP, NVV,
!     *               PJK, fwcurr, bpcurr, vvcurr )
c
c********************************************************************

!         WRITE(NOUT,*) '********************************************'
!         WRITE(NOUT,*) ' NUMBER TIME-LEVEL  KSTEP =', KSTEP
!         WRITE(NOUT,*) '********************************************'

cccc      include 'printb.inc'

          helout = helinp
c

!!!!!  WRITE(NTER,*) 'START  OF BASIC FREE BOUNDARY EQUILIBRIUM'

         call wrcoil(nc,ncpfc,rc,zc,pc,necon,wecon)

         ngrid=1

         call auto(rc,zc,pc,nc,nstep,ngrid,
     *                 rc1,zc1, rc2,zc2,
     *                 rc3,zc3, rc4,zc4,
     *                 ntype, necon, wecon )


!(((((((((((((((((((((((((((((((((((((((((
!       if(i_diag.eq.1)then
!       call  diag_params()
!       call diag_MATR(NC, NCPFC, nequi,
!     *                   RC,ZC,NECON,WECON,
!     * NLOO, RLOO, ZLOO ,
!     * NPRO, RPRO, ZPRO, FIPRO
!     *  )
!))))))))))))))))))))))))))))))))))))))))))


!((((((((((((((((((((((((((
!       call  write_diag()
!))))))))))))))))))))))))))


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 2005   continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         if(kastr.eq.0 .AnD. i_eqdsk.eq.0) then 
        write(fname,'(a,a)') path(1:kname),'inpol.dat'
        open(1,file=fname,form='formatted')
          !open(1,file='inpol.dat')
             read(1,*) i_bsh
          close(1)    
         else
             i_bsh=1
         endif

        if(kastr.eq.0 ) then 
        write(fname,'(a,a)') path(1:kname),'durs_d.dat'
        open(1,file=fname,form='formatted')
        !open(1,file='durs_d.dat')
            read(1,*) n_tht,n_psi,igdf,nurs,keyctr,i_eqdsk,i_betp
            read(1,*) epsro,betplx,tokf,psax,b0,r0,rax,zax
            read(1,*) alf0,alf1,alf2,bet0,bet1,bet2
        close(1)
        endif


!!:::::::::::::::::__________________________________:::::::::::::::

        call noauto

        call rd_ppind

        call rd_prob( nprob, rprob, zprob,  fiprob )

        call rd_loop( nloop, rloop, zloop )


        write(fname,'(a,a)') path(1:kname),'currents.wr'
        open(1,file=fname,form='formatted')
       !open(1,file='currents.wr',form='formatted')
        read(1,*) nequi,ncequi
        read(1,*)(pjk(j),j=1,ncequi)
       close(1)

        write(fname,'(a,a)') path(1:kname),'res_mat.wr'
        open(1,file=fname,form='formatted')
       !open(1,file='res_mat.wr')
        read(1,*) (res(j),j=1,ncequi)
       close(1)

       !fr_cam=0.d0
       !do j=nequi+1,ncequi
       ! fr_cam=fr_cam+1.d0/res(j)
       !enddo
       ! fr_cam=1.d0/fr_cam

       do l=1,ncequi
           pjkp(l) = pjk(l)
       enddo


!(((((((((((((((((((((((((((((((((((((((((((
!       if(i_diag.eq.1)then 
!               call  diag_params()
!               call  read_diag()
!               call curr_to_dina(ncequi,pjk,tokf)
!       end if
!)))))))))))))))))))))))))))))))))))))))))))

cccc   voltage for zero time-level
c
       DO i=1,nequi
          volk(i) = voltpf(i)*bbb
       END DO
       DO i=nequi+1,ncequi
          volk(i) = 0.0d0
       END DO



         platok = tokf
         kstep  = 0
           k_step=0
           ngav=keyctr



         call eqb( alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2,
     *             betplx, i_betp,
     *             ngav, k_step, platok, rax,zax, b0,r0, psax, igdf,
     *             n_tht, n_psi, epsro, nurs, i_eqdsk,i_bsh,
     *             psi_bnd,psi0_bnd)



c       call pau()

         ngrid=1


              nursb=nurs
              psicen=psax
              ngav1=0
              ftok=tokf


         call eq_0(  pjk,psip,ncequi,kstep,ngrid,
     *               alf0,alf1,alf2,bet0,bet1,bet2,
     *               betpol,betplx,zli3,
     *                  ngav1,
     *                 ftok, tokout, psiax, psiout, 
     *                 nursb,psi_bnd,alp_b,rax,zax,n_ctrl,b0,r0 )

         call rdexf(ncequi)

         call eq( pjk, psk, ncequi, kstep, ngrid,
     *                 alf0, alf1, alf2, bet0, bet1, bet2,
     *                 betpol, betplx, zli3,
     *                 ngav1,
     *                 ftok, tokout, psiax, psiout, 
     *                 nursb,psi_bnd,alp_b,rax,zax )


!!!!!!!!!!!!!!!!!!!!!!!!!!ttt!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      !WRITE(NOUT,*) '    '
      !WRITE(NOUT,*) 'PRINT AFTER BASIC FREE BOUNDARY EQUILIBRIUM'

          CALL eq_par( z0cen, alp, alpnew, qcen, nctrl, numlim, up,
     *                   rm,    zm,  rx0,    zx0    )
c***********************************************************************

c --- for diagnostic

          rmax0  = rm
          zmax0  = zm
          zcen0  = z0cen
c     ---------------------
          beold  = betpol
          zlold  = zli3
c     ---------------------
          raxpr  = rm
          zaxpr  = zm
          rmaold = rm
          zmaold = zm
          rxppr  = rx0
          zxppr  = zx0
          rxpold = rx0
          zxpold = zx0
c-----------------------------------------------------------------------
c --- definition of input parameters from "basic" equilibrium

          betplx = betpol
          betap0 = betpol
          ftok   = tokout
          psiax  = psiout
          helinp = helout
c ----------------------------------------------------------------------
c  for psi_axis value time scenario(for ngav1=2 or 3)

          time_sta = tstart
          time_fin = tstop

          psax_sta = psiout

c***********************************************************************

          psibou = up
          psidel = psiout - psibou

cccc      include 'printa.inc' ! ---  PRINTING

c***********************************************************************
c--- printing and writing of numbers and positions of limiter points:

        !CALL PRTLIM( NOUT,  NTER, NGRA1 )
c***********************************************************************
c***********************************************************************
c--------------------- writing in the disk files -----------------------
                 ! include 'wrt1.inc'
c-----------------------------------------------------------------------
       !WRITE(NOUT,*) '--------------------------------------------'
       !WRITE(NOUT,*) '      THE  END  OF  FILE:  "EVOL0.PRT"      '
       !WRITE(NOUT,*) '--------------------------------------------'

       !CLOSE(NOUT)

c***********************************************************************
cccc   OPEN (nout, file='EVOL1.PRT', form='FORMATTED')

c --- If we need to compute the "basic" equilibrium only -> kstop=0

       IF( kstop.EQ.0 ) THEN
           go to 222
       END IF
c***********************************************************************
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


1       continue

         i_bsh  = -1
         i_eqdsk  = 0
         kstep  =  0
         keyctr =  0
         e_psi  =  0.d0


         call eqb( alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2,
     *             betplx, i_betp,
     *             keyctr, kstep, platok, rax,zax, b0,r0, psax, igdf,
     *             n_tht, n_psi, epsro, nurs, i_eqdsk,i_bsh,
     *             psi_bnd,psi0_bnd)


         call get_par(psi_bnd)

         call bongri
         call psib_pla(pspl_av)
         call psib_ext(psex_av)

         psi0_bnd=pspl_av+psex_av


         ztok_n  = platok
         zpsim_n = psax


         return

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

cccc       WRITE(nout,*) ' NEW TIME-LEVEL  K = KSTEP =', kstep

 1111   CONTINUE
               !call test_urs

         keyctr = key_dmf
         !keyctr = 500

         call get_par(psi_bnd)
         call bongri
 !2241     continue
         call psib_pla(pspl_av)
         call psib_ext(psex_av)

         psi0_bnd=pspl_av+psex_av


         ztok_n  = platok
         zpsim_n = psax


!  100 KSTEP = KSTEP + 1       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!      TIMEV = TIMEV + TSTEP   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

       ngav1 = 0

c********************************************************************
c  definition of "KEYPRI" - PARAMETER for printing

!       IF((KSTEP/NFRPR1)*NFRPR1 .EQ. KSTEP)  THEN
!           KEYPRI = 1
!       ELSE
           keypri = 0
!       END IF


c====================================================================
c  definition of input equil. parameters: betplx, psiax, ftok, helinp

        IF(ngav1.EQ.0 .OR. ngav1.EQ.2 .OR. ngav1.EQ.4)
     *     betplx = betpol
        IF( ngav1.EQ.0 .OR. ngav1.EQ.1 )  THEN
           psiax  = psiout
           helinp = helout
        END IF
        IF( ngav1.EQ.2 .OR. ngav1.EQ.3 )  THEN
           ftok   = tokout
           helinp = helout
        END IF
        IF( ngav1.EQ.4 .OR. ngav1.EQ.5 )  THEN
           ftok   = tokout
           psiax  = psiout
        END IF
c
c***********************************************************************
c  computing of "VOLKP1"  voltage  for circuit equation 
c  for new time level kstep
c

         DO i=1,nequi
          volkp1(i) = voltpf(i)*bbb
         END DO
         DO i=nequi+1,ncequi
            volkp1(i) = 0.0d0
         END DO
c-----------------------------------------------------------------------
cccc        include 'printb.inc' !Printing
c***********************************************************************
              ! call test_urs

         call bongri
         call psib_pla(pspl_av)
         call psib_ext(psex_av)
         call get_flfi(flx_fi)

         psi0_bnd=pspl_av+psex_av
         psi0_ax=pspl_av+psex_av+psax
           ztok=platok
           zpsim=psax
           i_bsh=-1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

               !call test_urs

         n_dmf=3
        if(keyctr.eq.0) n_dmf=1

      do it_dmf=1,n_dmf !!!   it_dmf <<< mag.field diffusion iter. loop


        if(kpr.eq.1)print *,' it_dfm==',it_dmf

        if(keyctr.ne.0) then
         call eqb( alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2,
     *             betplx, i_betp,
     *             keyctr, kstep, platok, rax,zax, b0,r0, psax, igdf,
     *             n_tht, n_psi, epsro, nurs, i_eqdsk,i_bsh,
     *             psi_bnd,psi0_bnd )
        endif

c       print *,' psax keyctr platok=',psax,keyctr,platok

c       read (*,*)
               !call test_urs

           call wrb
               !call test_urs

           ftok=platok


        !if(key_fixbon.eq.1) exit


       IF( kstep.EQ.1 ) THEN

          ngrid  = 1

          DO l=1,ncequi
           pskp1(l) = psk(l)
          enddo

          nreg = 0
          nles = 0
               !call test_urs

          CALL evslv( nles,   nreg,  tstep,  tstep, sigm,
     *                ncequi, volk,  volkp1, res,
     *                psk,    pskp1, pjk,    pjkp1, pjkp,
     *                nout,   nter,  keypri, ereve )

               !call test_urs
       END IF    !!! ===>  FOR  IF( KSTEP.EQ.1 )
c***********************************************************************
       IF(kstep.NE.1) THEN

c----------------------------------------------------------------
c  the initial approximation for the case of closed-loop evolution
c
           DO l=1,ncequi
             pjkp(l) = pjk(l)
          enddo
c----------------------------------------------------------------
c
           nreg = 0
           nles = 1
           
               !call test_urs

           CALL evslv( nles,   nreg, tstepr, tstep, sigm,
     *                 ncequi, volk, volkp1, res,
     *                 pskm1,  psk,  pjk,    pjkp1, pjkp,
     *                 nout,   nter, keypri, ereve )


       END IF   !!! ===>  FOR  IF( KSTEP.NE.1 )
c***********************************************************************
               !call test_urs
         CALL eq_ax( pjkp1, pskp1, ncequi,   kstep, ngrid,
     *               alf0, alf1, alf2, bet0, bet1, bet2,
     *               betpol,   betplx,   zli3,
     *                   ngav1,
     *               ftok,     tokout, psiax, psiout, 
     *               enels,    erps,
     *               psi_bnd,alp_b,rax,zax,0 )




!(((((((((((((((((((((((((((((((((((((((((((
!       if(i_diag.eq.1)then
!       call psi_diag()
!       call curr_to_dina(ncequi,pjkp1,ftok)
!       call loopflux()
!       call probefield()
!       end if
!)))))))))))))))))))))))))))))))))))))))))))


         CALL eq_par( z0cen, alp, alpnew, qcen, nctrl, numlim, up,
     *                rm,    zm,  rx0,    zx0   )

c***********************************************************************
c ---  start of the iteration loop 
c***********************************************************************
 
             knel  = 0
 1000        knel = knel + 1


          DO l=1,ncequi
             pskp(l) = pskp1(l)
             pjkp(l) = pjkp1(l)
          ENDDO

          nreg = 0
          nles = 2

          CALL evslv( nles,   nreg,  tstep,  tstep, sigm,
     *                ncequi, volk,  volkp1, res,
     *                psk,    pskp1, pjk,    pjkp1, pjkp,
     *                nout,   nter,  keypri, ereve )




          CALL  differ( pjkp, pjkp1, ncequi, errcu1, errcu2,
     *                  curmax, curmin, nout, nter )

c         sgmcur = 1.00d0
          sgmcur = 0.75d0

          DO l=1,ncequi
             pjkp1(l) = sgmcur*pjkp1(l) + (1.0d0 - sgmcur)*pjkp(l)
          END DO
c-----------------------------------------------------------------------

          CALL eq_ax( pjkp1, pskp1, ncequi,   kstep, ngrid,
     *                alf0, alf1, alf2, bet0, bet1, bet2,
     *                betpol,   betplx,   zli3,
     *                    ngav1,
     *                ftok,     tokout, psiax, psiout, 
     *                enels,    erps,
     *                psi_bnd,alp_b,rax,zax,0 )

            write(*,*) 'stepon:EQ_AX done, erru=',erps 


!(((((((((((((((((((((((((((((((((((((((((
!       if(i_diag.eq.1)then
!       call psi_diag()
!       call curr_to_dina(ncequi,pjkp1,ftok)
!       call loopflux()
!       call probefield()
!       end if
!))))))))))))))))))))))))))))))))))))))))) 


          CALL eq_par( z0cen, alp, alpnew, qcen, nctrl, numlim, up,
     *                 rm,    zm,  rx0,    zx0    )

            !write(*,*) 'stepon:eq_par done, erru=',ERPS 

          IF(ngav1.EQ.2 .OR. ngav1.EQ.3)  ftok = tokout
          IF(ngav1.EQ.4 .OR. ngav1.EQ.5)  ftok = tokout
 
          !DELBET = BETPOL - BEOLD
          !DELZLI = ZLI3   - ZLOLD
          !DELRMB = RM     - RMAOLD
          !DELZMB = ZM     - ZMAOLD
          !DELRXB = RX0    - RXPOLD
          !DELZXB = ZX0    - ZXPOLD
          !DELRMA = RM     - RAXPR
          !DELZMA = ZM     - ZAXPR
          !DELRXP = RX0    - RXPPR
          !DELZXP = ZX0    - ZXPPR
              !dvz=dabs( DELZMB/(DELZMA+1.d-10) )




!          CALL  DIFGEO( RM, ZM, RAXPR, ZAXPR, RMAOLD, ZMAOLD,
!     *                  DELAXA, DELAXR, NOUT, NTER )

            !write(*,*) 'stepon:DIFGEO, erru=',ERPS 

          !BEOLD  = BETPOL
          !ZLOLD  = ZLI3
          !RMAOLD = RM
          !ZMAOLD = ZM
          !RXPOLD = RX0
          !ZXPOLD = ZX0
c
!          CALL  DIFFER( PSKP, PSKP1, NCEQUI, ERRPS1, ERRPS2,
!    *                  PSIMAX, PSIMIN, NOUT, NTER )
!          CALL  DIFTIM( PSK, PSKP, PSKP1, NCEQUI, DELPS1, DELPS2,
!     *                  DELMAX, DELMIN, NOUT, NTER )
 
            !write(*,*) 'stepon:DIFTIM, erru=',ERPS 

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c  -----  control of EXIT from iteration loop  ------
c         n_pr = 5
c           a_print(1) = knel
c         a_print(2) = erps
c           a_print(3) = delaxr
c           a_print(4) = delzma
c           a_print(5) = delzmb
c           num = 30
c           apr = 'KNEL ERPS DELAXR DELZst DELZit'
c           call out42(n_pr,a_print,num,apr)

c 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
            write(*,*) 'stepon:erru',erps 

               !IF( ERPS .LT. ENELS .AnD. dvz .LT. 5.d-3 ) THEN
               IF( erps .LT. enels  ) THEN

               !IF( dvz .LT. 5.d-3 ) THEN
                   go to 200
               END IF
               IF( knel .GE. knels ) THEN
                   keypri = 1
            write(*,*) 'spider,sstepon: no covergence ' 
            write(*,*) .GE.'KNEL  KNELS  ',knel,knels 
            pause 'pause'            
                   go to 200
               ELSE
                   go to 1000
               END IF

c 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
 
  200                          CONTINUE

        !call get_par(psi_bnd)

        diftok=dabs(platok-ztok)/(dabs(platok-ztok_n)+1.d-8)
        difpsi=dabs(psax-zpsim)/(dabs(psax-zpsim_n)+1.d-8)


c       print *,' platok ztok=',platok,ztok
c       print *,' psax zpsim=',psax,zpsim

c        diftok=dabs(platok-ztok)/(dabs(platok)+1.d-4)
c        difpsi=dabs(psax-zpsim)/(dabs(psax)+1.d-4)

                
        !write(*,*)'***it_dmf',it_dmf !
        !write(*,*)'diftok=',diftok  ! 
        !write(*,*)'difpsi=',difpsi  ! 

       !if(diftok.lt.3.5d-3 .AnD. it_dmf.gt.2) then
      !if(diftok.lt.5.0d-3 .AnD. difpsi.lt.5.0d-3 .AnD.it_dmf.gt.0) then
       if(diftok.lt.5.0d-3 .AnD. it_dmf.gt.0) then

         call eqb( alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2,
     *             betplx, i_betp,
     *             keyctr, kstep, platok, rax,zax, b0,r0, psax, igdf,
     *             n_tht, n_psi, epsro, nurs, i_eqdsk,i_bsh,
     *             psi_bnd,psi0_bnd )
         exit

         endif

         ztok  = platok
         zpsim = psax

        enddo ! it_dmf <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !print *,'stepon:eqb'

             ztok_n=platok
             zpsim_n=psax
        !print *,'ztok_n',ztok_n
        !print *,'zpsim_n',zpsim_n
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c --- printing of iteration loop accuracy
c
c      IF( keypri.EQ.1 ) THEN
c          WRITE(nout,*) '_______________________________________'
c          WRITE(nout,1915) kstep, knel, timev
c          WRITE(nout,*) 'ERPS   =', erps
c          WRITE(nout,*) 'ERRCU1 =', errcu1,' ERRCU2 =', errcu2
c...      WRITE(nout,*) 'CURMAX =', curmax,' CURMIN =', curmin
c          WRITE(nout,*) 'ERRPS1 =', errps1,' ERRPS2 =', errps2
c...      WRITE(nout,*) 'PSIMAX =', psimax,' PSIMIN =', psimin
c          WRITE(nout,*) 'DELPS1 =', delps1,' DELPS2 =', delps2
c...      WRITE(nout,*) 'DELMAX =', delmax,' DELMIN =', delmin
c          WRITE(nout,*) 'DELAXA =', delaxa,' DELAXR =', delaxr
c      END IF
c
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       CALL  bp_prob( nout,   nter,   keypri,
c     *                npro,   rpro,   zpro,   fipro,
c     *                brpro2, bzpro2, bpcom2 )
c***********************************************************************
c   compute full poloidal flux in given "loop" points
c
c       CALL  fl_loo( nout,   nter,   keypri,
c     *               nloo,   rloo,   zloo,
c     *               psloo2 )
c***********************************************************************
c
c.....  toroidal currents in passive conductor structures 
c
!        CALL PASCUR( NOUT, NTER, NEQUI, NFW, NBP, NVV,
!     *               PJKP1, fwcurr, bpcurr, vvcurr )
        !print *,'stepon:pascur'
c-----------------------------------------------------------------------
c --- printing
c 
          psibou = up
          psidel = psiout - psibou


c
cccc          include 'printa.inc'
c
c***********************************************************************
c--- time level parameters are written in  disk  files
c
         !!! include 'wrt2.inc'
c***********************************************************************
  222        CONTINUE
c
             errcu1 = 0.d0
             errcu2 = 0.d0
             errps1 = 0.d0
             errps2 = 0.d0
             delps1 = 0.d0
             delps2 = 0.d0
             ereve  = 0.d0
             erps   = 0.d0
c
             raxpr  = rm
             zaxpr  = zm
            ! RXPPR  = RX0
            ! ZXPPR  = ZX0
            kstepr  = kstep
            tstepr  = tstep
c
          DO 293 l=1,npfc
             pfcur1(l) = pfcur2(l)
             pfcw1(l)  = pfcw2(l)
             pfcd1(l)  = pfcd2(l)
  293     CONTINUE

          do i=1,nequi
             d_pf_mat(i) = pjkp1(i)*1.0d6  !!!  in [A]
          enddo
 
             camtok=0.d0
         do i=nequi+1,ncequi
             d_tcam_mat(i-nequi) = pjkp1(i)*1.0d6  !!!  in [A]
             camtok=camtok+pjkp1(i)
        enddo


          DO 260 l=1,ncequi
             pjk(l)   = pjkp1(l)
             volk(l)  = volkp1(l)
             pskm1(l) = psk(l)
             psk(l)   = pskp1(l)
  260     CONTINUE
c
c***********************************************************************
        if(kstep.lt.nstep_p)then
           i_tim=kstep

         time_t(i_tim)=time
         torcur(i_tim)=platok
         rm_t(i_tim)=rax
         zm_t(i_tim)=zax
         rxp_t(i_tim)=rxpnt
         zxp_t(i_tim)=zxpnt
         betpol_t(i_tim)=betpol
         bettor_t(i_tim)=betful
         psim_t(i_tim)=psax
         psib_t(i_tim)=psbo
        endif
c-----------------------------------------------------------------------
c for writing for "motion picture"
c
c       psiplb(kstep)=pspl_av
c       psiexb(kstep)=psex_av
c       psimag(kstep)=psi0_ax
c       flu_tor(kstep)=flx_fi
c
c       if(kstep.gt.0 .AnD. kstep.le.2000) then
c       kskw=2
c       if(kstep/kskw*kskw .eq. kstep) then
c         numwr=numwr+1
c         call wrdfmv(numwr,timev)
c         call wrdump(numwr,timev,kstep,psiplb,psiexb,psimag,flu_tor)
c       endif
c       endif
c 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
!!!!                   GO TO 100  
c 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100
c-----------------------------------------------------------------------
  103  FORMAT(2x,8e12.5)
 1915  FORMAT(3x,'LEVEL K =',i3,2x,'ITERATION KNEL =',i3,2x,
     *        'TIME(K) =',e12.5)
c********************************


      RETURN
      END