ReadNEMEC – GVEC

public subroutine ReadNEMEC(fileName, itype, ok)

Uses

- MODgvec_Globals
Help

Graph Key

Nodes of different colours represent the following:

Solid arrows point from a submodule to the (sub)module which it is descended from. Dashed arrows point from a module or program unit to modules which it uses.

READ VMEC "wout" datafile generated by NEMEC, routine provided by Erika Strumberger, IPP Garching can be either ascii or binary ! !test output !! WRITE(UNIT_StdOut,"(6x,'hsve',10x'hiota',9x,'hmass',9x,'hpres',9x,'hphip', & !! 9x,'hbuco',9x,'hbvco',10x,'fphi',11x,'hvp',9x,'hoverr', & !! 8x,'hspecw')") ! DO j=1,nsin ! WRITE(UNIT_StdOut,"(13(2x,e12.4))") hsve(j),hiota(j),hmass(j),hpres(j),hphip(j), & ! hbuco(j),hbvco(j),fphi(j),hvp(j),hoverr(j),hspecw(j) ! END DO ! WRITE(UNIT_StdOut,"(6x,'fsve',10x,'fjcuru',8x,'fjcurv')") ! DO j=1,nsin ! WRITE(UNIT_StdOut,"(3(2x,e12.4))") fsve(j),fjcuru(j),fjcurv(j) ! END DO

Arguments

Type	Intent		Name
character(len=*),	intent(in)	::	fileName
integer,	intent(in)	::	itype	=0: binary, =1: ascii
integer,	intent(out)	::	ok

Calls

Help

Called by

Help

Source Code

SUBROUTINE ReadNEMEC(fileName,itype,ok)
  USE MODgvec_globals,ONLY:GETFREEUNIT
  IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT/OUTPUT VARIABLES
  CHARACTER(LEN = *), INTENT(IN) :: fileName
  INTEGER, INTENT(IN)           :: itype  !! =0: binary, =1: ascii
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
  INTEGER, INTENT(OUT)           :: ok
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
  INTEGER              :: nrho,mpnt,nsin !nfp
  INTEGER              :: mpol1,iasym
  INTEGER              :: inUnit,iMode
  INTEGER              :: ierr,j,m,n,nmin0
  REAL(wp)             :: enfp,enrho,empol,entor,empnt,eiasym
  REAL(wp)             :: ds
  REAL(wp)             :: gam,phiedge

  REAL(wp),ALLOCATABLE :: hiota(:),hpres(:),hbuco(:),hbvco(:)
  REAL(wp),ALLOCATABLE :: hmass(:),hphip(:),fphi(:),hvp(:)
  REAL(wp),ALLOCATABLE :: hoverr(:),fjcuru(:),fjcurv(:),hspecw(:)

  REAL(wp),ALLOCATABLE :: frmnc(:,:,:),frmns(:,:,:)
  REAL(wp),ALLOCATABLE :: fzmnc(:,:,:),fzmns(:,:,:)
  REAL(wp),ALLOCATABLE :: hbsmnc_dw(:,:,:),hbsmns_dw(:,:,:)
  REAL(wp),ALLOCATABLE :: hbumnc_dw(:,:,:),hbumns_dw(:,:,:)
  REAL(wp),ALLOCATABLE :: hbvmnc_dw(:,:,:),hbvmns_dw(:,:,:)
  REAL(wp),ALLOCATABLE :: hbumnc_up(:,:,:),hbumns_up(:,:,:)
  REAL(wp),ALLOCATABLE :: hbvmnc_up(:,:,:),hbvmns_up(:,:,:)
  REAL(wp),ALLOCATABLE :: flmnc(:,:,:),flmns(:,:,:)

  REAL(wp),ALLOCATABLE :: fsve(:),hsve(:)    ! radial mesh
!===================================================================================================================================
  inUnit=GETFREEUNIT()
  ok = 0

  ! test input
  IF(itype.LT.0)THEN
    WRITE(UNIT_stdOut,*) '********** USER error **********',itype
    WRITE(UNIT_stdOut,*)'problem reading file:',TRIM(filename)
    ok = -1 ; RETURN
  END IF

  ! open equilibrium file
  IF(itype.EQ.0) THEN
    OPEN(inUnit,file=trim(filename),form='unformatted', status='old',iostat=ierr)
  ELSE
    OPEN(inUnit,file=trim(filename),form='formatted', status='old',iostat=ierr)
  END IF
  IF(ierr.NE.0) THEN
    WRITE(UNIT_stdOut,*) '********** USER error **********'
    WRITE(UNIT_stdOut,*)'ierr = ',ierr
    WRITE(UNIT_stdOut,*)'could not open file:',TRIM(filename)
    ok = -1 ; RETURN
  END IF

