interface_to_amns.F90

Go to the documentation of this file.

module interface_to_amns

  use amns_types
  type (amns_cpo_list), pointer :: first => null()

contains

! zn:                   nuclear charge
! za:                   ionization state
! code                  reaction type
! grid                  output data (has a structure of the type "grid")
! inpdata               input data: amns cpo
! nproc:                number of resulting states
! nz:                   number of initial states (if bundle: number of initial states)
! tdim(2,nproc):        dimesionality (1,nproc) and start positions (2,nproc) of tables (e.g. 1D: only T dependence, 2D: T and n dependence, ...)

  subroutine get_amns_data(reaction_type,reactants,grid,source,shot,run,backend,user,ds_version,ierr,error_description,debug)
    use itm_types
    use amns_utility
    use f90_kind
! do not remove the comment!!!
    use euitm_schemas   ! IGNORE
    use euitm_routines  ! IGNORE
    use read_structures ! IGNORE
    implicit none

    type (grid_t), intent(out) :: grid
    character*(answer_length), intent(out) :: source
    type(amns_reaction_type), intent(in) :: reaction_type
    type(amns_reactants_type), intent(in) :: reactants
    real (kind=r8) :: zn, za
    character (len=*) :: backend, user, ds_version
    integer :: ierr
    character*128 :: error_description
    logical :: debug
    integer :: iza, nz, nproc, i, k, j, isearch, ndim, ncoord
    type (amns_cpo_list), pointer :: last, work => null()
! for calling AMNS CPO
    character(len=5)::      treename
    integer ::              idx, shot, run
    logical :: found, done, match
    character (len=256) :: file
    integer :: nreac, nprod, nlhs, nrhs, ireac
    real (kind=r8) :: zn_cpo, za_cpo
    ierr=0
    error_description=''
    iza = -1

    work => first
    last => first
    found = .false.
!!    write(*,*) associated(first), associated(last), associated(work)
    do while(associated(work))
!!       write(*,*) work%shot, shot
       if(work%shot.eq.shot.and.work%run.eq.run) then
          found=.true.
          exit
       endif
       last => work
       work => work%next
!!       write(*,*) associated(first), associated(last), associated(work)
    enddo
!!    write(*,*) associated(first), associated(last), associated(work)
    if(.not.found) then
       if(.not.associated(first)) then
          allocate(first)
!!       write(*,*) first%shot, first%run, associated(first%prev), associated(first%next)
          work => first
!!       write(*,*) associated(first), associated(work), work%shot
       else if(.not.associated(last%next)) then
          allocate(last%next)
          work => last%next
          work%prev => last
       endif
       if(work%shot.eq.0) then
          treename = 'euitm'
          if(backend.eq.'ascii') then
             write(file,'(a,''_'',i6.6,''_'',i6.6,''.CPO'')') 'amns',shot,run
             if(debug) write(*,*) 'Reading data from ', trim(file), ' for ', shot, run
             call open_read_file(1, trim(file))
             call read_cpo(work%amns_cpo, 'amns')
             call close_read_file
          else
#ifdef UAL
             if(backend.eq.'mdsplus') then
                if(debug) write(*,*) 'Using mdsplus backend'
                call euitm_open_env(treename,shot,run,idx,trim(user),'amns',ds_version)
             elseif(backend.eq.'hdf5') then
                if(debug) write(*,*) 'Using hdf5 backend'
                call euitm_open_hdf5(treename,shot,run,idx)
             else
                error_description = 'Backend not yet supported: ' // trim(backend)
                write(*,*) error_description
                ierr=1
                return
             endif
             call euitm_get(idx,'amns',work%amns_cpo)
             call euitm_close(idx)
#else
             error_description = 'The AMNS routines were compiled without UAL support. ' // &
                  'Use the ascii backend'
             write(*,*) error_description
             ierr=2
             return
#endif
          endif
          work%shot=shot
          work%run=run
       endif
    endif

    nproc= SIZE(work%amns_cpo(1)%process)

    k=0
    do isearch = 1, nproc
! first look for the right process
       if(reaction_type%string .ne. work%amns_cpo(1)%process(isearch)%proc_label(1)) cycle
