add_resistivity.f90

 
      SUBROUTINE add_resistive_e
      USE stel_kinds
      USE shared_data, ONLY: l_getfsq, l_init_state, lverbose,          &
                             l_update_state, unit_out, nprecon,         &
                             fsq_total1, xc, gc
      USE shared_functions, ONLY: funct_island
      USE siesta_namelist, ONLY: ftol, eta_factor
      USE descriptor_mod, ONLY: iam
      USE siesta_bfield, ONLY: update_bfield
      USE siesta_displacement, ONLY: update_upperv
      USE siesta_init, ONLY: init_state
      USE siesta_state, ONLY: update_state
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      INTEGER, PARAMETER   :: nits = 10
      REAL(dp), PARAMETER  :: zero = 0
      REAL(dp)             :: save_eta, fmhd, eps
      INTEGER              :: itime, j
      LOGICAL, PARAMETER   :: lcurrent_only=.true.
      LOGICAL              :: lverb_save
      
      INTEGER, SAVE        :: nCheck = 0
!-----------------------------------------------
      IF (fsq_total1 .LE. ftol) RETURN
 
      lverb_save = lverbose
      lverbose = .false.
      
      fmhd = fsq_total1
 
      l_update_state = .false.  
      
!     No perturbation, vsupX = 0
      xc = 0
      CALL init_state(.false.)
      CALL update_upperv
 
      save_eta = eta_factor
      eta_factor = eta_factor/max(1,nits)
 
!     Diffuse B-field (eta*j) from last update_state call
      DO itime=1,nits
 
         CALL init_state (lcurrent_only)                                 !Use new B-field to get KsubXij currents
         CALL update_bfield (.true.)                                     !Resistive update to B-field
         CALL update_state (.false., zero, zero)                         !Updates B-field with resistive piece
 
      END DO
 
      ncheck = ncheck+1
      IF (ncheck .EQ. 0) CALL checkforces(xc, gc)                        !set iteration value > 0 at which force check is done
 
      l_getfsq = .true.
      l_init_state = .true.
      CALL funct_island
      
      
      eta_factor = save_eta
      lverbose = lverb_save
      
      IF (iam .EQ. 0) THEN
         DO j = 6, unit_out, unit_out-6
         IF (.NOT.lverbose .AND. j.EQ.6) cycle
            WRITE (j,'(/,a,i3)') ' UPDATING RESISTIVE E-FIELD: ITERATIONS=',nits
            WRITE (j, '(a,1p2e12.3)')                                   &
          ' JUMP IN FSQ DUE TO RESISTIVE DIFFUSION: ', fmhd, fsq_total1
         END DO
      END IF
 
!SPH101116: limit jump next time
      IF (fsq_total1 .GT. 4*fmhd) THEN
         eta_factor = eta_factor*(4*fmhd/fsq_total1)
      ELSE IF (fsq_total1 .LT. 1.1_dp*fmhd) THEN
         eta_factor = eta_factor*(1.1_dp*fmhd/fsq_total1)
      END IF
      
      END SUBROUTINE add_resistive_e