coronal.f90

Go to the documentation of this file.

MODULE coronal

CONTAINS
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   
  SUBROUTINE  set_coronal(COREIMPUR_IN, COREPROF_IN, COREIMPUR_OUT, INTERPOL, ICORONAL) 


! +++ Declaration of variables: 
    USE euitm_schemas
    USE euitm_routines
    USE copy_structures
    USE allocate_deallocate
    USE itm_types


    IMPLICIT NONE

    INTEGER,              PARAMETER  :: nslice = 1             !number of CPO ocurancies in the work flow
    INTEGER                          :: nrho, nrho2, irho
    INTEGER                          :: nnucl,inucl
    INTEGER                          :: nion, nzimp1
    INTEGER                          :: nimp, iimp  
    INTEGER,             ALLOCATABLE :: nzimp(:)
    INTEGER                          ::           izimp
    INTEGER                          :: nneut
    INTEGER,             ALLOCATABLE :: ncomp(:)
    INTEGER,             ALLOCATABLE :: ntype(:)


    INTEGER                          :: interpol           !interpolation index "0"-based on RHO_TOR; all other - based on RHO_TOR_NORM 
    INTEGER                          :: icoronal           !interpolation index: "0"-OFF; "1" - replace boundary conditions by coronal;   "2" - replace boundary conditions and profiles by coronal 

! +++ CPO derived types:
    TYPE (type_coreimpur), POINTER   :: coreimpur_in(:)    !input CPO with impurity 
    TYPE (type_coreprof),  POINTER   :: coreprof_in(:)     !input CPO with plasma profiles
    TYPE (type_coreimpur), POINTER   :: coreimpur_out(:)   !output CPO with sources uploaded from the data base 


! +++ Internal derived types:
    REAL (R8),          ALLOCATABLE  :: rho(:),  rho2(:), ne(:), te(:), n_impurity(:,:), nzt_impurity(:)   
    REAL (R8)                        :: nimp_tot, amn, zn   



! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   
! +++ allocate and define grid of output CPO:
    nrho               = SIZE (coreimpur_in(1)%rho_tor, dim=1)
    nrho2              = SIZE (coreprof_in(1)%rho_tor,  dim=1)

    CALL get_comp_dimensions(coreimpur_in(1)%COMPOSITIONS, nnucl, nion,  nimp,  nzimp, nneut, ntype, ncomp)
    IF(.NOT.(ASSOCIATED(coreimpur_out)))&
!!    CALL ALLOCATE_COREIMPUR_CPO   (NSLICE,      NRHO, NNUCL, NION,  NIMP,  NZIMP, NNEUT, NTYPE, NCOMP, COREIMPUR_OUT)     !! overwritten by next line
    CALL copy_cpo(coreimpur_in, coreimpur_out)


    IF(icoronal .EQ. 0) goto 10


    ALLOCATE ( rho(nrho))
    ALLOCATE ( rho2(nrho2))

    IF (interpol.EQ.0) THEN
       rho             = coreimpur_in(1)%rho_tor
       rho2            = coreprof_in(1)%rho_tor
    ELSE
       rho             = coreimpur_in(1)%rho_tor / coreimpur_in(1)%rho_tor(nrho)
       rho2            = coreprof_in(1)%rho_tor  / coreprof_in(1)%rho_tor(nrho2)
    END IF




! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   
    ALLOCATE ( te(nrho))
    ALLOCATE ( ne(nrho))
    
    IF (ASSOCIATED(coreprof_in(1)%ne%value))                                                      &
    CALL l3interp(coreprof_in(1)%ne%value,  rho2, nrho2,       ne,  rho, nrho)
    IF (ASSOCIATED(coreprof_in(1)%te%value))                                                      &
    CALL l3interp(coreprof_in(1)%te%value,  rho2, nrho2,       te,  rho, nrho)
 
    IF (coreprof_in(1)%te%boundary%type .EQ. 1)  te(nrho)    =  coreprof_in(1)%te%boundary%value(1)



! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   
    DO iimp = 1, nimp
       inucl          =      coreimpur_in(1)%COMPOSITIONS%IMPURITIES(iimp)%nucindex

       zn             =      coreimpur_in(1)%COMPOSITIONS%NUCLEI(inucl)%zn
       amn            =      coreimpur_in(1)%COMPOSITIONS%NUCLEI(inucl)%amn
       nimp_tot       =      sum(coreimpur_in(1)%IMPURITY(iimp)%boundary%value(1,:))
       write(*,*) 'set_coronal: iimp, NIMP_TOT = ', iimp, nimp_tot

       IF (ALLOCATED(n_impurity)) DEALLOCATE(n_impurity)
       ALLOCATE     (n_impurity(nrho,nzimp(iimp)))

       IF (ALLOCATED(nzt_impurity)) DEALLOCATE(nzt_impurity)
       ALLOCATE     (nzt_impurity(nrho))

                                                                

       nzimp1        =       nzimp(iimp)
       CALL coronal_distribution(nrho, ne, te, amn, zn, nzimp1, n_impurity)

       DO irho = 1,nrho
          DO izimp = 1,nzimp(iimp)
             IF (n_impurity(nrho,izimp).LE.1.0d-200) n_impurity(nrho,izimp) = 0._r8 
          ENDDO
       ENDDO

!      Replace boundary conditions by coronal
       IF(icoronal .EQ. 1) THEN
          DO izimp = 1,nzimp(iimp)
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(1,izimp)    = n_impurity(nrho,izimp)*nimp_tot 
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(2,izimp)    = 0.0_r8
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(3,izimp)    = 0.0_r8
          END DO
          
          