! found the right process, now search for the right reaction pattern
! first calculate the number of reactants and products in the database
       if(associated(work%amns_cpo(1)%process(isearch)%reactant)) then
          nreac = size(work%amns_cpo(1)%process(isearch)%reactant)
       else
          nreac = 0
       endif
       if(associated(work%amns_cpo(1)%process(isearch)%product)) then
          nprod = size(work%amns_cpo(1)%process(isearch)%product)
       else
          nprod = 0
       endif
! then calculate the number of reactants and products for the requested reaction
       nlhs=0
       nrhs=0
       if(allocated(reactants%components)) then
          do ireac = 1, size(reactants%components)
             if(reactants%components(ireac)%lr .eq. 0) then
                if (nrhs .ne. 0) then
                   error_description = 'Reactants should precede products'
                   write (*,*) error_description
                   ierr=3
                   return
                endif
                nlhs = nlhs + 1
             else if(reactants%components(ireac)%lr .eq. 1) then
                nrhs = nrhs + 1
             else
                write(error_description,*) 'Unknown lr option ', reactants%components(ireac)%lr, &
                     ' for component ', ireac
                write(*,*) error_description
                ierr=4
                return
             endif
          enddo
       endif
! the numbers should match
       if (debug) write(*,*) 'nreac/nprod', isearch, nreac, nprod, nlhs, nrhs
       if(nreac .ne. nlhs) cycle
       if(nprod .ne. nrhs) cycle
       match = .true.
! check the reactants
       do ireac = 1, nreac
! now check that the zn's match
          zn_cpo = 0
          if(associated(work%amns_cpo(1)%process(isearch)%reactant(ireac)%constituents)) then
             if(size(work%amns_cpo(1)%process(isearch)%reactant(ireac)%constituents) .eq. 1) then
!                if(work%amns_cpo(1)%process(isearch)%reactant(ireac)%constituents(1)%multiplicity .eq. 1) then
                   zn_cpo = work%amns_cpo(1)%process(isearch)%reactant(ireac)%constituents(1)%zn
!                else
!                   write(*,*) 'Case with multiplicity != 1 not yet coded'
!                   write(*,*) ireac, work%amns_cpo(1)%process(isearch)%reactant(ireac)%constituents(1)%multiplicity 
!                   stop  'ERROR' ! for the moment
!                endif
             else
                error_description = 'Case with number of constituents != 1 not yet coded'
                write(*,*) error_description
                ierr=5
                return
             endif
          else
             error_description = 'Must have at least one constituent'
             write(*,*) error_description
             ierr=6
             return
          endif
          if (debug) write(*,*) 'zn', isearch, ireac, zn_cpo, reactants%components(ireac)%zn
          match = match .and. (nint(reactants%components(ireac)%zn) .eq. nint(zn_cpo))
          if (.not. match) cycle
! now check that the za's match
          za_cpo = 0
          if(work%amns_cpo(1)%process(isearch)%reactant(ireac)%relative .eq. -1) then
! charge is not relevant
             if (debug) write(*,*) 'za', isearch, ireac, 'za not relavant'
          else if(work%amns_cpo(1)%process(isearch)%reactant(ireac)%relative .eq. 0) then
! charge is absolute
             if (debug) write(*,*) 'za', isearch, ireac, 'absolute charge', work%amns_cpo(1)%process(isearch)%reactant(ireac)%za, reactants%components(ireac)%za
             match = match .and. (work%amns_cpo(1)%process(isearch)%reactant(ireac)%za .eq. reactants%components(ireac)%za)
          else if(work%amns_cpo(1)%process(isearch)%reactant(ireac)%relative .eq. 1 ) then
! charge is relative
! ???
             iza = reactants%components(ireac)%za + 1
             if (debug) write(*,*) 'za', isearch, ireac, 'relative charge', work%amns_cpo(1)%process(isearch)%reactant(ireac)%za, reactants%components(ireac)%za
          else
             write(error_description,*) 'Unrecognized value for relative ', work%amns_cpo(1)%process(isearch)%reactant(ireac)%relative
             write(*,*) error_description
             ierr=7
             return
          endif
          if (.not. match) cycle
