cytran_mod.f90

Go to the documentation of this file.

MODULE cytran_mod
!-------------------------------------------------------------------------------
!CYTRAN-CYclotron radiation TRANsport
!
!CYTRAN_MOD is an F90 module that calculates cyclotran radiation transport in
!  toroidal plasmas
!
!References:
!
!  S.Tamor, SAIC report SAI-023-81-110LJ/LAPS-71 (1981)
!  S.Tamor, SAIC report SAI-023-81-189LJ/LAPS-72 (1981)
!  D.C.Baxter, S.Tamor, SAIC report SAI-023-81-215LJ/LAPS-74 (1981)
!  D.C.Baxter, S.Tamor, Nucl Technol/Fusion 3 (1983) 181
!  W.A.Houlberg, F90 free format 8/2004
!
!Contains PUBLIC routine:
!
!  CYTRAN              -cyclotron radiation transport
!
!Comments:
!
!  CYTRAN is an approximation to the cyclotron radiation transport process that
!    requires benchmarking against more complete transport analysis to improve
!    confidence in its results.
!
!  The magnetic geometry is simplified by assuming a circular cross section in a
!    plane of constant toroidal angle with an average magnetic field that is
!    only a function of minor radius.
!
!  Tamor benchmarked CYTRAN against the more complete cyclotron transport
!    analysis of the SNECTR code.
!
!  More recent benchmarking of CYTRAN against another cyclotron radiation code
!    has been carried out by:
!    F.Albajar, M.Bornatici, F.Engelmann, submitted to Nucl Fusion (2001)
!
!  Because the transport is expressed in terms of plasma surface area and
!    volumes, it can be applied to non-circular axisymmetric plasmas while
!    preserving the energy balance; the same is true if applied to
!    non-axisymmetric plasmas, but these must be viewed with more caution
!    because there have been no benchmark calculations for plasmas in which
!    field curvature in the toroidal direction may be more important.
!
!  In many places, physical and empirical constants are merged - this was done
!    in the original coding - and have not been unraveled in the effort to
!    convert to a more modern and cleaner code.
!
!  The modernization of the code structure into an F90 module takes advantage of
!    some of the more attractive features of F90:
!    -use of KIND for precision declarations
!    -optional arguments for I/O
!    -generic names for all intrinsic functions
!    -compilation using either free or fixed form
!    -no common blocks or other deprecated Fortran features
!    -dynamic and automatic alocation of variables
!    -array syntax for vector operations
!-------------------------------------------------------------------------------
USE itm_types
USE itm_constants

IMPLICIT NONE

!-------------------------------------------------------------------------------
! Private procedures
!-------------------------------------------------------------------------------
PRIVATE :: &
  cytran_opacity         !cyclotron opacity coefficients
                         !  called from CYTRAN

!-------------------------------------------------------------------------------
! Private data
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
! Public procedures
!-------------------------------------------------------------------------------
CONTAINS

SUBROUTINE cytran(ree,reo,roo,roe,nr_r,bmod_r,den_r,te_r,area_rm,dvol_r, &
                  psync_r, &
                  k_cyt_res)
!-------------------------------------------------------------------------------
!CYTRAN estimates cyclotron radiation energy transport in toroidal plasmas
!
!References:
!  S.Tamor, SAIC report SAI-023-81-110LJ/LAPS-71 (1981)
!  S.Tamor, SAIC report SAI-023-81-189LJ/LAPS-72 (1981)
!  D.C.Baxter, S.Tamor, SAIC report SAI-023-81-215LJ/LAPS-74 (1981)
!  D.C.Baxter, S.Tamor, Nucl Technol/Fusion 3 (1983) 181
!  W.A.Houlberg, F90 free format 8/2004
!-------------------------------------------------------------------------------

!Declaration of input variables
INTEGER, INTENT(IN) :: &
  nr_r                   !no. of radial plasma cells [-]

REAL(KIND=R8), INTENT(IN) :: &
  ree,                 & !refl of X mode from incident X mode [0-1]
  reo,                 & !refl of X mode from incident O mode [0-1]
  roe,                 & !refl of O mode from incident X mode [0-1]
  roo                    !refl of O mode from incident O mode [0-1]

