!
!......................................................................
      SUBROUTINE pstnode2
!......................................................................
!       version of nodes to permit arbitrary distribution of fin. elements
!       the plasma is divided into regions by the singular surfaces.
!       maximum allowed is 5 regions ie. 4 sing. surfaces.
!       each region is divided into zones where the nodeds are either
!       distributed uniformly, or packed with a power law around the
!       singular surfaces.
!
!       modified version of nodes to allow packing of mesh at several
!        resonant surfaces. mesh distribution is specified by:
!       alfmsh(5)  exponent of mesh scaling around magnetic axis and
!                  rational surfaces.
!       widmsh(5)  width of scaled zone.
!       distribution function of mesh nodes is
!  -/+zg(i)*|(z-zts(i))/widmsh(i)|**alfmsh(i)+zd(i)    in scaled zone.
!       zb*z + zc(i)       in linear zone,
!       for z = inode-1/m,
!       and i=1,nsing+ 1._r8   zts, zg, zd, zb and zc are computed
!       to ensure continuity of distribution function and its
!       derivative. a. pletzer, 2 july  90._r8 
!
! mofifications:
!
!         uniform distribution of z within each section. aplet  20.4_r8 .93
!         introduce distortion function to increase node density
!         near plasma edge.                              aplet  6.5_r8 .93
!......................
 USE pstcom

 USE temps

 USE comggf
 IMPLICIT NONE
  integer,parameter :: r8 = SELECTED_REAL_KIND(12,100)
      INTEGER N2
      INTEGER MTOT
      INTEGER JTEST
      INTEGER J
      INTEGER JMAX
      INTEGER J2
      INTEGER MMINUS
      INTEGER I
      INTEGER JFLAG2
      INTEGER JFMAX
      real*8 SUM
      real*8 B
      real*8 Z
      INTEGER JFLAG
      real*8 Z1
      real*8 Z2
      INTEGER I3
      INTEGER I0
      INTEGER I1
      INTEGER I2
      real*8 ZM
      INTEGER JS
      real*8 ZDXS
      INTEGER IN
      INTEGER IN2
      real*8 ZX
      real*8 DELPSL
      real*8 DELPSR
      real*8 DELPSI
      real*8 DELPSA
      INTEGER KN
      INTEGER KM
      real*8 DPS
      INTEGER IK
      real*8 ZL
      real*8 ZR



!
 LOGICAL	 	 	 :: nodchk
 REAL*8, DIMENSION(nfe) 	 :: width
 REAL*8, DIMENSION(nsg1) 	 :: zc
 REAL*8, DIMENSION(nsg1) 	 :: zg
 REAL*8, DIMENSION(nsg1) 	 :: zd
 REAL*8, DIMENSION(nsg1) 	 :: zts
 REAL*8, DIMENSION(nsg1) 	 :: zsing
 REAL*8, DIMENSION(nsg1,nsg1) 	 :: zonew
 REAL*8	 	 	 :: arad
 REAL*8, DIMENSION(nsg) 	 :: xrat
 INTEGER, DIMENSION(55) 	 :: iord
 INTEGER, DIMENSION(5) 	 :: lzone
 INTEGER, DIMENSION(2) 	 :: iop
 INTEGER, DIMENSION(3,nsg1) 	 :: mzoned
 INTEGER, DIMENSION(2) 	 :: minner
 COMMON /c2f90nodchk/ nodchk 
 COMMON /c2f90zonew/ zonew 
 COMMON /c2f90mzoned/ mzoned 
 COMMON /c2f90minner/ minner 
 COMMON /c2f90arad/ arad 
 COMMON /c2f90xrat/ xrat 

 integer index_ratsurf
 real*8 zdiav

      if( imode == 1 ) then
      write(itty,8756) nosurf,dlayb, dlay
      write(outmod,8756) nosurf,dlayb, dlay
!
!
 8756 format(/' equilibrium grid size         nosurf = ',i7/  &
       ' rel. width of big element      dlayb = ',f13.6,/  &
       ' rel. width for frobenius coeff. dlay = ',f13.6/)
      end if
!
      write(itty,8755) m
      write(outmod,8755 ) m
 8755 format(/' number of finite elements          m = ',i7/)
!
      msing(1) = 1
      nusurf   = nosurf
      mtot = m
!
      psidif = psia(nosurf) - psia(1)
      mminus = m - 2*nosing
!
! test for presence of singular surfaces...
!
      if( nosing ==  0 ) then
         !
         ! linear mesh...
         !
         psinod = (/ (psidif*(i-1)/m, i=1,mtot+1) /) + psia(1)
         !
         go to 101
      end if
!
! there are singular surface(s)...

