synchrotron_cytran.f90

Go to the documentation of this file.

  SUBROUTINE synchrotron_cytran(EQUILIBRIUM, COREPROF, CORESOURCE, CODEPARAMS)

!---------------------------------------------------------------------------------------
!  This routine calculate sybchrotron radiation      
!
!     Source:       NTCC Library                        
!     Developers:   W. Houlberg, P. Strand              
!     Contacts:     par.strand@chalmers.se             
!----------------------------------------------------------------------------------------

    USE itm_types
    USE itm_constants
    USE euitm_schemas
    USE copy_structures
    USE allocate_deallocate
    USE cytran_mod      

!PIS    USE l3interps
      
    IMPLICIT NONE
 
! CPO derived types:
    TYPE (type_equilibrium), POINTER :: equilibrium(:)! input CPO equilibrium
    TYPE (type_coreprof),    POINTER :: coreprof(:)   ! input CPO main plasma
    TYPE (type_coresource),  POINTER :: coresource(:) ! output CPO with syncrhotron radiation
    TYPE (type_param)  :: codeparams                  ! CODE PARAMETERS

! local quantities  

! Dimensions: CPO related
    INTEGER :: nrho,                      & !number of radial points     (input, determined from COREPROF CPO)
               nnucl,                     & !number of ion species       (input, determined from COREPROF CPO)
               nion,                      & !number of ion species       (input, determined from COREPROF CPO)
               nimp,                      & !number of impurity species  (input, determined from COREPROF CPO)
               nneut                        !number of neutrals species  (input, determined from COREPROF CPO)

    INTEGER, ALLOCATABLE :: ntype(:),     & !number of impurity ionization states (input)     
                            ncomp(:),     & !max_number of distinct atoms enter the composition-"1" wich is neutral
                            nzimp(:)

    INTEGER, PARAMETER :: nocur = 1         !number of CPO ocurancies in the work flow

! Model parameters
    REAL(R8),ALLOCATABLE :: psync_ref(:), & ! Synchrotron term original formula [W/m**3] 
                            bmod_r(:),    & ! <|B|> Magnetic field strengt      [T]
                            ne_r(:),      & ! electron density                  [/m**3]    
                            te_r(:),      & ! electron temperature              [keV] 
                            dvol_r(:),    & ! Incremental volume                [m**3]
                            area_rm(:),   & ! Surface area                      [m**2]
                            psync_r(:)      ! Cytran estimate                   [W/m**3] 

    REAL(R8), SAVE :: f_self,             & ! x(o) to x(o) mode conversion      [-]
                      f_abs                 ! Wall absorption                   [-]   
   INTEGER, SAVE :: kmodel,               & !Internal model selectr = 0 CYTRAN  =1 Original calculation 
                    kflip,                & ! flip sign of original calculation (originally given as a loss) = 0 no, 1 yes
                    kverbosity              ! Output flag

   REAL(R8)        ::        ree, roo, roe, reo          ! REflection coefficients

! Grids 
    REAL(R8),ALLOCATABLE :: rho_r(:),     & ! Radial grid normalized toroidal flux [m]
                            rho_rm(:)       ! Surface area grid cell boundaries first point axis and nrho point at boundary
! Temporary placeholder
    REAL(R8),ALLOCATABLE :: temp(:)        

! Counters
    INTEGER                          :: i , return_status
! flags

      logical, save :: initial = .true.

! +++ Set dimensions:
    nrho = SIZE (coreprof(1)%rho_tor)
    CALL get_comp_dimensions(coreprof(1)%compositions, nnucl, nion, nimp, nzimp, nneut, ntype, ncomp)

! +++ Allocate output CPOs:
    CALL allocate_coresource_cpo(nocur, nrho, nnucl, nion, nimp, nzimp, nneut, ntype, ncomp, coresource)
    CALL copy_cpo(coreprof(1)%compositions, coresource(1)%compositions)

    coresource(1)%datainfo%cocos = 13
    coresource(1)%time = coreprof(1)%time
    coresource(1)%VALUES(1)%rho_tor = coreprof(1)%rho_tor
    coresource(1)%VALUES(1)%rho_tor_norm = coreprof(1)%rho_tor/coreprof(1)%rho_tor(nrho)
 