REAL(KIND=R8), INTENT(IN) :: &
  bmod_r(:),           & !<|B|> [T]
  den_r(:),            & !electron density in cell i [/m**3]
  dvol_r(:),           & !cell volume [m**3]
  te_r(:),             & !electron temperature in cell i [keV]
  area_rm(:)             !surface area at inner boundary of cell [m**2]
                         !plasma outer surface area is at nr_r+1


!Declaration of optional input variables
INTEGER, INTENT(IN), OPTIONAL :: &
  k_cyt_res              !frequency resolution level [-]
                         !=1 delta(omega)=default of 100 frequncy intervals
                         !10 recommended to keep graininess < a couple %


!Declaration of output variables
REAL(KIND=R8), INTENT(OUT) :: &
  psync_r(:)             !net power source/loss in cell i [keV/m**3/s]

!-------------------------------------------------------------------------------
!Declaration of local variables
INTEGER :: &
  i,j,k_res,nze,nzo

REAL(KIND=R8) :: &
  areae,areao,blackb,delw,dwght,fwght,emass,omass,extray,ordray, &
  qee,qeo,qoe,qoo,srce,srco,taucrt,taue,tauo,testi,tse,tso,wfreq, &
  wght,wmax,wmin,xe,xo

REAL(KIND=R8) :: &
  alphae(1:nr_r),alphao(1:nr_r),delta(1:nr_r),temp(1:nr_r)

!-------------------------------------------------------------------------------
!Initialization
!-------------------------------------------------------------------------------
!Frequency resolution
k_res=1

IF(present(k_cyt_res)) THEN

  IF(k_cyt_res > 1) k_res=k_cyt_res

ENDIF

!Minimum and maximum frequencies and interval
wmin=5.0e11_r8*minval(bmod_r(1:nr_r))
wmax=2.4e12_r8*maxval(bmod_r(1:nr_r))
delw=0.08_r8*(wmax-wmin)/k_res

!Other
dwght=1
fwght=3
taucrt=0.6_r8
testi=0
extray=1.0e36_r8

psync_r(1:nr_r)=0
temp(1:nr_r)=te_r(1:nr_r)*dvol_r(1:nr_r)

DO i=1,nr_r-1 !Over radial zones

  delta(i)=2*dvol_r(i)/(area_rm(i)+area_rm(i+1))

ENDDO !Over radial zones

delta(nr_r)=dvol_r(nr_r)/area_rm(nr_r)

DO j=1,100*k_res !Over frequency

  IF(extray < 0.03_r8*testi) EXIT

  wfreq=wmin+delw*REAL(j-1,r8)
  dwght=-dwght
  wght=fwght+dwght
  IF(j == 1) wght=1
  wght=wght*delw/3

  !Blackbody radiation factor (/m**2/s)
  blackb=4.48e-20_r8*wfreq*wfreq

  !Get opacity factors
  CALL cytran_opacity(nr_r,te_r,bmod_r,wfreq,alphae,alphao)
  alphae(1:nr_r)=alphae(1:nr_r)*6.03e-17_r8*den_r(1:nr_r)/bmod_r(1:nr_r)
  alphao(1:nr_r)=alphao(1:nr_r)*6.03e-17_r8*den_r(1:nr_r)/bmod_r(1:nr_r)

  taue=0
  tauo=0
  nze=nr_r
  nzo=nr_r

  DO i=nr_r,1,-1 !Over radial zones

    taue=taue+alphae(i)*delta(i)
    IF(taue < taucrt) nze=i-1
    tauo=tauo+alphao(i)*delta(i)
    IF(tauo < taucrt) nzo=i-1

  ENDDO !Over radial zones 

  !Extraordinary mode
  IF(nze > 0) THEN

    areae=area_rm(nze+1)
    srce=areae*te_r(nze)
    emass=sum(den_r(1:nr_r)*dvol_r(1:nr_r))

  ELSE

    areae=0
    srce=0
    emass=0

  ENDIF

  !Ordinary mode
  IF(nzo > 0) THEN

    areao=area_rm(nzo+1)
    srco=areao*te_r(nzo)
    omass=sum(den_r(1:nr_r)*dvol_r(1:nr_r))

  ELSE

    areao=0
    srco=0
    omass=0

  ENDIF