! now check that the am's match
          if(debug) then
             write(*,*) 'am', isearch, ireac, work%amns_cpo(1)%process(isearch)%reactant(ireac)%amn, reactants%components(ireac)%mi
          endif
          if(nint(work%amns_cpo(1)%process(isearch)%reactant(ireac)%amn) .ne. 0 .and. nint(reactants%components(ireac)%mi) .ne. 0) then
             match = match .and. (nint(work%amns_cpo(1)%process(isearch)%reactant(ireac)%amn) .eq. nint(reactants%components(ireac)%mi))
          endif
          if (.not. match) cycle
       enddo
       if (.not. match) cycle
! check the products
       do ireac = 1, nprod
! now check that the zn's match
          zn_cpo = 0
          if(associated(work%amns_cpo(1)%process(isearch)%product(ireac)%constituents)) then
             if(size(work%amns_cpo(1)%process(isearch)%product(ireac)%constituents) .eq. 1) then
!                if(work%amns_cpo(1)%process(isearch)%product(ireac)%constituents(1)%multiplicity .eq. 1) then
                   zn_cpo = work%amns_cpo(1)%process(isearch)%product(ireac)%constituents(1)%zn
!                else
!                   write(*,*) 'Case with multiplicity != 1 not yet coded'
!                   stop  'ERROR' ! for the moment
!                endif
             else
                error_description = 'Case with number of constituents != 1 not yet coded'
                write(*,*) error_description
                ierr=8
                return
             endif
          else
             error_description = 'Must have at least one constituent'
             write(*,*) error_description
             ierr=9
             return
          endif
          if (debug) write(*,*) 'zn', isearch, ireac, zn_cpo, reactants%components(ireac+nlhs)%zn
          match = match .and. (nint(reactants%components(ireac+nlhs)%zn) .eq. nint(zn_cpo))
! now check that the za's match
          za_cpo = 0
          if(work%amns_cpo(1)%process(isearch)%product(ireac)%relative .eq. -1) then
! charge is not relevant
             if (debug) write(*,*) 'za', isearch, ireac, 'za not relavant'
          else if(work%amns_cpo(1)%process(isearch)%product(ireac)%relative .eq. 0) then
! charge is absolute
             if (debug) write(*,*) 'za', isearch, ireac, 'absolute charge', work%amns_cpo(1)%process(isearch)%product(ireac)%za, reactants%components(ireac+nlhs)%za
             match = match .and. (work%amns_cpo(1)%process(isearch)%product(ireac)%za .eq.  &
                  reactants%components(ireac+nlhs)%za)
          else if(work%amns_cpo(1)%process(isearch)%product(ireac)%relative .eq. 1 ) then
! charge is relative
! ???
             if (debug) write(*,*) 'za', isearch, ireac, 'relative charge', work%amns_cpo(1)%process(isearch)%product(ireac)%za, reactants%components(ireac+nlhs)%za
          else
             write(error_description,*) 'Unrecognized value for relative ', work%amns_cpo(1)%process(isearch)%reactant(ireac)%relative
             write(*,*) error_description
             ierr=10
             return
          endif
          if (.not. match) cycle
! now check that the am's match
          if(debug) write(*,*) 'am', isearch, ireac, work%amns_cpo(1)%process(isearch)%product(ireac)%amn, reactants%components(ireac+nlhs)%mi
          if(nint(work%amns_cpo(1)%process(isearch)%product(ireac)%amn) .ne. 0 .and. nint(reactants%components(ireac+nlhs)%mi) .ne. 0) then
             match = match .and. (nint(work%amns_cpo(1)%process(isearch)%product(ireac)%amn) .eq. nint(reactants%components(ireac+nlhs)%mi))
          endif
          if (.not. match) cycle
       enddo
       if (.not. match) cycle
       match = .true.
       k = isearch
       exit
    enddo
    if(k .eq. 0) then
       write(error_description,*) 'GET_AMNS_DATA: Case ',trim(reaction_type%string),' not yet implemented'
       if(debug) then
          write(*,*) error_description
          do isearch = 1, nproc
             write(*,*) 'Process ', isearch, ' = ' , work%amns_cpo(1)%process(isearch)%proc_label(1)
          enddo
       endif
       ierr=11
       return
    endif

    ndim = work%amns_cpo(1)%tables(k)%ndim         ! dimensionality of the requested matrix
    ncoord=work%amns_cpo(1)%tables(k)%coord_index
    if(iza .eq. -1) iza=1
    if(size(work%amns_cpo(1)%tables(k)%table) .eq. 1) iza=1
    if(iza .gt. size(work%amns_cpo(1)%tables(k)%table)) then
       error_description = 'iza out of bounds'
       write(*,*) error_description
       ierr=12
       return
    endif