! +++ Set default variables

    if (initial ) THEN 
      kmodel = 0
      kflip = 0
      kverbosity = 1 
      f_self = 0.95_r8 ! Reference value is 0.95
      f_abs  = 0.80_r8 ! Reference values is in the range of 0.7 - 0.95
    
      CALL assign_code_parameters(codeparams, return_status)
      initial = .false.
    END IF 

    ree = f_self*(1.0_r8-f_abs)
    roo = ree
    roe = (1.0_r8-f_self)*(1.0_r8-f_abs)
    reo = roe

    ALLOCATE (area_rm(1:nrho), rho_rm(1:nrho), temp(1:nrho)) ! _rm grid sizes
    ALLOCATE (psync_ref(1:nrho-1), bmod_r(1:nrho-1), ne_r(1:nrho-1), te_r(1:nrho-1), & ! _r grid sizes
              dvol_r(1:nrho-1), rho_r(1:nrho-1), psync_r(1:nrho-1))

    rho_rm(:) = coreprof(1)%rho_tor(:)

    ! Map surface area
    call l3interp(equilibrium(1)%profiles_1d%surface,equilibrium(1)%profiles_1d%rho_tor, &
                   size(equilibrium(1)%profiles_1d%rho_tor), area_rm, rho_rm, nrho) 
    ! Map volume
    call l3interp(equilibrium(1)%profiles_1d%volume,equilibrium(1)%profiles_1d%rho_tor,  &
                   size(equilibrium(1)%profiles_1d%rho_tor), temp, rho_rm, nrho) 
    DO i = 1, nrho-1
      dvol_r(i) = temp(i+1)-temp(i)
      rho_r(i) = 0.5_r8*(rho_rm(i+1) + rho_rm(i)) 
    end do   

    ! Map Bmod
    call l3interp(equilibrium(1)%profiles_1d%gm5,equilibrium(1)%profiles_1d%rho_tor,     &
                   size(equilibrium(1)%profiles_1d%rho_tor), temp, rho_r, nrho-1) 
    bmod_r(1:nrho-1) =sqrt(temp(1:nrho-1)) 
      
   ! Calculate ne, Te on _rm grid
   DO i = 1, nrho-1
     ne_r(i) = 0.5_r8*(coreprof(1)%ne%value(i+1)+coreprof(1)%ne%value(i))
     te_r(i) = 0.5_r8*(coreprof(1)%te%value(i+1)+coreprof(1)%te%value(i))/1.0e3_r8
   END DO

! Ccalculate synchrotron radiation   
   SELECT CASE (kmodel)

      CASE (0) ! Cytran calculation

        psync_ref(:)  = -0.01_r8*6.2e-20_r8*coreprof(1)%toroid_field%B0**2.0*ne_r(:)*te_r(:)*1.0e3_r8
        CALL cytran(ree, reo, roo, roe, nrho-1, bmod_r, ne_r, te_r, area_rm, dvol_r, psync_r, k_cyt_res= 10 )
        ! Convert from kev/m**3/s to W/m**3
        psync_r(:) = psync_r(:) * itm_ev* 1.0e3_r8 

      CASE (1) !Original ETS formula

        psync_r(:)  = (-1.0_r8)**kflip*0.01_r8*6.2e-20_r8*coreprof(1)%toroid_field%B0**2.0*ne_r(:)*te_r(:)*1.0e3_r8

  END SELECT


