analytics1.f90

Go to the documentation of this file.

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
MODULE analytics1

CONTAINS




! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
  SUBROUTINE analytical_plasma (                                            &
       time,coreprof_in,                                                    &
       equilibrium,coreprof_analytic,coretransp,coresource,coreimpur, code_parameters)       


!-------------------------------------------------------!
!     This routine manufactures the solution for the    !
!     set of transport equations describing the main    !
!     plasma. It provides sources and transport         !
!     coefficients for all transport equations.         !
!-------------------------------------------------------!
!     Source:       ---                                 !
!     Developers:   R.Stankiewicz, D.Kalupin            !
!     Kontacts:     Roman.Stankiewich@gmail.com         !
!                   D.Kalupin@fz-juelich.de             !
!                                                       !
!     Comments:     ---                                 !
!                                                       !
!-------------------------------------------------------!


    USE itm_types
    USE euitm_schemas
    USE ets_plasma
    USE itm_test_routines
    USE allocate_deallocate
    USE copy_cpo_ets
    USE is_set_structures
    USE copy_structures
    use deallocate_structures
!!    use euitm_printcpo

    IMPLICIT NONE


    INTEGER, SAVE                    :: shot_no             !shot mumber
    INTEGER, SAVE                    :: run_no              !run mumber
    
    INTEGER, SAVE                    :: nrho                !number of radial points    (input)
    INTEGER, SAVE                    :: nnucl               !number of ion species      (input)
    INTEGER, SAVE                    :: nneut               !number of ion species      (input)
    INTEGER, SAVE                    :: nion                !number of ion species      (input)
    INTEGER, SAVE                    :: nimp                !number of impurity species (input)
    INTEGER, ALLOCATABLE, SAVE       :: nzimp(:)            !number of charge states for each impurity (input)
    INTEGER, ALLOCATABLE, SAVE       :: ncomp(:)            !number of charge states for each impurity (input)
    INTEGER, ALLOCATABLE, SAVE       :: ntype(:)            !number of charge states for each impurity (input)
    INTEGER, SAVE                    :: npoints             !number of charge states for each impurity (input)
    INTEGER, SAVE                    :: ntime               !number of time points      (input)
    INTEGER,               PARAMETER :: nocur = 1           !number of CPO ocurancies in the work flow
    
    INTEGER                          :: irho                !current radial knot
    INTEGER                          :: iion                !current ion type
    INTEGER                          :: itime               !current time step
    
    INTEGER, SAVE                    :: nsol                !Number of analytical example
    INTEGER, SAVE                    :: solver_type         !representation of transport equations 
    !1-"standard"; 2-"integral"(default)
    INTEGER, SAVE                    :: sigma_source         !origin of Plasma electrical conductivity
    !0: plasma collisions;1: transport module;2: source module
    
    REAL (R8), SAVE                  :: convrec             !required convergency 
    REAL (R8)                        :: time                !Time
    REAL (R8), SAVE                  :: tau                 !time step, and mixing coefficient
    REAL (R8), SAVE                  :: amix                !mixing factor
    INTEGER                          :: iter                !iteration index
    INTEGER, PARAMETER               :: maxiter=1000        !maximum number of convergence iterations
    
    INTEGER, SAVE                    :: psi_bnd_type        !Type of boundary conditions current
    INTEGER, SAVE                    :: ni_bnd_type         !Type of boundary conditions ion density 
    INTEGER, SAVE                    :: ti_bnd_type         !Type of boundary conditions ion temperature
    INTEGER, SAVE                    :: te_bnd_type         !Type of boundary conditions electron temperature
    INTEGER, SAVE                    :: vtor_bnd_type       !Type of boundary conditions toroidal rotation
    
    TYPE (run_control)               :: control             !contains all parameters required by run

    INTEGER                          :: ifail


! +++ CPO derived types:
    TYPE (type_equilibrium),POINTER        :: equilibrium(:)       !input CPO with geometry
    TYPE (type_coreprof),   POINTER        :: coreprof_in(:)       !CPO with internal ETS parameters profiles from previous time step
    TYPE (type_coreprof),   POINTER        :: coreprof_analytic(:) !CPO with profiles
    TYPE (type_coretransp), POINTER        :: coretransp(:)        !CPO with transport coefficients
    TYPE (type_coresource), POINTER        :: coresource(:)        !CPO with sources
    TYPE (type_coreimpur),  POINTER        :: coreimpur(:)         !CPO with impurities
    TYPE (type_param) :: code_parameters


! +++ Local derived types:
    TYPE (magnetic_geometry)               :: geometry1         !contains all geometry quantities
    TYPE (plasma_profiles)                 :: profiles1         !contains profiles of plasma parameters
    TYPE (transport_coefficients)          :: transport1        !contains profiles of trasport coefficients
    TYPE (sources_and_sinks)               :: sources1          !contains profiles of sources
    TYPE (global_param)                    :: global            !contains global

    CHARACTER (len=32)               :: database_format

    LOGICAL, SAVE :: first=.true.

! +++ Set up dimensions:

    WRITE(6,*) time
!!    call is_set_cpo(coreprof_in(1),'coreprof_in')
!!!    flush(6)

    IF(first) THEN
       IF (.NOT. ASSOCIATED(code_parameters%parameters)) THEN
          WRITE(6, *) 'ERROR: code parameters not associated!'
          stop
       END IF
       CALL process_xml(solver_type,sigma_source,tau,amix,convrec,nrho,nion,nimp,nzimp,ntime,nsol,      &
            psi_bnd_type,ni_bnd_type,ti_bnd_type,te_bnd_type,vtor_bnd_type,                &
            shot_no, run_no, code_parameters, database_format)
       first=.false.
    ENDIF
    
    npoints = 101
    nnucl = nion + nimp
    nneut = nion + nimp
    if(.not. allocated(ntype)) ALLOCATE (ntype(nneut)) 
    if(.not. allocated(ncomp)) ALLOCATE (ncomp(nneut))
    ntype = 1
    ncomp = 1

    CALL allocate_equilibrium_cpo(nocur,  nrho,  nrho,  101,   npoints,                           equilibrium      )         
        
    CALL allocate_coreprof_cpo(nocur,  nrho,  nnucl, nion,  nimp,  nzimp, nneut, ntype, ncomp, coreprof_analytic    )
    CALL allocate_coretransp_cpo(nocur,  nrho,  nnucl, nion,  nimp,  nzimp, nneut, ntype, ncomp, coretransp  )        
    CALL allocate_coresource_cpo(nocur,  nrho,  nnucl, nion,  nimp,  nzimp, nneut, ntype, ncomp, coresource  )

    if(nimp .GT. 0) &
    CALL allocate_coreimpur_cpo(nocur,  nrho,  nnucl, nion,  nimp,  nzimp, nneut, ntype, ncomp, coreimpur   )        


    CALL copy_codeparam(coreprof_in(1)%codeparam, coreprof_analytic(1)%codeparam)
    call deallocate_cpo(coreprof_analytic(1)%COMPOSITIONS)
    CALL copy_cpo(coreprof_in(1)%COMPOSITIONS,        coreprof_analytic(1)%COMPOSITIONS)
    call deallocate_cpo(coretransp(1)%COMPOSITIONS)
    CALL copy_cpo(coreprof_in(1)%COMPOSITIONS,        coretransp(1)%COMPOSITIONS)
    call deallocate_cpo(coresource(1)%COMPOSITIONS)
    CALL copy_cpo(coreprof_in(1)%COMPOSITIONS,        coresource(1)%COMPOSITIONS)