! some of the 1d functions require 2d tables
    if(ndim.eq.1) then
       if(.not.associated(work%amns_cpo(1)%tables(k)%table(iza)%table_1d)) then
          if (debug) write(*,*) 'Special treatment for 1d functions requiring a 2d table'
          if(associated(work%amns_cpo(1)%tables(k)%table(iza)%table_2d)) then
             ndim=2
          endif
       endif
    endif

    done=.true.
    select case (ndim)
    case(1)
       if(ncoord.gt.0) then
          call new_grid(grid, &
               work%amns_cpo(1)%tables(k)%table(iza)%table_1d(:), &
               work%amns_cpo(1)%tables_coord(ncoord)%coords(1)%coord(:))
       else
          call new_grid(grid, &
               work%amns_cpo(1)%tables(k)%table(iza)%table_1d(:))
       endif
    case(2)
       if(ncoord.gt.0) then
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_2d(:,:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(1)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(2)%coord(:))
       else
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_2d(:,:))
       endif
    case(3)
       if(ncoord.gt.0) then
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_3d(:,:,:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(1)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(2)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(3)%coord(:))
       else
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_3d(:,:,:))
       endif
    case(4)
       if(ncoord.gt.0) then
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_4d(:,:,:,:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(1)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(2)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(3)%coord(:), &
            work%amns_cpo(1)%tables_coord(ncoord)%coords(4)%coord(:))
       else
          call new_grid(grid, &
            work%amns_cpo(1)%tables(k)%table(iza)%table_4d(:,:,:,:))
       endif
    case(5)
       done=.false.
       write(*) '0D and 5D can be implemented later'
    end select

    if(done) then
       grid%ndim=work%amns_cpo(1)%tables(k)%ndim
       if(work%amns_cpo(1)%tables(k)%result_trans == 0) then
          grid%is_lin= .true.
       elseif(work%amns_cpo(1)%tables(k)%result_trans == 1) then
          grid%is_log= .true.
       else
          grid%interpol_function = work%amns_cpo(1)%tables(k)%result_trans
       endif
       if(associated(work%amns_cpo(1)%tables(k)%state_label)) then
          if(size(work%amns_cpo(1)%tables(k)%state_label).ge.iza) then
             grid%state_label=work%amns_cpo(1)%tables(k)%state_label(iza)
          endif
       else
          grid%state_label=''
       endif
       if(associated(work%amns_cpo(1)%tables(k)%result_label)) then
          grid%result_label=work%amns_cpo(1)%tables(k)%result_label(1)
       else
          grid%result_label=''
       endif
       if(associated(work%amns_cpo(1)%tables(k)%result_unit)) then
          grid%result_unit=work%amns_cpo(1)%tables(k)%result_unit(1)
       else
          grid%result_unit=''
       endif
       if(ncoord.gt.0) then
          do j=1,ndim
             if(work%amns_cpo(1)%tables_coord(ncoord)%coords(j)%transform == 0) then
                grid%axe(j)%is_lin= .true.
             elseif(work%amns_cpo(1)%tables_coord(ncoord)%coords(j)%transform == 1) then
                grid%axe(j)%is_log= .true.
             else
                write(*,*) 'The coordinate transformation is unknown, the process', trim(reaction_type%string)
             endif
             grid%axe(j)%lbound= lbound(grid%axe(j)%x, 1)
          enddo
       endif
    endif

    if(associated(work%amns_cpo(1)%source)) then
       if(size(work%amns_cpo(1)%source).ge.1) then
          source=work%amns_cpo(1)%source(1)
       else
          source="Not specified in the AMNS CPO"
       endif
    else
       source="Not specified in the AMNS CPO"
    endif

  end subroutine get_amns_data

  subroutine end_amns_data

    use deallocate_structures ! IGNORE
    implicit none

    type (amns_cpo_list), pointer :: work, next => null()

    work => first

    do while(associated(work))
       next => work%next
       call deallocate_cpo(work%amns_cpo)
       deallocate(work)
       work => next
    enddo

    first => null()

  end subroutine end_amns_data

end module interface_to_amns