!      Replace boundary conditions and concentration profiles by coronal
       ELSE IF(icoronal .EQ. 2) THEN
          DO izimp = 1,nzimp(iimp)
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(1,izimp)    = n_impurity(nrho,izimp)*nimp_tot 
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(2,izimp)    = 0.0_r8
             coreimpur_out(1)%IMPURITY(iimp)%boundary%value(3,izimp)    = 0.0_r8

             coreimpur_out(1)%IMPURITY(iimp)%nz(:,izimp)                = n_impurity(:,izimp)*nimp_tot 

             coreimpur_out(1)%IMPURITY(iimp)%z(:,izimp)                 = izimp 
             coreimpur_out(1)%IMPURITY(iimp)%zsq(:,izimp)               = izimp**2

          END DO

       ! experimental
       ELSE IF(icoronal .EQ. 3) THEN
          IF (ASSOCIATED(coreimpur_in(1)%impurity(iimp)%nz) .and.   &
              ASSOCIATED(coreimpur_in(1)%impurity(iimp)%z))  THEN 
                !NZT_IMPURITY(:) = SUM(COREIMPUR_IN(1)%impurity(IIMP)%z(:,:) * COREIMPUR_IN(1)%impurity(IIMP)%nz(:,:),   &
                !                      dim=2)
                nzt_impurity(:) = sum(coreimpur_in(1)%impurity(iimp)%nz(:,:),   &
                                      dim=2)
                
               ! N_IMPURITY(:,IZIMP)*NZT_IMPURITY(:)
               ! NZ_IMPURITY(:) = SUM(COREIMPUR_IN(1)%impurity(IIMP)%z(:,:) * COREIMPUR_IN(1)%impurity(IIMP)%nz(:,:),   &

                DO izimp = 1,nzimp(iimp)
                   coreimpur_out(1)%IMPURITY(iimp)%nz(:,izimp)            = n_impurity(:,izimp)*nzt_impurity(:)

                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(1,izimp)    = coreimpur_out(1)%IMPURITY(iimp)%nz(nrho,izimp) 
                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(2,izimp)    = 0.0_r8
                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(3,izimp)    = 0.0_r8

                   coreimpur_out(1)%IMPURITY(iimp)%z(:,izimp)                 = izimp 
                   coreimpur_out(1)%IMPURITY(iimp)%zsq(:,izimp)               = izimp**2

                END DO
           ELSE 
                DO izimp = 1,nzimp(iimp)
                   coreimpur_out(1)%IMPURITY(iimp)%nz(:,izimp)                = 0.0_r8

                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(1,izimp)    = 0.0_r8 
                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(2,izimp)    = 0.0_r8
                   coreimpur_out(1)%IMPURITY(iimp)%boundary%value(3,izimp)    = 0.0_r8

                   coreimpur_out(1)%IMPURITY(iimp)%z(:,izimp)                 = izimp 
                   coreimpur_out(1)%IMPURITY(iimp)%zsq(:,izimp)               = izimp**2

                END DO
           ENDIF
          
       END IF

    END DO




    IF (ALLOCATED(n_impurity)) DEALLOCATE (n_impurity)
    IF (ALLOCATED(nzt_impurity)) DEALLOCATE (nzt_impurity)
    DEALLOCATE (te)
    DEALLOCATE (ne)
    DEALLOCATE (rho)
    DEALLOCATE (rho2)


 10 RETURN


    END SUBROUTINE  set_coronal 
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   






! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   
SUBROUTINE coronal_distribution (NRHO, NE, TE, AMN,  ZN, NZIMP1, N_IMPURITY)
  USE itm_types
  USE amns_types
  USE amns_module


  IMPLICIT NONE

  INTEGER                                 :: nzimp1
  INTEGER                                 :: nrho, irho
  TYPE (amns_handle_type)                 :: amns
  type (amns_handle_rx_type), ALLOCATABLE :: amns_ei(:), amns_rc(:), amns_lr(:), amns_br(:)
  type (amns_reaction_type)               :: xx_rx
  TYPE (amns_reactants_type)              :: species
  TYPE (amns_query_type)                  :: query
  TYPE (amns_answer_type)                 :: answer
  TYPE (amns_set_type)                    :: set
  REAL (kind=R8)                          :: te(nrho), ne(nrho), n_impurity(nrho,nzimp1)
  REAL (kind=R8),             ALLOCATABLE :: na(:), rhs(:)
  REAL (kind=R8),             ALLOCATABLE :: rate_ei(:,:), rate_rc(:,:), rate_lr(:,:), rate_br(:,:)
  REAL (kind=R8),             ALLOCATABLE :: l(:), d(:), u(:)
  REAL (kind=R8)                          :: zn, mi, amn
  REAL (kind=R8)                          :: test, norm, line_radiation, recombination_radiation
  CHARACTER (len=12),         ALLOCATABLE :: state_labels(:)
  INTEGER                                 :: ns, is, isref
!================================================================================



  mi = amn
  mi = 0
  ns = nzimp1

! some initializations
  ALLOCATE(amns_ei(0:ns), amns_rc(0:ns), amns_lr(0:ns), amns_br(0:ns))
  ALLOCATE(state_labels(0:ns))
  ALLOCATE(l(0:ns), d(0:ns), u(0:ns), na(0:ns), rhs(0:ns))
  ALLOCATE(rate_ei(1:nrho,0:ns), rate_rc(1:nrho,0:ns), rate_lr(1:nrho,0:ns), rate_br(1:nrho,0:ns))

! set up the AMNS system

  CALL itm_amns_setup(amns)
  WRITE(*,*) 'Done ITM_AMNS_SETUP'

! set up tables for ionization
  xx_rx%string='EI'
  DO is=0,ns
     IF(is .NE. zn) THEN
        allocate(species%components(4))
        species%components =  &
             (/ amns_reactant_type(zn, is, mi, 0), &
                amns_reactant_type(0, -1, 0, 0), &
                amns_reactant_type(zn, is+1, mi, 1), &
                amns_reactant_type(0, -1, 0, 1) &
              /)
        CALL itm_amns_setup_table(amns, xx_rx, species, amns_ei(is))
        deallocate(species%components)
        query%string='state_label'
        CALL itm_amns_query_table(amns_ei(is),query,answer)
        state_labels(is)=answer%string
     ENDIF
  ENDDO

! set up tables for line radiation
  xx_rx%string='LR'
  DO is=0,ns
     IF(is .NE. zn) THEN
        allocate(species%components(2))
        species%components = &
             (/ amns_reactant_type(zn, is, mi, 0), &
                amns_reactant_type(zn, is, mi, 1) &
             /)
        CALL itm_amns_setup_table(amns, xx_rx, species, amns_lr(is))
        deallocate(species%components)
        query%string='state_label'
        CALL itm_amns_query_table(amns_lr(is),query,answer)
        state_labels(is)=answer%string
     ENDIF
  ENDDO

! set up tables for recombination
  xx_rx%string='RC'
  DO is=0,ns
     IF(is .NE. 0) THEN
        allocate(species%components(4))
        species%components =  &
             (/ amns_reactant_type(zn, is, mi, 0), &
                amns_reactant_type(0, -1, 0, 0), &
                amns_reactant_type(zn, is-1, mi, 1), &
                amns_reactant_type(0, -1, 0, 1) &
              /)
        CALL itm_amns_setup_table(amns, xx_rx, species, amns_rc(is))
        deallocate(species%components)
        query%string='state_label'
        CALL itm_amns_query_table(amns_rc(is),query,answer)
        state_labels(is)=answer%string
     ENDIF
  ENDDO

! set up tables for recombination radiation
  xx_rx%string='BR'
  DO is=0,ns
     IF(is .NE. 0) THEN
        allocate(species%components(2))
        species%components = &
             (/ amns_reactant_type(zn, is, mi, 0), &
                amns_reactant_type(zn, is, mi, 1) &
             /)
        CALL itm_amns_setup_table(amns, xx_rx, species, amns_br(is))
        deallocate(species%components)
        query%string='state_label'
        CALL itm_amns_query_table(amns_br(is),query,answer)
        state_labels(is)=answer%string
     ENDIF
  ENDDO

  WRITE(*,*) 'Done ITM_AMNS_SETUP_TABLE'

! interpolate ionization and line radiation data
  set%string='nowarn'
  DO is=0,ns
     IF(is .NE. zn) THEN
         CALL itm_amns_set_table(amns_ei(is),set)
         CALL itm_amns_rx(amns_ei(is),rate_ei(:,is),te,ne)
         CALL itm_amns_set_table(amns_lr(is),set)
         CALL itm_amns_rx(amns_lr(is),rate_lr(:,is),te,ne)
     ENDIF
  ENDDO

! interpolate recombination and recombination radiation data
  set%string='nowarn'
  DO is=0,ns
     IF(is .NE. 0) THEN
         CALL itm_amns_set_table(amns_rc(is),set)
         CALL itm_amns_rx(amns_rc(is),rate_rc(:,is),te,ne)
         CALL itm_amns_set_table(amns_br(is),set)
         CALL itm_amns_rx(amns_br(is),rate_br(:,is),te,ne)
     ENDIF
  ENDDO

  WRITE(*,*) 'Done ITM_AMNS_RX'

! finalize the tables
  DO is=0,ns
     IF(is .NE. zn) THEN
        CALL itm_amns_finish_table(amns_ei(is))
        CALL itm_amns_finish_table(amns_lr(is))
     ENDIF
  ENDDO
  DO is=0,ns
     IF(is .NE. 0) THEN
        CALL itm_amns_finish_table(amns_rc(is))
        CALL itm_amns_finish_table(amns_br(is))
     ENDIF
  ENDDO

  WRITE(*,*) 'Done ITM_AMNS_FINISH_TABLE'

! finalize the system
  CALL itm_amns_finish(amns)

  WRITE(*,*) 'Done ITM_AMNS_FINISH'

  OPEN(10,file='coronal.out')
  WRITE(10,'(100a15)') '#te            ',state_labels, &
                       '   LR          ',  &
                       '   BR          '

  DO irho=1, nrho

! set up the matrix
     l(0)=0.0_r8
     d(0)=-rate_ei(irho,0)
     u(0)=rate_rc(irho,1)
     rhs(0)=0.0_r8
     DO is=1,ns-1
        l(is)=rate_ei(irho,is-1)
        d(is)=-rate_ei(irho,is)-rate_rc(irho,is)
        u(is)=rate_rc(irho,is+1)
        rhs(is)=0.0_r8
     ENDDO
     l(ns)=rate_ei(irho,ns-1)
     d(ns)=-rate_rc(irho,ns)
     u(ns)=0.0_r8
     rhs(ns)=0.0_r8

! we need to set the value of 1 charge state
     isref=0
     DO is=1, ns-1
        IF(rate_ei(irho,is) .LT. rate_rc(irho,is)) THEN
           EXIT
        ELSE
           isref=isref+1
        ENDIF
     ENDDO
!     write(*,*) te(IRHO), isref
     u(isref)=0.0_r8
     d(isref)=1.0_r8
     l(isref)=0.0_r8
     rhs(isref)=1.0_r8

! solve (l,d,u) na = rhs
     CALL solve_tridiag(l,d,u,rhs,na,ns+1)
     na(:) = na(:) / sum(na)

     line_radiation = sum(na(0:ns-1)*rate_lr(irho,0:ns-1))
     recombination_radiation = sum(na(1:ns)*rate_br(irho,1:ns))

! OUTPUT DENSITIES:
     DO is=1,ns
        n_impurity(irho,is) = na(is)
     ENDDO 


! check the answer

     norm=max(maxval(rate_ei(irho,0:ns-1) * na(0:ns-1)), maxval(rate_rc(irho,1:ns) * na(1:ns)))
     is=0
     test = (-rate_ei(irho,is) * na(is) + rate_rc(irho,is+1) * na(is+1))/norm
     IF(abs(test).GT.1e-15_r8) WRITE(*,*) 'LARGE ERROR: ',irho, is, test
     DO is=1, ns-1
        test = (rate_ei(irho,is-1) * na(is-1) - (rate_rc(irho,is) + rate_ei(irho,is)) * na(is) + rate_rc(irho,is+1) * na(is+1))/norm
        IF(abs(test).GT.1e-15_r8) WRITE(*,*) 'LARGE ERROR: ',irho, is, test
     ENDDO
     is=ns
     test = (rate_ei(irho,is-1) * na(is-1) - rate_rc(irho,is) * na(is))/norm
     IF(abs(test).GT.1e-15_r8) WRITE(*,*) 'LARGE ERROR: ',irho, is, test


  ENDDO
  CLOSE(10)


  RETURN


END SUBROUTINE coronal_distribution


SUBROUTINE solve_tridiag(a,b,c,v,x,n)
  USE itm_types
  IMPLICIT NONE
!      a - sub-diagonal (means it is the diagonal below the main diagonal)
!      b - the main diagonal
!      c - sup-diagonal (means it is the diagonal above the main diagonal)
!      v - right part
!      x - the answer
!      n - number of equations
 
  INTEGER,INTENT(in) :: n
  REAL(kind=R8),DIMENSION(n),INTENT(in) :: a,b,c,v
  REAL(kind=R8),DIMENSION(n),INTENT(out) :: x
  REAL(kind=R8),DIMENSION(n) :: bp,vp
  REAL(kind=R8) :: m
  INTEGER i
 
! Make copies of the b and v variables so that they are unaltered by this sub
  bp(1) = b(1)
  vp(1) = v(1)
  
  !The first pass (setting coefficients):
  firstpass: DO i = 2,n
     m = a(i)/bp(i-1)
     bp(i) = b(i) - m*c(i-1)
     vp(i) = v(i) - m*vp(i-1)
  END DO firstpass
  
  x(n) = vp(n)/bp(n)
  !The second pass (back-substition)
  backsub:DO i = n-1, 1, -1
     x(i) = (vp(i) - c(i)*x(i+1))/bp(i)
  END DO backsub
  
END SUBROUTINE solve_tridiag
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +  
! + + + + + + + + + + + + + + + + + + + + + + + + + + + +   



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

END MODULE coronal