! +++ Set up parameters for the test work flow:
!    CALL perfon ('setcpo')
    CALL set_cpo(nrho,nion,nimp,nzimp,ntime,nsol,                               &
         psi_bnd_type,ni_bnd_type,ti_bnd_type,te_bnd_type,vtor_bnd_type,          &
         equilibrium, coreprof_in, coretransp)
!    CALL perfoff
    

! +++ Allocation of local derived types:
    CALL allocate_magnetic_geometry(nrho,       geometry1,  ifail)
    CALL allocate_plasma_profiles(nrho,nion,  profiles1,  ifail)
    CALL allocate_transport_coefficients(nrho,nion,  transport1, ifail)
    CALL allocate_sources_and_sinks(nrho,nion,  sources1,   ifail)
    CALL allocate_global_param(            global,     ifail)



! +++ Copy CPOs in local derived types:
    CALL convert_cpo_to_local_types(coreprof_in,  profiles1, transport1)



! +++ Call MAIN_PLASMA with internal ETS derived types:
    CALL analytics_full(time, geometry1, profiles1, transport1, sources1) 



! +++ Copy local derived types in CPOs:
    CALL convert_local_to_cpo_types(geometry1,   profiles1, transport1, sources1,    &
         equilibrium, coreprof_analytic,  coretransp, coresource)