!Extraordinary mode
  IF(nr_r > nze) THEN

    xe=4*sum(dvol_r(nze+1:nr_r)*alphae(nze+1:nr_r))
    srce=srce+4*sum(temp(nze+1:nr_r)*alphae(nze+1:nr_r))

  ELSE

    xe=0

  ENDIF

!Ordinary mode
  IF(nr_r > nzo) THEN

    xo=4*sum(dvol_r(nzo+1:nr_r)*alphao(nzo+1:nr_r))
    srco=srco+4*sum(temp(nzo+1:nr_r)*alphao(nzo+1:nr_r))

  ELSE

    xo=0

  ENDIF

  qee=areao+area_rm(nr_r+1)*(1.0_r8-roo)+xo
  qoo=areae+area_rm(nr_r+1)*(1.0_r8-ree)+xe
  qeo=reo*area_rm(nr_r+1)
  qoe=roe*area_rm(nr_r+1)
  extray=(qee*srce+qeo*srco)*blackb/(qee*qoo-qeo*qoe)
  ordray=(qoe*srce+qoo*srco)*blackb/(qee*qoo-qeo*qoe)
  testi=max(testi,extray)

  DO i=1,nr_r !Over radial zones

    !Extraordinary mode
    IF(i <= nze) THEN

      tse=itm_pi*areae*(blackb*te_r(nze)-extray)*(den_r(i)/emass)*wght

    ELSE

      tse=4*itm_pi*alphae(i)*(blackb*te_r(i)-extray)*wght

    ENDIF

    psync_r(i)=psync_r(i)-tse

    !Ordinary mode
    IF(i <= nzo) THEN

      tso=itm_pi*areao*(blackb*te_r(nzo)-ordray)*(den_r(i)/omass)*wght

    ELSE

      tso=4*itm_pi*alphao(i)*(blackb*te_r(i)-ordray)*wght

    ENDIF

    psync_r(i)=psync_r(i)-tso

  ENDDO !Over radial zones

ENDDO !Over frequency

END SUBROUTINE cytran

SUBROUTINE cytran_opacity(nr_r,te_r,bmod_r,wfreq,alphae,alphao)
!-------------------------------------------------------------------------------
!CYTRAN_OPACITY determines the opacity coefficients for cyclotron radiation
!
!References:
!  S.Tamor, SAIC report SAI-023-81-110LJ/LAPS-71 (1981)
!  S.Tamor, SAIC report SAI-023-81-189LJ/LAPS-72 (1981)
!  D.C.Baxter, S.Tamor, SAIC report SAI-023-81-215LJ/LAPS-74 (1981)
!  D.C.Baxter, S.Tamor, Nucl Technol/Fusion 3 (1983) 181
!  W.A.Houlberg, F90 free format 8/2004
!-------------------------------------------------------------------------------

!Declaration of input variables
INTEGER, INTENT(IN) :: &
  nr_r                   !no. of radial cells [-]

REAL(KIND=R8), INTENT(IN) :: &
  wfreq                  !wave angular frequency [radians/s]

REAL(KIND=R8), INTENT(IN) :: &
  te_r(:),             & !electron temperature [keV]
  bmod_r(:)

!Declaration of output variables
REAL(KIND=R8), INTENT(OUT) :: &
  alphae(:),           & !extra-ordinary mode opacity [-]
  alphao(:)              !ordinary mode opacity [-]

!-------------------------------------------------------------------------------
!Declaration of local variables
INTEGER :: &
  i

REAL(KIND=R8) :: &
  arg,tden,wloc

DO i=1,nr_r !Over zones

  !wloc=wfreq/electron cyclotron frequency
  wloc=wfreq/(1.7588e+11*bmod_r(i))
  tden=1.0_r8/(4.0+te_r(i)+25.0/te_r(i))
  arg=max(0.0_r8,0.045_r8+(wloc-2.0_r8)*tden)
  alphae(i)=10.0**(1.45-7.8*sqrt(arg))/(wloc**2)
  arg=max(0.0_r8,0.180_r8+(wloc-1.0_r8)*tden)
  alphao(i)=10.0_r8**(2.45_r8-8.58_r8*sqrt(arg))/(wloc**2)

ENDDO !Over zones

END SUBROUTINE cytran_opacity

END MODULE cytran_mod