!     reset arrays...
!
      alfmsh = 1.0_r8
      widmsh(nsing1+1) =  0.0_r8 

! compute mesh exponents
      
      if(lpack) then
         alfmsh(1) = 2.0_r8
         if (nsing1>2) then
            do i = 2, nsing1 
               if(abs(xmu(i-1)-0.5_r8)>1.e-6_r8) alfmsh(i) = 1.0_r8/xmu(i-1)
            enddo
         endif
      endif
!!$         do i = 1, nosing
!!$            index_ratsurf = 1
!!$            do in = 1, nosurf
!!$               if(psisin(i+1) > psia(in)) &
!!$                    & index_ratsurf = in
!!$            enddo
!!$            zdiav = (di(index_ratsurf)+di(index_ratsurf+1))/2._r8
!!$            if(abs(zdiav+0.25_r8) > 1.e-6_r8) then
!!$               alfmsh(i+1) = 1.0_r8/sqrt(-zdiav)
!!$            else
!!$               lzerop = .true.
!!$               alfmsh(i+1) = 1.0_r8
!!$            endif
!!$         enddo
!!$
!!$      endif
         
!
!...............................................
!
       jflag2 = 0
       jfmax = 10
!
 100   continue
!
!       check alfmsh and widmsh...
!
       sum = - widmsh(1)
       do 80 i=1,nsing1
       alfmsh(i) = max(1.0_r8, abs( alfmsh(i) ) )
       widmsh(i) = abs( widmsh(i) )
       sum = sum +  2._r8 *widmsh(i)
 80     continue
!
       if( sum >=    1._r8  ) then
       write(itty,*)  ' *** widmsh too large'
       write(outmod,*)' *** widmsh too large'
       do 81 i=1,nsing1
       widmsh(i) = widmsh(i)/(sum+ 0.1_r8 )
 81     continue
       write(itty,*)'corrected values are:'
       write(outmod,*)'corrected values are:'
       write(itty,*)'widmsh = ',(widmsh(i),i=1,nsing1) 
       write(outmod,*)'widmsh = ',(widmsh(i),i=1,nsing1) 
       end if
!
!       position of rational surfaces...
!
       do 82 i=1,nsing1+1
       zsing(i) = (psisin(i) - psia(1))/psidif
 82     continue
!
!       slope of mesh distribution function in linear zones...
!
       b =  1._r8  - (alfmsh(1)- 1._r8 )*widmsh(1)/alfmsh(1)
       do i=2,nsing1
          b = b -  2._r8 *(alfmsh(i)- 1._r8 )*widmsh(i)/alfmsh(i)
       end do
       b =  1._r8 /b
!
!       first zone near magnetic axis...
!
       zd(1) =  0._r8 
       zg(1) = b*widmsh(1)/alfmsh(1)
!
!       first linear mesh zone...
!
       zc(1) = - b*widmsh(1)*(alfmsh(1)-1)/alfmsh(1)
!
!       iterate over regions...
!
       zts(1) =  0._r8 
       do i=2,nsing1
          zc(i) = zc(i-1) -  2._r8 *b*(alfmsh(i)-1)*widmsh(i)/alfmsh(i)
          zg(i) = b*widmsh(i)/alfmsh(i)
          !
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$          if(  1._r8  - zsing(i) >   1.E-8_r8  ) then
!!$          zd(i) =  1._r8  - (  1._r8  - zsing(i))**( 1._r8 /betmsh)
!!$          else
!!$          zd(i) =  1._r8 
!!$          end if
          zd(i) = zsing(i)
          !
          !aplet  22.4_r8 .98 zts(i) = ( zd(i) - (zc(i)+zc(i-1))/ 2._r8  )/b
          zts(i) = zts(i-1) - zsing(i-1)/b + zsing(i)/b &
               + widmsh(i-1)*(alfmsh(i-1)- 1._r8 )/alfmsh(i-1) &
               + widmsh(i  )*(alfmsh(i  )- 1._r8 )/alfmsh(i  )
       end do
       zts(nsing1 + 1 ) =  1._r8 
!
!       check...
!
       z = b + zc(nsing1)
       if( abs(z- 1._r8 ) >    1.E-10_r8 ) then
       write(itty, * )  ' ***error in node2: zc(nosing)+b = ; ',z
       write(outmod, * )' ***error in node2: zc(nosing)+b = ; ',z
       end if