! +++ Deallocation of ETS derived types:
    CALL deallocate_magnetic_geometry(geometry1,  ifail)
    CALL deallocate_plasma_profiles(profiles1,  ifail)
    CALL deallocate_transport_coefficients(transport1, ifail)
    CALL deallocate_sources_and_sinks(sources1,   ifail)


    WRITE(6,*) 'EQUILIBRIUM'
    IF(.NOT.ASSOCIATED(equilibrium(1)%codeparam%codename)) THEN
       ALLOCATE(equilibrium(1)%codeparam%codename(1))   ! For a string of 132 characters max.
       equilibrium(1)%codeparam%codename(1)   = 'ANALYTIC_SOLVER_TEST'
    ENDIF
    IF(.NOT.ASSOCIATED(equilibrium(1)%codeparam%codeversion)) THEN
       ALLOCATE(equilibrium(1)%codeparam%codeversion(1))   ! For a string of 132 characters max.
       equilibrium(1)%codeparam%codeversion(1)   = '?'
    ENDIF
    IF(.NOT.ASSOCIATED(equilibrium(1)%codeparam%parameters)) THEN
       ALLOCATE(equilibrium(1)%codeparam%parameters(1))   ! For a string of 132 characters max.
       equilibrium(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
    ENDIF
    IF(.NOT.ASSOCIATED(equilibrium(1)%codeparam%output_diag)) THEN
       ALLOCATE(equilibrium(1)%codeparam%output_diag(1))   ! For a string of 132 characters max.
       equilibrium(1)%codeparam%output_diag(1) = 'my_output_diag'
    ENDIF
    equilibrium(1)%codeparam%output_flag = 0   ! Integer output flag, 0 means the run was successful and can be used in the rest of the workflow, <0 means failure
    equilibrium(1)%time=time
!!    call printcpo(EQUILIBRIUM(1))
    WRITE(6,*) 'COREPROF_ANALYTIC'
    IF(.NOT.ASSOCIATED(coreprof_analytic(1)%codeparam%codename)) THEN
       ALLOCATE(coreprof_analytic(1)%codeparam%codename(1))   ! For a string of 132 characters max.
       coreprof_analytic(1)%codeparam%codename(1)   = 'ANALYTIC_SOLVER_TEST'
    ENDIF
    IF(.NOT.ASSOCIATED(coreprof_analytic(1)%codeparam%codeversion)) THEN
       ALLOCATE(coreprof_analytic(1)%codeparam%codeversion(1))   ! For a string of 132 characters max.
       coreprof_analytic(1)%codeparam%codeversion(1)   = '?'
    ENDIF
    IF(.NOT.ASSOCIATED(coreprof_analytic(1)%codeparam%parameters)) THEN
       ALLOCATE(coreprof_analytic(1)%codeparam%parameters(1))   ! For a string of 132 characters max.
       coreprof_analytic(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
    ENDIF
    IF(.NOT.ASSOCIATED(coreprof_analytic(1)%codeparam%output_diag)) THEN
       ALLOCATE(coreprof_analytic(1)%codeparam%output_diag(1))   ! For a string of 132 characters max.
       coreprof_analytic(1)%codeparam%output_diag(1) = 'my_output_diag'
    ENDIF
    coreprof_analytic(1)%codeparam%output_flag = 0   ! Integer output flag, 0 means the run was successful and can be used in the rest of the workflow, <0 means failure
    coreprof_analytic(1)%time=time
!!    call printcpo(COREPROF_ANALYTIC(1))
    WRITE(6,*) 'CORETRANSP'
    IF(.NOT.ASSOCIATED(coretransp(1)%codeparam%codename)) THEN
       ALLOCATE(coretransp(1)%codeparam%codename(1))   ! For a string of 132 characters max.
       coretransp(1)%codeparam%codename(1)   = 'ANALYTIC_SOLVER_TEST'
    ENDIF
    IF(.NOT.ASSOCIATED(coretransp(1)%codeparam%codeversion)) THEN
       ALLOCATE(coretransp(1)%codeparam%codeversion(1))   ! For a string of 132 characters max.
       coretransp(1)%codeparam%codeversion(1)   = '?'
    ENDIF
    IF(.NOT.ASSOCIATED(coretransp(1)%codeparam%parameters)) THEN
       ALLOCATE(coretransp(1)%codeparam%parameters(1))   ! For a string of 132 characters max.
       coretransp(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
    ENDIF
    IF(.NOT.ASSOCIATED(coretransp(1)%codeparam%output_diag)) THEN
       ALLOCATE(coretransp(1)%codeparam%output_diag(1))   ! For a string of 132 characters max.
       coretransp(1)%codeparam%output_diag(1) = 'my_output_diag'
    ENDIF
    coretransp(1)%codeparam%output_flag = 0   ! Integer output flag, 0 means the run was successful and can be used in the rest of the workflow, <0 means failure
    coretransp(1)%time=time
!!    call printcpo(CORETRANSP(1))
    WRITE(6,*) 'CORESOURCE'
    IF(.NOT.ASSOCIATED(coresource(1)%codeparam%codename)) THEN
       ALLOCATE(coresource(1)%codeparam%codename(1))   ! For a string of 132 characters max.
       coresource(1)%codeparam%codename(1)   = 'ANALYTIC_SOLVER_TEST'
    ENDIF
    IF(.NOT.ASSOCIATED(coresource(1)%codeparam%codeversion)) THEN
       ALLOCATE(coresource(1)%codeparam%codeversion(1))   ! For a string of 132 characters max.
       coresource(1)%codeparam%codeversion(1)   = '?'
    ENDIF
    IF(.NOT.ASSOCIATED(coresource(1)%codeparam%parameters)) THEN
       ALLOCATE(coresource(1)%codeparam%parameters(1))   ! For a string of 132 characters max.
       coresource(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
    ENDIF
    IF(.NOT.ASSOCIATED(coresource(1)%codeparam%output_diag)) THEN
       ALLOCATE(coresource(1)%codeparam%output_diag(1))   ! For a string of 132 characters max.
       coresource(1)%codeparam%output_diag(1) = 'my_output_diag'
    ENDIF
    coresource(1)%codeparam%output_flag = 0   ! Integer output flag, 0 means the run was successful and can be used in the rest of the workflow, <0 means failure
    coresource(1)%time=time
!!    call printcpo(CORESOURCE(1))
    if(nimp.gt.0) then
       WRITE(6,*) 'COREIMPUR'
       IF(.NOT.ASSOCIATED(coreimpur(1)%codeparam%codename)) THEN
          ALLOCATE(coreimpur(1)%codeparam%codename(1))   ! For a string of 132 characters max.
          coreimpur(1)%codeparam%codename(1)   = 'ANALYTIC_SOLVER_TEST'
       ENDIF
       IF(.NOT.ASSOCIATED(coreimpur(1)%codeparam%codeversion)) THEN
          ALLOCATE(coreimpur(1)%codeparam%codeversion(1))   ! For a string of 132 characters max.
          coreimpur(1)%codeparam%codeversion(1)   = '?'
       ENDIF
       IF(.NOT.ASSOCIATED(coreimpur(1)%codeparam%parameters)) THEN
          ALLOCATE(coreimpur(1)%codeparam%parameters(1))   ! For a string of 132 characters max.
          coreimpur(1)%codeparam%parameters(1) = 'my_code_specific_parameters'
       ENDIF
       IF(.NOT.ASSOCIATED(coreimpur(1)%codeparam%output_diag)) THEN
          ALLOCATE(coreimpur(1)%codeparam%output_diag(1))   ! For a string of 132 characters max.
          coreimpur(1)%codeparam%output_diag(1) = 'my_output_diag'
       ENDIF
       coreimpur(1)%codeparam%output_flag = 0   ! Integer output flag, 0 means the run was successful and can be used in the rest of the workflow, <0 means failure
       coreimpur(1)%time=time
!!    call printcpo(COREIMPUR(1))
    endif
    WRITE(6,*) 'Finished in ANALYTICAL_PLASMA'
!!!    flush(6)

    RETURN



  END SUBROUTINE analytical_plasma

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
  SUBROUTINE convert_cpo_to_local_types  (COREPROF,  PROFILES1, TRANSPORT1)

!-------------------------------------------------------!
!     This routine converts CPOs into the local         !
!     derived types.                                    !
!-------------------------------------------------------!
!     Source:       ---                                 !
!     Developers:   D.Kalupin                           !
!     Kontacts:     D.Kalupin@fz-juelich.de             !
!                                                       !
!     Comments:     ---                                 !
!                                                       !
!-------------------------------------------------------!


    USE euitm_schemas

    USE ets_plasma

    IMPLICIT NONE

    INTEGER :: iion, nion

! +++ CPO derived types:
    TYPE (type_coreprof),          POINTER     :: coreprof(:)      !input CPO with plasma profiles
! +++ Local derived types:
    TYPE (plasma_profiles),        INTENT(INOUT) :: profiles1        !contains profiles of plasma parameters
    TYPE (transport_coefficients), INTENT(INOUT) :: transport1       !contains profiles of trasport coefficients



! +++ Convert profiles:
    nion = size(coreprof(1)%compositions%ions)
    DO iion =1, nion
       profiles1%ZION(iion)   = coreprof(1)%compositions%ions(iion)%zion 
       profiles1%ZION2(iion)  = coreprof(1)%compositions%ions(iion)%zion**2 
       profiles1%MION(iion)   = coreprof(1)%compositions%nuclei(coreprof(1)%compositions%ions(iion)%nucindex)%amn
    END DO

    profiles1%PSI_BND_TYPE    = coreprof(1)%psi%boundary%type

    profiles1%NI_BND_TYPE     = coreprof(1)%ni%boundary%type

    profiles1%TI_BND_TYPE     = coreprof(1)%ti%boundary%type

    profiles1%TE_BND_TYPE     = coreprof(1)%te%boundary%type

    profiles1%VTOR_BND_TYPE   = coreprof(1)%vtor%boundary%type



! +++ Convert transport:
    transport1%C1(1)          = 0.0e0_r8                               
    transport1%C1(2)          = 1.5e0_r8                               
    transport1%C1(3)          = 2.5e0_r8                               




    RETURN


  END SUBROUTINE convert_cpo_to_local_types


! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  




! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
  SUBROUTINE convert_local_to_cpo_types  (GEOMETRY1,   PROFILES1, TRANSPORT1, SOURCES1,    &
       equilibrium, coreprof,  coretransp, coresource)

!-------------------------------------------------------!
!     This routine converts local into the CPOs         !
!     derived types.                                    !
!-------------------------------------------------------!
!     Source:       ---                                 !
!     Developers:   D.Kalupin                           !
!     Kontacts:     D.Kalupin@fz-juelich.de             !
!                                                       !
!     Comments:     ---                                 !
!                                                       !
!-------------------------------------------------------!


    USE euitm_schemas

    USE ets_plasma

    USE itm_constants

    IMPLICIT NONE

! +++ Local derived types:
    TYPE (magnetic_geometry),      INTENT(IN)  :: geometry1        !contains all geometry quantities
    TYPE (plasma_profiles),        INTENT(IN)  :: profiles1        !contains profiles of plasma parameters
    TYPE (transport_coefficients), INTENT(IN)  :: transport1       !contains profiles of trasport coefficients
    TYPE (sources_and_sinks),      INTENT(IN)  :: sources1         !contains profiles of sources

! +++ CPO derived types:
    TYPE (type_equilibrium),       POINTER     :: equilibrium(:)   !output CPO with geometry
    TYPE (type_coreprof),          POINTER     :: coreprof(:)      !output CPO with plasma profiles
    TYPE (type_coretransp),        POINTER     :: coretransp(:)    !output CPO with transport coefficients
    TYPE (type_coresource),        POINTER     :: coresource(:)    !output CPO with sources



! +++ Convert geometry:
    equilibrium(1)%profiles_1d%rho_tor     = geometry1%RHO
!!!DPC-EQ-4.09a
    equilibrium(1)%profiles_1d%vprime      = geometry1%VPR
    equilibrium(1)%profiles_1d%gm3         = geometry1%G1
    equilibrium(1)%profiles_1d%gm8         = geometry1%G2     ! added DPC 2012-04-27
    equilibrium(1)%profiles_1d%gm2         = geometry1%G3
    equilibrium(1)%profiles_1d%F_dia       = geometry1%FDIA

    coreprof(1)%rho_tor                    = geometry1%RHO
    coreprof(1)%toroid_field%r0            = geometry1%RGEO 
    coreprof(1)%toroid_field%b0            = geometry1%BGEO 

    coretransp(1)%values(1)%rho_tor                  = geometry1%RHO
    coresource(1)%values(1)%rho_tor                  = geometry1%RHO



! +++ Convert profiles:
!    allocate(COREPROF(1)%compositions%ions(NION))
!    DO IION =1, NION
!       COREPROF(1)%compositions%ions(iion)%zion           = PROFILES1%ZION(IION)  
!       COREPROF(1)%composition%amn            = PROFILES1%MION 
!    END DO

    coreprof(1)%psi%value                  = profiles1%PSI 
    coreprof(1)%profiles1d%q%value         = profiles1%QSF
    coreprof(1)%psi%boundary%value         = profiles1%PSI_BND
    coreprof(1)%psi%boundary%type          = profiles1%PSI_BND_TYPE

    coreprof(1)%ni%value                   = profiles1%NI
    coreprof(1)%ni%boundary%value          = profiles1%NI_BND
    coreprof(1)%ni%boundary%type           = profiles1%NI_BND_TYPE
!DPC
!    COREPROF(1)%ni%transp_coef%diff        = sum(TRANSPORT1%DIFF_NI,dim=3)
!    COREPROF(1)%ni%transp_coef%vconv       = sum(TRANSPORT1%VCONV_NI,dim=3)
!DPC.end

    coreprof(1)%ne%value                   = profiles1%NE
!DPC
!    COREPROF(1)%ne%transp_coef%diff        = TRANSPORT1%DIFF_NI
!    COREPROF(1)%ne%transp_coef%vconv       = TRANSPORT1%VCONV_NI
!DPC.end

    coreprof(1)%ti%value                   = profiles1%TI
    coreprof(1)%ti%boundary%value          = profiles1%TI_BND
    coreprof(1)%ti%boundary%type           = profiles1%TI_BND_TYPE
!DPC
    coreprof(1)%ti%transp_coef%diff        = transport1%DIFF_TI
    coreprof(1)%ti%transp_coef%vconv       = transport1%VCONV_TI
!DPC.end

    coreprof(1)%te%value                   = profiles1%TE
    coreprof(1)%te%boundary%value          = profiles1%TE_BND
    coreprof(1)%te%boundary%type           = profiles1%TE_BND_TYPE
!DPC
    coreprof(1)%te%transp_coef%diff        = transport1%DIFF_TE
    coreprof(1)%te%transp_coef%vconv       = transport1%VCONV_TE
!DPC.end

    coreprof(1)%vtor%value                 = profiles1%VTOR
    coreprof(1)%vtor%boundary%value        = profiles1%VTOR_BND
    coreprof(1)%vtor%boundary%type         = profiles1%VTOR_BND_TYPE
!DPC
    coreprof(1)%vtor%transp_coef%diff      = transport1%DIFF_VTOR
    coreprof(1)%vtor%transp_coef%vconv     = transport1%VCONV_VTOR
!DPC.end

!DPC set fluxes to 0 (otherwise we write garbage)
    coreprof(1)%ni%flux%flux_dv            = 0.0_r8
    coreprof(1)%ne%flux%flux_dv            = 0.0_r8
    coreprof(1)%ti%flux%flux_dv            = 0.0_r8
    coreprof(1)%te%flux%flux_dv            = 0.0_r8
    coreprof(1)%vtor%flux%flux_dv          = 0.0_r8
!DPC.end

! +++ Convert transport:
    coretransp(1)%values(1)%sigma                    = transport1%SIGMA

    coretransp(1)%values(1)%te_transp%diff_eff       = transport1%DIFF_TE
    coretransp(1)%values(1)%te_transp%vconv_eff      = transport1%VCONV_TE

    coretransp(1)%values(1)%ni_transp%diff_eff       = transport1%DIFF_NI 
    coretransp(1)%values(1)%ni_transp%vconv_eff      = transport1%VCONV_NI

    coretransp(1)%values(1)%ti_transp%diff_eff       = transport1%DIFF_TI 
    coretransp(1)%values(1)%ti_transp%vconv_eff      = transport1%VCONV_TI 

    coretransp(1)%values(1)%vtor_transp%diff_eff     = transport1%DIFF_VTOR 
    coretransp(1)%values(1)%vtor_transp%vconv_eff    = transport1%VCONV_VTOR



! +++ Convert sources:
    coresource(1)%VALUES(1)%j                        = sources1%CURR_EXP

    coresource(1)%VALUES(1)%qe%exp                   = sources1%QE_EXP * itm_ev    !! DPC 2009-07-02
    coresource(1)%VALUES(1)%qe%imp                   = sources1%QE_IMP

    coresource(1)%VALUES(1)%si%exp                   = sources1%SI_EXP
    coresource(1)%VALUES(1)%si%imp                   = sources1%SI_IMP

    coresource(1)%VALUES(1)%qi%exp                   = sources1%QI_EXP * itm_ev    !! DPC 2009-07-02
    coresource(1)%VALUES(1)%qi%imp                   = sources1%QI_IMP

    coresource(1)%VALUES(1)%ui%exp                   = sources1%UI_EXP
    coresource(1)%VALUES(1)%ui%imp                   = sources1%UI_IMP



    RETURN


  END SUBROUTINE convert_local_to_cpo_types


! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
  SUBROUTINE analytics_full  (TIME, GEOMETRY1, PROFILES1, TRANSPORT1, SOURCES1)   


!-------------------------------------------------------!
!     This routine manufactures the solution for the    !
!     set of transport equations describing the main    !
!     plasma. It provides sources and transport         !
!     coefficients for all transport equations.         !
!-------------------------------------------------------!
!     Source:       ---                                 !
!     Developers:   R.Stankiewicz, D.Kalupin            !
!     Kontacts:     Roman.Stankiewich@gmail.com         !
!                   D.Kalupin@fz-juelich.de             !
!                                                       !
!     Comments:     ---                                 !
!                                                       !
!-------------------------------------------------------!


    USE itm_types
    USE ets_plasma
    USE analytical_functions
    use plasma_collisionality

    IMPLICIT NONE

    REAL (R8)           :: time                                  !Time


    INTEGER             :: nrho, irho                            !number of radial points (input) and current radial knot
    INTEGER             :: nion, iion                            !number of ion species   (input) and current ion type
    INTEGER             :: sigma_source
    INTEGER             :: ifail


! +++ Local derived types:
    TYPE (magnetic_geometry),      INTENT(INOUT)    :: geometry1   !contains all geometry quantities
    TYPE (plasma_profiles),        INTENT(INOUT)    :: profiles1   !contains profiles of plasma parameters
    TYPE (transport_coefficients), INTENT(INOUT)    :: transport1  !contains profiles of trasport coefficients
    TYPE (sources_and_sinks),      INTENT(INOUT)    :: sources1    !contains profiles of sources
    TYPE (collisionality)                           :: collisions1 !contains all terms determined by plasma collisions


! +++ Local variables:
!     Geometry:
    REAL (R8)           :: rho(profiles1%nrho),       hrho(profiles1%nrho)
    REAL (R8)           :: fdia(profiles1%nrho),      vpr(profiles1%nrho)
    REAL (R8)           :: g1(profiles1%nrho),        g2(profiles1%nrho)
    REAL (R8)           :: g3(profiles1%nrho)

!     Plasma composition:
    REAL (R8)           :: mion(profiles1%nion),      zion(profiles1%nion)

!     Plasma profiles:
    REAL (R8)           :: psi(profiles1%nrho) 
    REAL (R8)           :: ne(profiles1%nrho)
    REAL (R8)           :: ni(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: te(profiles1%nrho) 
    REAL (R8)           :: ti(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: vtor(profiles1%nrho,profiles1%nion)

!     Transport coefficients:
    REAL (R8)           :: sigma(profiles1%nrho)
    REAL (R8)           :: diff_ni(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: vconv_ni(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: diff_ti(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: vconv_ti(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: diff_te(profiles1%nrho)
    REAL (R8)           :: vconv_te(profiles1%nrho)
    REAL (R8)           :: diff_vtor(profiles1%nrho,profiles1%nion)
    REAL (R8)           :: vconv_vtor(profiles1%nrho,profiles1%nion)

!     Sources:
    REAL (R8)           :: vie(profiles1%nrho),       vei(profiles1%nrho)
    REAL (R8)           :: qie(profiles1%nrho),       qei(profiles1%nrho)
    REAL (R8)           :: vzi(profiles1%nrho),       uzi(profiles1%nrho)
    REAL (R8)           :: qzi(profiles1%nrho),       wzi(profiles1%nrho)

!     Other quantities:
    REAL (R8)           :: dne(profiles1%nrho), dtne(profiles1%nrho)
    REAL (R8)           :: c1
    REAL (R8)           :: bt, r                                  


! +++ Functions:
    REAL (R8)           :: x, t
    REAL (R8)           :: aval1, aval2, aval3, aval4, aval5, aval6, aval7
    REAL (R8)           :: gamma(profiles1%nion,profiles1%nrho)
    REAL (R8)           :: dgamma(profiles1%nion,profiles1%nrho)
    REAL (R8)           :: gammae(profiles1%nrho)
    REAL (R8)           :: dgammae(profiles1%nrho)


! +++ Set up dimensions:
    nrho         = SIZE(profiles1%NI)
    nion         = SIZE(profiles1%MION)
    c1           = transport1%C1(1)
    sigma_source = 0


    CALL allocate_collisionality(nrho, nion, collisions1, ifail)




!=================================================================================      
! +++ MAGNETIC GEOMETRY:       

!--------------------------------------------------------------------------------!
!     This part sets up the geometry used by analytical solution.                !
!--------------------------------------------------------------------------------!

    t                       = time 
    r                       = 1.0e0_r8
    bt                      = abt(t)


    DO irho=1,nrho
       x                     = arho(irho,nrho)
       rho(irho)             = x
       vpr(irho)             = avpr(x,t)
       g1(irho)              = ag1(x,t)
       g2(irho)              = ag2(x,t)
       g3(irho)              = ag3(x,t)
       fdia(irho)            = afdia(x,t)

       geometry1%RHO(irho)   = x
       geometry1%VPR(irho)   = vpr(irho)  
       geometry1%G1(irho)    = g1(irho)  
       geometry1%G2(irho)    = g2(irho)  
       geometry1%G3(irho)    = g3(irho)  
       geometry1%FDIA(irho)  = fdia(irho)
    ENDDO

    geometry1%RGEO           = r
    geometry1%BGEO           = bt






!=================================================================================      
! +++ PROFILES OF PLASMA PARAMETERS:

!--------------------------------------------------------------------------------!
!     This part describes analytical profiles for all plasma parameters.         !
!--------------------------------------------------------------------------------!

    t                             = time

    DO iion=1,nion
       zion(iion)                  = profiles1%ZION(iion)
       mion(iion)                  = profiles1%MION(iion)!*MP
    ENDDO


    DO irho=1,nrho
       x                           = rho(irho)
       te(irho)                    = ate(x,t)
       psi(irho)                   = apsi(x,t)

       profiles1%TE(irho)          = te(irho)
       profiles1%PSI(irho)         = psi(irho)
       profiles1%NE(irho)          = 0.e0_r8

       DO iion=1,nion
          ti(irho,iion)             = ati(iion,x,t)
          ni(irho,iion)             = ani(iion,x,t)
          vtor(irho,iion)           = avtor(iion,x,t)

          profiles1%NE(irho)        = profiles1%NE(irho)+profiles1%ZION(iion)*ni(irho,iion)
          profiles1%NI(irho,iion)   = ni(irho,iion)
          profiles1%TI(irho,iion)   = ti(irho,iion)
          profiles1%VTOR(irho,iion) = vtor(irho,iion)
       ENDDO

    ENDDO






!=================================================================================       
! +++ CALCULATING COLLISIONAL TERMS:
    CALL plasma_collisions(geometry1,profiles1,collisions1,ifail) 






!=================================================================================       
! +++ TRANSPORT_COEFFICIENT

!--------------------------------------------------------------------------------!
!     This part describes analytical profiles of transport coefficients.         !
!--------------------------------------------------------------------------------!

    DO irho=1,nrho
       x                                   = rho(irho)


       IF(sigma_source.EQ.1) THEN
          sigma(irho)                       = asigma(x,t)

          transport1%SIGMA(irho)            = sigma(irho)
       ENDIF


       IF(sigma_source.EQ.0) THEN
          sigma(irho)                       = collisions1%SIGMA(irho)
          transport1%SIGMA(irho)            = sigma(irho)
       ENDIF

    ENDDO




    DO irho=1,nrho
       x                                   = rho(irho)
       diff_te(irho)                       = akate(x,t)
       vconv_te(irho)                      = avte(x,t)

       transport1%DIFF_TE(irho)            = diff_te(irho)
       transport1%VCONV_TE(irho)           = vconv_te(irho)


       DO iion=1,nion
          diff_ni(irho,iion)                = ad(iion,x,t)
          vconv_ni(irho,iion)               = av(iion,x,t)
          diff_ti(irho,iion)                = akappae(iion,x,t)
          vconv_ti(irho,iion)               = avti(iion,x,t)
          diff_vtor(irho,iion)              = advtor(iion,x,t)
          vconv_vtor(irho,iion)             = aconvtor(iion,x,t)


          IF(c1.EQ.0.0e0_r8) THEN
             transport1%DIFF_NI(irho,iion,1) = diff_ni(irho,iion)
             transport1%VCONV_NI(irho,iion,1)= vconv_ni(irho,iion)
          ELSE IF(c1.EQ.1.5e0_r8) THEN
             transport1%DIFF_NI(irho,iion,2) = diff_ni(irho,iion)
             transport1%VCONV_NI(irho,iion,2)= vconv_ni(irho,iion)
          ELSE IF(c1.EQ.2.5e0_r8) THEN
             transport1%DIFF_NI(irho,iion,3) = diff_ni(irho,iion)
             transport1%VCONV_NI(irho,iion,3)= vconv_ni(irho,iion)
          END IF
          transport1%DIFF_TI(irho,iion)     = diff_ti(irho,iion)
          transport1%VCONV_TI(irho,iion)    = vconv_ti(irho,iion)
          transport1%DIFF_VTOR(irho,iion)   = diff_vtor(irho,iion)        
          transport1%VCONV_VTOR(irho,iion)  = vconv_vtor(irho,iion)
       ENDDO

    ENDDO



!=================================================================================        
! +++ SOURCES_AND_SINKS 
!--------------------------------------------------------------------------------!
!     CALCULATION OF EXPLICIT SOURCES FROM TRANSPORT EQUATION USING THE          !
!     ANALATICAL FORMULAE FOR SOLUTIONS AND TRANSPORT COEFFICIENTS               !
!--------------------------------------------------------------------------------!





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ CURRENT TRANSPORT:
    t                        = time

    DO irho=2,nrho
       x                      = rho(irho)

       aval1                  = sigma(irho)*dtapsi(x,t)

       aval2                  = -x*dtabt(t)/2.0e0_r8/abt(t)*dapsi(x,t)*sigma(irho)
       aval3                  = fdia(irho)**2/itm_mu0/abt(t)/x

       aval4                  = vpr(irho)/(4.0e0_r8*itm_pi**2)*g3(irho)/fdia(irho)
       aval5                  = davpr(x,t)/(4.0e0_r8*itm_pi**2)*g3(irho)/fdia(irho)      
       aval5                  = aval5+avpr(x,t)/(4.0e0_r8*itm_pi**2)*dag3(x,t)/fdia(irho)
       aval5                  = aval5-avpr(x,t)/(4.0e0_r8*itm_pi**2)*g3(irho)/fdia(irho)**2*dafdia(x,t)
       aval1                  = aval1+aval2-aval3*aval4*ddapsi(x,t)-aval3*aval5*dapsi(x,t)


       aval1                  = aval1*(2.e0_r8*itm_pi*x)/vpr(irho)+acurr_imp(x,t)*psi(irho)


       sources1%CURR_EXP(irho)= -aval1
       sources1%CURR_IMP(irho)= acurr_imp(x,t)

    ENDDO



    sources1%CURR_EXP(1)     = sources1%CURR_EXP(2)
    sources1%CURR_IMP(1)     = acurr_imp(0.e0_r8,t)



! +++ BOUNDARY CONDITION
    IF(profiles1%PSI_BND_TYPE.EQ.0) THEN
       DO irho=2,nrho
          profiles1%QSF(irho)  = 2.e0_r8*itm_pi*rho(irho)*abt(t)/dapsi(rho(irho),t)
       ENDDO

       profiles1%QSF(1)       = 2.e0_r8*itm_pi*abt(t)/ddapsi(rho(1),t)                
    ENDIF


    IF(profiles1%PSI_BND_TYPE.EQ.1) THEN
       profiles1%PSI_BND(1)    = psi(nrho)
    ENDIF

    IF(profiles1%PSI_BND_TYPE.EQ.2) THEN
       profiles1%PSI_BND(1)    = vpr(nrho)*g3(nrho)*dapsi(rho(nrho),t)/4.e0_r8/itm_pi**2/itm_mu0
    ENDIF

    IF(profiles1%PSI_BND_TYPE.EQ.3) THEN
       profiles1%PSI_BND(1)    = dtapsi(rho(nrho),t)
    ENDIF





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ PARTICLE TRANSPORT

    DO iion=1,nion
       t                             = time
       gamma(iion,1)                 = 0.0e0_r8               !DPC 2009-01-19
       dgamma(iion,1)                = 0.0e0_r8               !DPC 2009-01-19
       DO irho=2,nrho
          x                           = rho(irho)
          aval1                       = dtavpr(x,t)*ani(iion,x,t)+avpr(x,t)*dtani(iion,x,t)
          aval2                       = -dtabt(t)/2.e0_r8/abt(t)*(avpr(x,t)* ani(iion,x,t)+ &
               x *(davpr(x,t)* ani(iion,x,t)+avpr(x,t)* dani(iion,x,t)))
          gamma(iion,irho)            = avpr(x,t)* ag1(x,t)*(-ad(iion,x,t)* dani(iion,x,t)+ani(iion,x,t)*av(iion,x,t))
          dgamma(iion,irho)           = avpr(x,t)* ag1(x,t)*(-ad(iion,x,t)* ddani(iion,x,t)+dani(iion,x,t)*av(iion,x,t))  

          dgamma(iion,irho)           = dgamma(iion,irho)+ &
               davpr(x,t)* ag1(x,t)*(-ad(iion,x,t)* dani(iion,x,t)+av(iion,x,t)*ani(iion,x,t))       
          dgamma(iion,irho)           = dgamma(iion,irho)+ &
               avpr(x,t)* dag1(x,t)*(-ad(iion,x,t)* dani(iion,x,t)+av(iion,x,t)*ani(iion,x,t))         
          dgamma(iion,irho)           = dgamma(iion,irho)+ &
               avpr(x,t)* ag1(x,t)*(-dad(iion,x,t)* dani(iion,x,t)+dav(iion,x,t)*ani(iion,x,t))           

          aval1                       = aval1+aval2+dgamma(iion,irho)+avpr(x,t)*ani(iion,x,t)*asi_imp(iion,x,t)
          sources1%SI_EXP(irho,iion)  = aval1/avpr(x,t)
          sources1%SI_IMP(irho,iion)  = asi_imp(iion,x,t)
       ENDDO


       sources1%SI_EXP(1,iion)       = sources1%SI_EXP(2,iion)
       sources1%SI_IMP(1,iion)       = asi_imp(iion,0.e0_r8,t)


! +++ BOUNDARY CONDITION
       IF (profiles1%NI_BND_TYPE(iion).EQ.1) THEN
          profiles1%NI_BND(1,iion)   = ani(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%NI_BND_TYPE(iion).EQ.2) THEN
          profiles1%NI_BND(1,iion)   = -dani(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%NI_BND_TYPE(iion).EQ.3) THEN
          profiles1%NI_BND(1,iion)   = -ani(iion,rho(nrho),t)/dani(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%NI_BND_TYPE(iion).EQ.4) THEN
          profiles1%NI_BND(1,iion)   = gamma(iion,nrho)   
       ENDIF

    ENDDO





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ ELECTRON DENSISTY AND FLUXES
    DO irho =1,nrho
       x                           = rho(irho)
       ne(irho)                    = 0.e0_r8
       dne(irho)                   = 0.e0_r8
       dtne(irho)                  = 0.e0_r8
       gammae(irho)                = 0.e0_r8
       dgammae(irho)               = 0.e0_r8

       DO iion=1,nion 
          ne(irho)                  = ne(irho) + zion(iion)*ni(irho,iion)
          dtne(irho)                = dtne(irho) + zion(iion)*dtani(iion,x,t)
          dne(irho)                 = dne(irho) + zion(iion)*dani(iion,x,t)
          gammae(irho)              = gammae(irho) + gamma(iion,irho)*zion(iion)
          dgammae(irho)             = dgammae(irho) + dgamma(iion,irho)*zion(iion)
       ENDDO

    ENDDO





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ ION ENERGY TRANSPORT
    DO iion=1,nion
       t                            = time

       DO irho=2,nrho
          x                          = rho(irho)
          aval1                      = dtavpr(x,t)*avpr(x,t)**(2.e0_r8/3.e0_r8)*ani(iion,x,t)*ati(iion,x,t)*5.e0_r8/3.e0_r8
          aval1                      = aval1+avpr(x,t)**(5.e0_r8/3.e0_r8)*dtani(iion,x,t)*ati(iion,x,t)
          aval1                      = aval1+avpr(x,t)**(5.e0_r8/3.e0_r8)*ani(iion,x,t)*dtati(iion,x,t)

          aval2                      = -dtabt(t)/2.e0_r8/abt(t)

          aval3                      = ani(iion,x,t)*ati(iion,x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
          aval3                      = aval3+ x*dani(iion,x,t)*ati(iion,x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
          aval3                      = aval3+ x*ani(iion,x,t)*dati(iion,x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
          aval3                      = aval3+ x*ani(iion,x,t)*ati(iion,x,t)*avpr(x,t)**(2.e0_r8/3.e0_r8)*davpr(x,t)*5.e0_r8/3.e0_r8

          aval1                      = aval1+aval2*aval3
          aval1                      = aval1*3.e0_r8/2.e0_r8

          aval2                      = avpr(x,t)*ag1(x,t)*ani(iion,x,t)                                    

          aval3                      = -akappae(iion,x,t)*dati(iion,x,t)+ati(iion,x,t)*avti(iion,x,t)

          aval4                      = davpr(x,t)*ag1(x,t)*ani(iion,x,t)+ &
               avpr(x,t)*dag1(x,t)*ani(iion,x,t)+avpr(x,t)*ag1(x,t)*dani(iion,x,t)

          aval5                      = -dakappae(iion,x,t)*dati(iion,x,t)+ati(iion,x,t)*davti(iion,x,t)
          aval5                      = aval5-akappae(iion,x,t)*ddati(iion,x,t)+dati(iion,x,t)*avti(iion,x,t)

          aval6                      = aval2*aval5+aval4*aval3

          aval1                      = aval1+aval6*avpr(x,t)**(2.e0_r8/3.e0_r8)


          aval2                      = dati(iion,x,t)*gamma(iion,irho)+ati(iion,x,t)*dgamma(iion,irho)

          aval1                      = aval1+c1*aval2*avpr(x,t)**(2.e0_r8/3.e0_r8)
          aval1                      = aval1/avpr(x,t)**(5.e0_r8/3.e0_r8)

          vei(irho)                  = collisions1%VEI(irho,iion)
          qei(irho)                  = collisions1%QEI(irho,iion)
          vzi(irho)                  = collisions1%VZI(irho,iion)
          qzi(irho)                  = collisions1%QZI(irho,iion)

          aval2                      = (aqi_imp(iion,x,t) + vei(irho) + vzi(irho) )                                         

          sources1%QI_EXP(irho,iion) = aval1- qei(irho) - qzi(irho) +aval2*ati(iion,x,t)
          sources1%QI_IMP(irho,iion) = aqi_imp(iion,x,t)
       ENDDO

       sources1%QI_EXP(1,iion)      = sources1%QI_EXP(2,iion)
       sources1%QI_IMP(1,iion)      = aqi_imp(iion,0.e0_r8,t)
    ENDDO




! +++ BOUNDARY CONDITION
    DO iion=1,nion
       IF (profiles1%TI_BND_TYPE(iion).EQ.1)THEN
          profiles1%TI_BND(1,iion)  = ati(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%TI_BND_TYPE(iion).EQ.2)THEN
          profiles1%TI_BND(1,iion)  = -dati(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%TI_BND_TYPE(iion).EQ.3)THEN
          profiles1%TI_BND(1,iion)  = -ati(iion,rho(nrho),t)/dati(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%TI_BND_TYPE(iion).EQ.4)THEN
          x                         = rho(nrho)
          aval2                     = avpr(x,t)*ag1(x,t)*ani(iion,x,t) 
          aval3                     = (-akappae(iion,x,t)*dati(iion,x,t)+ati(iion,x,t)*avti(iion,x,t))* &
               aval2+c1*ati(iion,x,t)*gamma(iion,nrho)
          profiles1%TI_BND(1,iion)  = aval3   
       ENDIF

    ENDDO






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ ELECTRON ENERGY

    t                              = time

    DO irho=2,nrho
       x                             = rho(irho)
       aval1                         = dtavpr(x,t)*avpr(x,t)**(2.e0_r8/3.e0_r8)*ne(irho)*ate(x,t)*5.e0_r8/3.e0_r8
       aval1                         = aval1+avpr(x,t)**(5.e0_r8/3.e0_r8)*dtne(irho)*ate(x,t)
       aval1                         = aval1+avpr(x,t)**(5.e0_r8/3.e0_r8)*ne(irho)*dtate(x,t)
       aval2                         = -dtabt(t)/2.e0_r8/abt(t)
       aval3                         = ne(irho)*ate(x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
       aval3                         = aval3+ x*dne(irho)*ate(x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
       aval3                         = aval3+ x*ne(irho)*date(x,t)*avpr(x,t)**(5.e0_r8/3.e0_r8)
       aval3                         = aval3+ x*ne(irho)*ate(x,t)*avpr(x,t)**(2.e0_r8/3.e0_r8)*davpr(x,t)*5.e0_r8/3.e0_r8
       aval1                         = aval1+aval2*aval3
       aval1                         = aval1*3.e0_r8/2.e0_r8
       aval2                         = avpr(x,t)*ag1(x,t)*ne(irho)                        
       aval3                         = -akate(x,t)*date(x,t)+ate(x,t)*avte(x,t)
       aval4                         = davpr(x,t)*ag1(x,t)*ne(irho)+avpr(x,t)*dag1(x,t)*ne(irho)+avpr(x,t)*ag1(x,t)*dne(irho)
       aval5                         = -dakate(x,t)*date(x,t)+ate(x,t)*davte(x,t)
       aval5                         = aval5-akate(x,t)*ddate(x,t)+date(x,t)*avte(x,t)
       aval6                         = aval2*aval5+aval4*aval3
       aval1                         = aval1+aval6*avpr(x,t)**(2.e0_r8/3.e0_r8)
       c1                            = 0.e0_r8
       aval2                         = c1*(avpr(x,t)**(2.e0_r8/3.e0_r8)*(date(x,t)*gammae(irho)+ ate(x,t)*dgammae(irho)))
       aval1                         = aval1+aval2
       aval1                         = aval1/avpr(x,t)**(5.e0_r8/3.e0_r8)
       vie(irho)                     = collisions1%VIE(irho)
       qie(irho)                     = collisions1%QIE(irho)
       aval2                         = (aqe_imp(x,t) + vie(irho)  )                                         
       sources1%QE_EXP(irho)         = aval1 +aval2*ate(x,t)-qie(irho)
       sources1%QE_IMP(irho)         = aqe_imp(x,t)
    ENDDO

    sources1%QE_EXP(1)              = sources1%QE_EXP(2)
    sources1%QE_IMP(1)              = sources1%QE_IMP(2)



! +++ BOUNDARY CONDITION
    IF (profiles1%TE_BND_TYPE.EQ.1)THEN
       profiles1%TE_BND(1)          = te(nrho)    
    ENDIF

    IF (profiles1%TE_BND_TYPE.EQ.2)THEN
       profiles1%TE_BND(1)          = -date(rho(nrho),t)   
    ENDIF

    IF (profiles1%TE_BND_TYPE.EQ.3)THEN
       profiles1%TE_BND(1)          = -te(nrho)/date(rho(nrho),t)    
    ENDIF

    IF (profiles1%TE_BND_TYPE.EQ.4)THEN
       x                            = rho(nrho)
       aval2                        = avpr(x,t)*ag1(x,t)*ne(nrho) 
       aval3                        = -akate(x,t)*date(x,t)+ate(rho(nrho),t)*avte(x,t)
       profiles1%TE_BND(1)          = aval2*aval3+c1*ate(rho(nrho),t)*gammae(nrho)
    ENDIF





! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! +++ ROTATION TRANSPORT
    DO iion=1,nion
       t                            = time

       DO irho=2,nrho
          x                          = rho(irho)
          aval1                      = dtavpr(x,t)*ani(iion,x,t)*avtor(iion,x,t)*ag2(x,t)
          aval1                      = aval1+avpr(x,t)*dtani(iion,x,t)*avtor(iion,x,t)*ag2(x,t)
          aval1                      = aval1+avpr(x,t)*ani(iion,x,t)*dtavtor(iion,x,t)*ag2(x,t)
          aval1                      = aval1+avpr(x,t)*ani(iion,x,t)*avtor(iion,x,t)*dtag2(x,t)

          aval2                      = -dtabt(t)/2.e0_r8/abt(t)

          aval3                      = ani(iion,x,t)*avtor(iion,x,t)*avpr(x,t)*ag2(x,t)
          aval3                      = aval3+ x*dani(iion,x,t)*avtor(iion,x,t)*avpr(x,t)*ag2(x,t)
          aval3                      = aval3+ x*ani(iion,x,t)*davtor(iion,x,t)*avpr(x,t)*ag2(x,t)
          aval3                      = aval3+ x*ani(iion,x,t)*avtor(iion,x,t)*davpr(x,t)*ag2(x,t)
          aval3                      = aval3+ x*ani(iion,x,t)*avtor(iion,x,t)*davpr(x,t)*dag2(x,t)

          aval1                      = aval1+aval2*aval3

          aval2                      = avpr(x,t)*ag2(x,t)*ani(iion,x,t)*ag1(x,t)                                    

          aval3                      = -advtor(iion,x,t)*davtor(iion,x,t)+avtor(iion,x,t)*aconvtor(iion,x,t)

          aval4                      = davpr(x,t)*ag2(x,t)*ani(iion,x,t)+ &
               avpr(x,t)*dag2(x,t)*ani(iion,x,t)+avpr(x,t)*ag2(x,t)*dani(iion,x,t)
          aval4                      = aval4*ag1(x,t)+avpr(x,t)*ag2(x,t)*ani(iion,x,t)*dag1(x,t)

          aval5                      = -dadvtor(iion,x,t)*davtor(iion,x,t)+avtor(iion,x,t)*daconvtor(iion,x,t)
          aval5                      = aval5-advtor(iion,x,t)*ddavtor(iion,x,t)+davtor(iion,x,t)*aconvtor(iion,x,t)

          aval6                      = aval2*aval5+aval4*aval3 +dag2(x,t)*gamma(iion,irho)*avtor(iion,x,t)!*AVAL4
          aval6                      = aval6+(ag2(x,t)*dgamma(iion,irho)*avtor(iion,x,t)+ &
               ag2(x,t)*gamma(iion,irho)*davtor(iion,x,t))!*AVAL4

          aval1                      = aval1+aval6
          aval1                      = aval1*profiles1%MION(iion)/avpr(x,t)



          uzi(irho)                  = collisions1%UZI(irho,iion)
          wzi(irho)                  = collisions1%WZI(irho,iion)
          aval2                      = aui_imp(iion,x,t)  + wzi(irho)                                        

          sources1%UI_EXP(irho,iion) = aval1 +aval2*avtor(iion,x,t) - uzi(irho)
          sources1%UI_IMP(irho,iion) = aui_imp(iion,x,t)
       ENDDO

       sources1%UI_EXP(1,iion)      = sources1%UI_EXP(2,iion)
       sources1%UI_EXP(1,iion)      = aui_imp(iion,0.e0_r8,t)
       sources1%UI_IMP(1,iion)      = aui_imp(iion,0.e0_r8,t)
    ENDDO




! +++ BOUNDARY CONDITION
    DO iion=1,nion
       IF (profiles1%VTOR_BND_TYPE(iion).EQ.1)THEN
          profiles1%VTOR_BND(1,iion) = avtor(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%VTOR_BND_TYPE(iion).EQ.2)THEN
          profiles1%VTOR_BND(1,iion) = -davtor(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%VTOR_BND_TYPE(iion).EQ.3)THEN
          profiles1%VTOR_BND(1,iion) = -avtor(iion,rho(nrho),t)/davtor(iion,rho(nrho),t)    
       ENDIF

       IF (profiles1%VTOR_BND_TYPE(iion).EQ.4)THEN
          x                          = rho(nrho)
          aval2                      = avpr(x,t)*ag2(x,t)*ani(iion,x,t)*ag1(x,t) 
          aval3                      = -advtor(iion,x,t)*davtor(iion,x,t)+avtor(iion,rho(nrho),t)*aconvtor(iion,x,t)
          profiles1%VTOR_BND(1,iion) = aval2*aval3 +avtor(iion,rho(nrho),t)*ag2(x,t)*gamma(iion,nrho)
          profiles1%VTOR_BND(1,iion) = profiles1%VTOR_BND(1,iion)*mion(iion)  
       ENDIF

    ENDDO



    CALL deallocate_collisionality(collisions1, ifail) 


    RETURN



  END SUBROUTINE analytics_full

! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  





END MODULE analytics1