! --- read dimensions
  IF(itype.EQ.0) THEN
    READ(inUnit) gam,enfp,enrho,empol,entor,empnt,eiasym,phiedge
  ELSE
    READ(inUnit,*) gam,enfp,enrho,empol,entor,empnt,eiasym,phiedge
  END IF
  nfp    = NINT(enfp)
  nrho   = NINT(enrho)
  mpol   = NINT(empol)
  ntor   = NINT(entor)
  mpnt   = NINT(empnt)
  iasym  = NINT(eiasym)

  WRITE(UNIT_stdOut,'(6X,A)')        '-----------------------------------'
  WRITE(UNIT_stdOut,'(6X,A)')        'NEMEC file parameters:'
  WRITE(UNIT_stdOut,'(6X,A,I6)')     '  nfp     : ',nfp
  WRITE(UNIT_stdOut,'(6X,A,I6)')     '  nrho    : ',nrho
  WRITE(UNIT_stdOut,'(6X,A,I6)')     '  mpol    : ',mpol
  WRITE(UNIT_stdOut,'(6X,A,I6)')     '  ntor    : ',ntor
  WRITE(UNIT_stdOut,'(6X,A,E23.15)') '  gamma   : ',gam
  WRITE(UNIT_stdOut,'(6X,A,E23.15)') '  phiEdge : ',phiEdge
  IF(iasym.EQ.0)THEN
    WRITE(UNIT_stdOut,'(6X,A,I6)')   '  iasym=0... SYMMETRIC (half fourier modes)'
  ELSE
    WRITE(UNIT_stdOut,'(6X,A,I6)')   '  iasym=1... ASYMMETRIC'
  END IF
  WRITE(UNIT_stdOut,'(6X,A)')        '-----------------------------------'

  nsin   = nrho-1
  mpol1  = mpol-1
  ds     = 1./REAL(nsin,wp)

  ALLOCATE(frmnc(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(frmns(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(fzmnc(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(fzmns(0:mpol1,-ntor:ntor,0:nsin))

  ALLOCATE(hbumnc_up(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbumns_up(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbvmnc_up(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbvmns_up(0:mpol1,-ntor:ntor,0:nsin))

  ALLOCATE(hbsmnc_dw(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbsmns_dw(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbumnc_dw(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbumns_dw(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbvmnc_dw(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(hbvmns_dw(0:mpol1,-ntor:ntor,0:nsin))

  ALLOCATE(flmnc(0:mpol1,-ntor:ntor,0:nsin))
  ALLOCATE(flmns(0:mpol1,-ntor:ntor,0:nsin))

  ALLOCATE(hiota(1:nsin),hpres(1:nsin),hbuco(1:nsin),hbvco(1:nsin))
  ALLOCATE(hmass(1:nsin),hphip(1:nsin),fphi(1:nsin),hvp(1:nsin))
  ALLOCATE(hoverr(1:nsin),fjcuru(1:nsin),fjcurv(1:nsin))
  ALLOCATE(hspecw(1:nsin))

  ALLOCATE(fsve(0:nsin),hsve(1:nsin))

! --- read NEMEC output
  IF(itype.EQ.0) THEN

    DO j=0,nsin
      DO m=0,mpol1
        nmin0=-ntor
        IF(m.EQ.0) nmin0=0
        DO n=nmin0,ntor
! --- full mesh
          READ(inUnit) frmnc(m,n,j),fzmns(m,n,j),         &
                      frmns(m,n,j),fzmnc(m,n,j),         &
! --- half mesh
                      hbumnc_up(m,n,j),hbvmnc_up(m,n,j), &
                      hbumns_up(m,n,j),hbvmns_up(m,n,j), &
! --- full mesh
                      flmns(m,n,j),flmnc(m,n,j),         &
! --- half mesh
                      hbumnc_dw(m,n,j),hbvmnc_dw(m,n,j), &
                      hbsmns_dw(m,n,j),                  &
                      hbumns_dw(m,n,j),hbvmns_dw(m,n,j), &
                      hbsmnc_dw(m,n,j)
        END DO
      END DO
    END DO

! --- half mesh
    READ(inUnit) (hiota(j),hmass(j),hpres(j),hphip(j),hbuco(j), &
                hbvco(j),fphi(j),hvp(j),hoverr(j),fjcuru(j),   &
                fjcurv(j),hspecw(j),j=1,nsin)
  ELSE !itype = 1

    DO j=0,nsin
      DO m=0,mpol1
        nmin0=-ntor
        IF(m.EQ.0) nmin0=0
        DO n=nmin0,ntor
! --- full mesh
          READ(inUnit,*) frmnc(m,n,j),fzmns(m,n,j),       &
                      frmns(m,n,j),fzmnc(m,n,j),         &
! --- half mesh
                      hbumnc_up(m,n,j),hbvmnc_up(m,n,j), &
                      hbumns_up(m,n,j),hbvmns_up(m,n,j), &
! --- full mesh
                      flmns(m,n,j),flmnc(m,n,j),         &
! --- half mesh
                      hbumnc_dw(m,n,j),hbvmnc_dw(m,n,j), &
                      hbsmns_dw(m,n,j),                  &
                      hbumns_dw(m,n,j),hbvmns_dw(m,n,j), &
                      hbsmnc_dw(m,n,j)
        END DO
      END DO
    END DO

! --- half mesh
    READ(inUnit,*) (hiota(j),hmass(j),hpres(j),hphip(j),hbuco(j),   &
                  hbvco(j),fphi(j),hvp(j),hoverr(j),fjcuru(j),     &
                  fjcurv(j),hspecw(j),j=1,nsin)
  END IF

  CLOSE(inUnit)

! --- grid in s
  fsve(0)     = 0.
  fsve(nsin)  = 1.
  DO j=1,nsin-1
! --- half mesh
    hsve(j)  = (j-0.5)*ds
! --- full mesh
    fsve(j)   = j*ds
  END DO
  hsve(nsin) = (nsin-0.5)*ds

!!!  !test output
!!!!  WRITE(UNIT_StdOut,"(6x,'hsve',10x'hiota',9x,'hmass',9x,'hpres',9x,'hphip', &
!!!!                      9x,'hbuco',9x,'hbvco',10x,'fphi',11x,'hvp',9x,'hoverr', &
!!!!                      8x,'hspecw')")
!!!  DO j=1,nsin
!!!    WRITE(UNIT_StdOut,"(13(2x,e12.4))") hsve(j),hiota(j),hmass(j),hpres(j),hphip(j), &
!!!                                        hbuco(j),hbvco(j),fphi(j),hvp(j),hoverr(j),hspecw(j)
!!!  END DO
!!!  WRITE(UNIT_StdOut,"(6x,'fsve',10x,'fjcuru',8x,'fjcurv')")
!!!  DO j=1,nsin
!!!    WRITE(UNIT_StdOut,"(3(2x,e12.4))") fsve(j),fjcuru(j),fjcurv(j)
!!!  END DO

! --- translate to GVEC datastructure

  IF(gam.GT.1.0e-04) &
    CALL abort(__STAMP__, &
                      "readNEMEC: currently, only gamma=0 is supported")

  IF(SUM(ABS(fsve(1:nsin)*phiEdge-fphi(1:nsin)))/REAL(nsin,wp).GT.1.0e-08) &
    CALL abort(__STAMP__, &
                      "readNEMEC: toroidal flux does not seem equidistant...")

  !not needed further
  DEALLOCATE(hbumnc_up,hbumns_up,hbvmnc_up,hbvmns_up)
  DEALLOCATE(hbsmnc_dw,hbsmns_dw,hbumnc_dw,hbumns_dw,hbvmnc_dw,hbvmns_dw)
  DEALLOCATE(hbuco,hbvco,hmass,hphip,hvp,hoverr,fphi,fjcuru,fjcurv,hspecw)


  nFluxVMEC=nrho

  mn_mode = ntor+1 + (mpol-1)*(2*ntor+1)

  mn_mode_nyq=mn_mode !nyquist do not exist here

  lasym = (iasym.EQ.1)

  CALL alloc_all() !needs nFluxVMEC,mn_mode,mn_mode_nyq and lasym

  ! parameters that are not read?! set to zero here
  b0=0.
  aMinor=0.
  rMajor=0.
  volume=0.
  gmnc=0.
  bmnc=0.

  !use same loop as for read in to set modes m,n
  iMode=0
  DO m=0,mpol1
    nmin0=-ntor
    IF(m.EQ.0) nmin0=0
    DO n=nmin0,ntor
      iMode=iMode+1
      xm(iMode)=REAL(m,wp) ; xn(iMode)=REAL(n,wp)
    END DO !n
  END DO !m
  IF(iMode.NE.mn_mode) STOP 'mn_mode not correct'

  lambda_grid="full"

  DO iMode=1,mn_mode
    rmnc(iMode,1:nFluxVMEC)=frmnc(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
    zmns(iMode,1:nFluxVMEC)=fzmns(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
    lmns(iMode,1:nFluxVMEC)=flmns(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
  END DO
  IF(lasym)THEN
    DO iMode=1,mn_mode
      rmns(iMode,1:nFluxVMEC)=frmns(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
      zmnc(iMode,1:nFluxVMEC)=fzmnc(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
      lmnc(iMode,1:nFluxVMEC)=flmnc(NINT(xm(iMode)),NINT(xn(iMode)),0:nsin)
    END DO !iMode
    gmns=0.
    bmns=0.
  END IF

  !dont forget nfp on toroidal modes
  xn=REAL(nfp)*xn

  xm_nyq=xm; xn_nyq=xn

  phi= fsve*(-phiedge) /TWOPI
  chi = 0.

  CALL halfToFull(nFluxVMEC,hiota,iotaf)

  CALL halfToFull(nFluxVMEC,hpres,presf)
  presf=presf/(2.0e-07_wp*TWOPI) !/mu_0

  DEALLOCATE(fsve,hsve,hiota,hpres)
  DEALLOCATE(frmnc,frmns,fzmnc,fzmns)
  DEALLOCATE(flmnc,flmns)

END SUBROUTINE ReadNEMEC

ReadNEMEC Subroutine

Contents

Source Code

public subroutine ReadNEMEC(fileName, itype, ok)

Uses

Arguments

Calls

Called by

Source Code