!
       jflag = 0
       do 88 i=1,nsing1
       z1 = zts(i+1) - zts(i)
       z2 = widmsh(i) + widmsh(i+1)
       if( z1 < z2) then
          write(itty,*)' *** i/M distance between rational surfaces = ', &
               z1, ' < widmsh(',i+1,') + widmsh(',i,') = ' , z2
          write(itty,*)' *** check msin and/or decrease widmsh '
          write(outmod,*)' *** i/M distance between rational surfaces = ', &
               z1, ' < widmsh(',i+1,') + widmsh(',i,') = ' , z2
          write(outmod,*)' *** check msin and/or decrease widmsh '
       widmsh(i) = widmsh(i)/(z2/z1 +  0.1_r8 )
       widmsh(i+1) = widmsh(i+1)/(z2/z1 +  0.1_r8 )
       jflag = jflag + 1
       end if
  88   continue
!
       jflag2 = jflag2 + 1 
       if( jflag /= 0  .AND.  jflag2 <= jfmax ) go to 100
!
!
!      construct psi values at nodes
!
!
!       scan over regions...
!
       i3 = 0
       msing(1) = 1
          do 85 i = 1,nsing1
          i0 = i3 + 1
          i1 = i0 + int( widmsh(i) * real(mminus) )
          i3 = int( zts(i + 1) *real(mminus) )  + 1
          msing(i+1) = i3 + 1 + 2*(i-1)
          if( i == nsing1 ) msing(i+1) = mtot + 1
          i2 = i3 - int( widmsh(i+1) *real(mminus))
          if( i ==  nsing1 ) i2 = i3
!
          mzoned(1,i) = i1 - i0 + 1
          if( i /=  1 ) mzoned(1,i) = i1 - i0 + 2
          mzoned(2,i) = i2 - i1
          mzoned(3,i) = i3 - i2
          if( i /= nsing1 ) mzoned(3,i) = i3 - i2 + 1
          lzone(i) = 3
          if( i ==  nsing1 ) lzone(i) = 2
!
       zm = real( i3-i0+1 )
       if( (i /= 1)  .AND.  (i /= nsing1) ) zm = zm +  1.0_r8 
       js = i0 - 1
       if( i == 1 ) js = js + 1
       zdxs = zts(i+1) - zts(i)
!
       do 86 in = i0,i1
          in2 = in + (i-1)*2
          z = zts(i) + zdxs*( real(in-js)/zm )
          zx = +zg(i  )*abs( (z-zts(i  ))/widmsh(i  ) )**alfmsh(i  ) &
               + zd(i  )
          !
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$          if(  1._r8  - zx >   1.E-8_r8  ) then
!!$          psinod(in2) = ( 1._r8 -( 1._r8 -zx)**betmsh )* psidif + psia(1)
!!$          else
!!$            psinod(in2) = psidif + psia(1)
!!$          end if
          psinod(in2) = zx *psidif + psia(1)
          !
86        continue
!
       do 87 in = i1+1,i2
       in2 = in + (i-1)*2
       z = zts(i) + zdxs*( real(in-js)/zm )
       zx =  b*z + zc(i)
!
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$       if(  1._r8  - zx >   1.E-8_r8  ) then
!!$       psinod(in2) = ( 1._r8 -( 1._r8 -zx)**betmsh )*psidif+psia(1)
!!$       else
!!$       psinod(in2) = psidif+psia(1)
!!$       end if
       psinod(in2) = zx *psidif + psia(1)
!
 87    continue
!
       do 89 in = i2+1,i3
       in2 = in + (i-1)*2
       z = zts(i) + zdxs*( real(in-js)/zm )
       zx = -zg(i+1)*abs( (z-zts(i+1))/widmsh(i+1) )**alfmsh(i+1) &
       + zd(i+1)
!
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!$       if(  1._r8  - zx >   1.E-8_r8  ) then
!!$       psinod(in2) = ( 1._r8 -( 1._r8 -zx)**betmsh)*psidif+psia(1)
!!$       else
!!$       psinod(in2) = psidif+psia(1)
!!$       end if
       psinod(in2) = zx * psidif + psia(1)
!
 89    continue
!
 85       continue
!
!----------------------------------------------------------------
!
 101   continue
!
!     construct new grid of psi surfaces for remapping onto later.
!
       msub = 2
!!$       if ( .not. lsub )  go to 62
!!$!
!!$!      default to one grid point between each node...
!!$!
!!$       do i = 1, mtot+1
!!$          msub(i) = 2
!!$       end do
!!$ 62    continue
!
       if( nosing /= 0 ) then