! Diagnostic output

   IF (kverbosity /= 0 ) THEN 

     IF (kmodel == 1 ) THEN
       write(6,*) 'CYTRAN model: '
       write(6,*)
       write(6,'(A38,f7.4)') '           Wall reflection coefficient',f_abs           
       write(6,'(A38,f7.4)') '    x (o) mode to x(o) mode conversion', f_self
       write(6,*)
       write(6,'(A38,f7.4)') '   reflection coefficients: Roo = Ree =', ree 
       write(6,'(A38,f7.4)') '   reflection coefficients: Roe = Reo =', roe 

       write(6,*) ' i    rho     ne     Te   Bmod   dvol Psync (MW)'
       DO i = 1, nrho-1
         write(6,'(i3, 7f7.3)') i, rho_r(i), ne_r(i)/1.e20, te_r(i), bmod_r(i), dvol_r(i), &
                                   psync_r(i)/1.0e6_r8
       END DO

     ELSE IF (kmodel == 0 ) THEN
       write(6,*) 'Original ETS Formula: '
       write(6,*) 
       if (kflip == 0 ) THEN 
         write(6,*) ' Treat term as a sink (positive sign on source - need to be negative in source combiner)'
       ELSE IF (kflip == 1) THEN
        write(6,*) ' Treat term as a source (negative  sign on source -  eed to be positive in source combiner)'
     END IF
     write(6,*) 
     write(6,*) ' i    rho     ne     Te   Bmod   dvol  Psync  P_old (MW)'
     DO i = 1, nrho-1
         write(6,'(i3, 7f7.3)') i, rho_r(i), ne_r(i)/1.e20, te_r(i), bmod_r(i), dvol_r(i), &
                                   psync_r(i)/1.0e6_r8, psync_ref(i)/1.0e6_r8 
       END DO
     END IF
   ENDIF 

   ! Map back to output grid
   call l3interp(psync_r, rho_r, nrho-1,temp, rho_rm, nrho)
   temp(nrho) = temp(nrho-1) ! make a conservative extrapolation!
!Store in output CPO on -rm grid!
    coresource(1)%VALUES(1)%qe%exp(1:nrho) = temp(1:nrho)

!Add identifiers to COURESOURCE
    ALLOCATE(coresource%values(1)%sourceid%id(1))
    ALLOCATE(coresource(1)%values(1)%sourceid%description(1))
    coresource(1)%values(1)%sourceid%id = 'syncrotronradiation'
    coresource(1)%values(1)%sourceid%flag = 17
    coresource(1)%values(1)%sourceid%description = 'Source from syncrotron radiation cytran'


    IF(ALLOCATED(ntype)) DEALLOCATE (ntype)
    IF(ALLOCATED(ncomp)) DEALLOCATE (ncomp)
    IF(ALLOCATED(nzimp)) DEALLOCATE (nzimp)

    DEALLOCATE (area_rm, rho_rm, temp)
    DEALLOCATE (psync_ref, bmod_r, ne_r, te_r, dvol_r, rho_r, psync_r)


    RETURN

CONTAINS

    SUBROUTINE assign_code_parameters(codeparameters, return_status)
      USE itm_types
      USE euitm_schemas
      USE euitm_xml_parser
      USE xmllib_pathquery, ONLY : xmlget
      IMPLICIT NONE


      TYPE(type_param)       :: codeparameters
      INTEGER(ITM_I4)        :: return_status

      CHARACTER(len = 132)   :: prefix

      TYPE(tree)             :: parameter_list 
      INTEGER(ITM_I4)        :: ierror, nparm

      logical :: l_flip 


      return_status          = 0


!-- parse xml-string codeparameters%parameters
      CALL euitm_xml_parse(codeparameters, nparm, parameter_list)

      ! Set model - check to see if CYTRAN  unique values have been set
      prefix = 'parameters/MODEL/CYTRAN/f_self'
      CALL xmlget(parameter_list, prefix, f_self, error_flag= ierror)
      if (ierror /= 0 ) write(6,*) 'Path; ', trim(prefix), ' = ', f_self, ' errorflag =', ierror
      IF (ierror == 0 ) THEN  ! CYTRAN
          kmodel = 0 
          prefix = 'parameters/MODEL/CYTRAN/f_abs'
          CALL xmlget(parameter_list, prefix, f_abs, error_flag= ierror)
          IF (ierror /= 0 ) write(6,*) 'Path; ', trim(prefix), ' = ', f_abs, ' errorflag =', ierror

      ELSE 
          kmodel = 1
          prefix = 'parameters/MODEL/ETS_Formula/Flip_sign'
          CALL xmlget(parameter_list, prefix, l_flip, error_flag= ierror)
          IF (ierror /= 0) write(6,*) 'Path; ', trim(prefix), ' = ', l_flip, ' errorflag =', ierror
          if (l_flip) kflip =1
      END IF


      prefix = 'parameters/FLAGS/verbosity'

      CALL xmlget(parameter_list, prefix, kverbosity, error_flag= ierror)
      IF (ierror /= 0 ) write(6,*) 'Path; ', trim(prefix), ' = ', kverbosity, ' errorflag =', ierror

!-- destroy tree
      CALL destroy_xml_tree(parameter_list)
    END SUBROUTINE
  END SUBROUTINE synchrotron_cytran