!
! add 2 equidistant surfaces on either side of rational surface.
!
      do 90 i = 1,nosing
      delpsl = psisin( i+1 ) - psinod( msing(i+1) - 1 )
      delpsr = psinod( msing(i+1) + 2 ) - psisin( i+1 )
      delpsi = delpsl/ 2._r8 
      delpsa = abs( dlay * psidif )
      if( delpsr < delpsl ) delpsi = delpsr/ 2._r8 
      if( delpsa  <=  delpsi ) then
         delpsi = delpsa
      else
         if( imode == 1 ) then
            write(itty,  *)' ***node2: delpsa  >   ' &
               ,'delpsi between psisin(',i,') and psisin(',i+1,')'
            write(itty,  *)' should decrease dlay'
            write(outmod,  *)' ***node2: delpsa  >   ' &
               ,'delpsi between psisin(',i,') and psisin(',i+1,')'
            write(outmod,  *)' should decrease dlay'
         end if
      end if
!
      psinod( msing(i+1) ) = psisin(i+1) - delpsi
      psinod( msing(i+1) + 1 ) = psisin(i+1) + delpsi
!
       if( imode == 1 ) then
      write(itty,8757)   delpsi 
      write(outmod,8757) delpsi 
 8757 format( ' half width of element near rat. surf = ',f13.8/)
       end if
!
 90   continue
!
      end if
!
 
!       reset m to reflect the new number of surfaces
!
       m = mtot
       mp = m + 1

!!$       if(mp >= nfe) then
!!$          print*,'***ERROR** Total number of finite elements ', mp,' >= grid resolution set by nfe=',nfe
!!$          stop
!!$       endif
!
       kn = 1
       km = 1
       psinew(1) = psinod(1)
        do 65 i = 1, mtot
       kn = kn + 1
       dps = ( psinod(kn)-psinod(kn-1) ) / msub(i)
       dtent(kn-1) = dps * msub(i)
       do 63 ik = 1, msub(i)
       km = km + 1
       psinew(km) = psinod(kn-1) + ik * dps
 63    continue
 65    continue
!
!      set nusurf
!
       nusurf = km
!
!      for diagnostics
!
      write(outmod,9007)(psisin(i),i=1,nsing1+1)
      write(outmod,9015)mtot
!
      write(itty  ,9017)
      write(outmod,9017)
      do i = 1, nsing1
      write(itty  ,9018) i, mzoned(1,i), mzoned(2,i), mzoned(3,i)
      write(outmod,9018) i, mzoned(1,i), mzoned(2,i), mzoned(3,i)
      end do	  
!
! set up table for plotting mesh characteristics
!
!!$      write(outnod,9008)
!!$      zl =  0._r8 
!!$      zr =  0._r8 
!!$      do 51 i = 1,mtot + 1
!!$      write(outnod,9020) i, psinod(i), zl, zr
!!$      zl =  2.0_r8  * abs( real((i+1)/2) - real(i+1)/ 2._r8  )
!!$      zr =  2.0_r8  * abs( real(i/2) - real(i)/ 2._r8  )
!!$ 51   continue
!
      write(outmod,9010)(psinew(i),i=1,nusurf)
 9007 format(/,10(' psisin = ',7e12.5))
 9008 format(2x,'node',4x,'psinod',5x,'#1',2x,'#2')
 9020 format(2x,i4,e14.5,f4.0,f4.0)
 9010 format(' psinew = ',8e12.5)
!
!
 9012 format('this region is divided into ',i3,' zones',/ &
 ,' zone no.',   2x,' zone width',8x,'in psi',6x,' fin. elements',2x,' dpsi',/)
 9013 format(4x,i3,4x,f10.3,3x,e15.5,9x,i3,4x,e12.5)
 9014 format(/,' region no ',i3,' from ',e15.5,' to ',e15.5,/)
 9015 format(/,' there are', i5, ' finite elements in the mesh',/)
!
 9016 format('error in distributing the finite elements, check mzoned'    ,/ &
,' mtot = ',i5,' mzoned= ',5(11i4))
!
 9017 format(//' number of mesh nodes in:'/ &
         ' section     zone-1     zone-2     zone-3')
 9018 format(3x,i2,3x,3(8x,i3))
!
       if ( check1 ) then
          print *,'psia'
          print *,psia
          print *,'psinod'
          print *,psinod
          print *,'psinew'
          print *,psinew
       end if
          
!
9005     format(2x,' surf    psi-old      psi-new',/, &
     130(1x,i5,2e14.5,/) )
!
! check if psinew is in increasing order...
!
      jtest = 0
      do j = 1, nusurf-1
      if( psinew(j+1) <= psinew(j) ) then
      jtest = jtest+1
      write(  itty,*)' *** node2: psinew not in increasing order' &
  , ' about node ',j,' psinew = ',psinew(j)
      write(outmod,*)' *** node2: psinew not in increasing order' &
  , ' about node ',j,' psinew = ',psinew(j)
      write(  itty,*)'will stop '
      write(outmod,*)'will stop '
      stop ' ERROR: psinew not in increasing order in node2'
      end if
      end do
!
       return
       end subroutine



