text:
module m_Interp
  use MPI
  use m_TypeDef
  use m_LogInfo
  use m_Decomp2d
  use m_Parameters
  use m_BC_and_Halo
  use m_MeshAndMetries
  implicit none
  private

  integer, dimension(:),allocatable::indxp,indyp,indzp,indxc,indyc,indzc
  real(RK),dimension(:),allocatable::Ratxp,Ratyp,Ratzp,Ratxc,Ratyc,Ratzc
  type(DECOMP_INFO),save:: decomp_uxNew,decomp_uyNew,decomp_uzNew,decomp_prNew

  public:: Interp_init,Interp_uvwp
contains

  !**********************************************************************
  ! Interp_init
  !**********************************************************************
  subroutine Interp_init()
    implicit none

    ! locals
    integer::i,j,k,ic,jc,kc,js,je,intRed,ierror
    integer::xpSet,xpEnd,xpSize,xcSet,xcEnd,xcSize
    integer::zpSet,zpEnd,zpSize,zcSet,zcEnd,zcSize
    real(RK)::rdxOld,rdzOld,xstCoord,xedCoord,zstCoord,zedCoord

    rdxOld=real(nxcOld,RK);     rdzOld=real(nzcOld,RK)
    xstCoord= xpOld(yStart(1)); xedCoord= xpOld(yend(1)+1)
    zstCoord= zpOld(yStart(3)); zedCoord= zpOld(yend(3)+1)

    ! xp grid
    xpSet=1
    do i=1,nxcNew
      if(xpNew(i)<xstCoord .and. xpNew(i+1)>=xstCoord) then
         xpSet=i+1;exit
      endif
    enddo
    xpEnd=nxcNew
    do i=xpSet,nxcNew
      if(xpNew(i)<xedCoord .and. xpNew(i+1)>=xedCoord) then
         xpEnd=i;exit
      endif
    enddo
    xpSize=xpEnd-xpSet+1
    call MPI_REDUCE(xpSize,intRed,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ierror);intRed=intRed/p_col
    if(intRed/=nxcNew .and. nrank==0) call MainLog%CheckForError(ErrT_Abort,"Interp_init","xpSize wrong !!!")

    ! xc grid
    xcSet=1
    do i=1,nxcNew
      if(xcNew(i)<xstCoord .and. xcNew(i+1)>=xstCoord) then
         xcSet=i+1;exit
      endif
    enddo
    xcEnd=nxcNew
    do i=xcSet,nxcNew
      if(xcNew(i)<xedCoord .and. xcNew(i+1)>=xedCoord) then
         xcEnd=i;exit
      endif
    enddo
    xcSize=xcEnd-xcSet+1
    call MPI_REDUCE(xcSize,intRed,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ierror);intRed=intRed/p_col
    if(intRed/=nxcNew .and. nrank==0) call MainLog%CheckForError(ErrT_Abort,"Interp_init","xcSize wrong !!!")

    ! zp grid
    zpSet=1
    do k=1,nzcNew
      if(zpNew(k)<zstCoord .and. zpNew(k+1)>=zstCoord) then
         zpSet=k+1;exit
      endif
    enddo
    zpEnd=nzcNew
    do k=zpSet,nzcNew
      if(zpNew(k)<zedCoord .and. zpNew(k+1)>=zedCoord) then
         zpEnd=k;exit
      endif
    enddo
    zpSize=zpEnd-zpSet+1
    call MPI_REDUCE(zpSize,intRed,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ierror);intRed=intRed/p_row
    if(intRed/=nzcNew .and. nrank==0) call MainLog%CheckForError(ErrT_Abort,"Interp_init","zpSize wrong !!!")

    ! zc grid
    zcSet=1
    do k=1,nzcNew
      if(zcNew(k)<zstCoord .and. zcNew(k+1)>=zstCoord) then
         zcSet=k+1;exit
      endif
    enddo
    zcEnd=nzcNew
    do k=zcSet,nzcNew
      if(zcNew(k)<zedCoord .and. zcNew(k+1)>=zedCoord) then
         zcEnd=k;exit
      endif
    enddo
    zcSize=zcEnd-zcSet+1
    call MPI_REDUCE(zcSize,intRed,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ierror);intRed=intRed/p_row
    if(intRed/=nzcNew .and. nrank==0) call MainLog%CheckForError(ErrT_Abort,"Interp_init","zcSize wrong !!!")

    ! DECOMP_INFO for ux,uy,uz and pr grid respectively.
    decomp_uxNew%yst(1)=xpSet; decomp_uxNew%yst(2)=1;      decomp_uxNew%yst(3)=zcSet
    decomp_uxNew%yen(1)=xpEnd; decomp_uxNew%yen(2)=nycNew; decomp_uxNew%yen(3)=zcEnd
    decomp_uxNew%ysz=decomp_uxNew%yen-decomp_uxNew%yst+1

    decomp_uyNew%yst(1)=xcSet; decomp_uyNew%yst(2)=1;      decomp_uyNew%yst(3)=zcSet
    decomp_uyNew%yen(1)=xcEnd; decomp_uyNew%yen(2)=nycNew; decomp_uyNew%yen(3)=zcEnd
    decomp_uyNew%ysz=decomp_uyNew%yen-decomp_uyNew%yst+1

    decomp_uzNew%yst(1)=xcSet; decomp_uzNew%yst(2)=1;      decomp_uzNew%yst(3)=zpSet
    decomp_uzNew%yen(1)=xcEnd; decomp_uzNew%yen(2)=nycNew; decomp_uzNew%yen(3)=zpEnd
    decomp_uzNew%ysz=decomp_uzNew%yen-decomp_uzNew%yst+1

    decomp_prNew%yst(1)=xcSet; decomp_prNew%yst(2)=1;      decomp_prNew%yst(3)=zcSet
    decomp_prNew%yen(1)=xcEnd; decomp_prNew%yen(2)=nycNew; decomp_prNew%yen(3)=zcEnd
    decomp_prNew%ysz=decomp_prNew%yen-decomp_prNew%yst+1

    ! interpolation coefficient
    allocate(indxp(xpSet:xpEnd),Ratxp(xpSet:xpEnd))
    allocate(indxc(xcSet:xcEnd),Ratxc(xcSet:xcEnd))
    allocate(indzp(zpSet:zpEnd),Ratzp(zpSet:zpEnd))
    allocate(indzc(zcSet:zcEnd),Ratzc(zcSet:zcEnd))
    allocate(indyp(1:nycNew),Ratyp(1:nycNew),indyc(1:nycNew),Ratyc(1:nycNew))
    do i=xpSet,xpEnd
       ic= floor(xpNew(i)*rdxOld)+1; ic=min(ic,yend(1)); ic=max(ic,ystart(1))
       indxp(i)=ic; Ratxp(i)= (xpOld(ic+1)-xpNew(i))*rdxOld
    enddo
    do i=xcSet,xcEnd
       ic= floor(xcNew(i)*rdxOld)+1; ic=min(ic,yend(1)); ic=max(ic,ystart(1))
       if(xcNew(i)<xcOld(ic)) ic=ic-1;
       indxc(i)=ic; Ratxc(i)= (xcOld(ic+1)-xcNew(i))*rdxOld
    enddo
    do k=zpSet,zpEnd
       kc= floor(zpNew(k)*rdzOld)+1; kc=min(kc,yend(3)); kc=max(kc,ystart(3))
       indzp(k)=kc; Ratzp(k)= (zpOld(kc+1)-zpNew(k))*rdzOld
    enddo
    do k=zcSet,zcEnd
       kc= floor(zcNew(k)*rdzOld)+1; kc=min(kc,yend(3)); kc=max(kc,ystart(3))
       if(zcNew(k)<zcOld(kc)) kc=kc-1;
       indzc(k)=kc; Ratzc(k)= (zcOld(kc+1)-zcNew(k))*rdzOld
    enddo

    do j=1,nycNew
      ! if pos%y is within [0,yly), jc will be within [1,nyc]
      js=0; je=nypOld+1
      do
        jc=(js+je)/2
        if(je-js==1) exit
        if(ypNew(j)< ypOld(jc)) then
          je =jc
        else
          js =jc
        endif
      enddo
      if(jc>nycOld) call MainLog%CheckForError(ErrT_Abort,"Interp_init","wrong jc 1")
      indyp(j)=jc; Ratyp(j)=(ypOld(jc+1)-ypNew(j))/(ypOld(jc+1)-ypOld(jc))
    enddo

    do j=1,nycNew
      ! if pos%y is within [0,yly), jc will be within [1,nyc]
      js=0; je=nypOld+1
      do
        jc=(js+je)/2
        if(je-js==1) exit
        if(ycNew(j)< ypOld(jc)) then
          je =jc
        else
          js =jc
        endif
      enddo
      if(jc>nycOld) call MainLog%CheckForError(ErrT_Abort,"Interp_init","wrong jc 2")

      if(ycNew(j)<ycOld(jc)) jc=jc-1
      indyc(j)=jc; Ratyc(j)=(ycOld(jc+1)-ycNew(j))/(ycOld(jc+1)-ycOld(jc))
    enddo

  end subroutine Interp_init

  !**********************************************************************
  ! Interp_uvwp
  !**********************************************************************
  subroutine Interp_uvwp(mbOld,ArrOld)
    implicit none
    type(MatBound),intent(in)::mbOld
    real(RK),dimension(mbOld%xmm:mbOld%xpm,mbOld%ymm:mbOld%ypm,mbOld%zmm:mbOld%zpm),intent(inout)::ArrOld

    ! locals
    character(256)::chFileOld,chFileNew
    integer::idx1,idy1,idz1,idx2,idy2,idz2
    integer(MPI_OFFSET_KIND)::dispOld,dispNew
    integer::ic,jc,kc,id,jd,kd,ierror,fhOld,fhNew
    real(RK),dimension(:,:,:),allocatable::ArrNew
    real(RK)::prx1,pry1,prz1,prx2,pry2,prz2,MeanValue,MeanValueR

    ! open "RestartOld" and "RestartNew"
    write(chFileOld,'(A)')OldRestartName
    write(chFileNew,'(A)')NewRestartName
    call MPI_FILE_OPEN(MPI_COMM_WORLD, chFileNew, MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, fhNew, ierror)
    call MPI_FILE_OPEN(MPI_COMM_WORLD, chFileOld, MPI_MODE_RDONLY,                 MPI_INFO_NULL, fhOld, ierror)
    if(ierror /= 0 .and. nrank==0) then
      call MainLog%CheckForError(ErrT_Abort,"Interp_uvwp","Cannot open RestartOld: "//trim(chFileOld))
    endif
    dispNew=0_MPI_OFFSET_KIND; dispOld=0_MPI_OFFSET_KIND

    ! ux grid  =============
    allocate(ArrNew(decomp_uxNew%yst(1):decomp_uxNew%yen(1),decomp_uxNew%yst(2):decomp_uxNew%yen(2),decomp_uxNew%yst(3):decomp_uxNew%yen(3)))
    asso_ux: associate( uxOld=>ArrOld, uxNew=>ArrNew); uxOld=zero
    call decomp2d_readVar(fhOld,dispOld,uxOld(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
    call SetBC_and_UpdateHalo_ux(uxOld,mbOld)
    do kc=decomp_uxNew%yst(3),decomp_uxNew%yen(3)
      idz1=indzc(kc);  prz1=Ratzc(kc)
      idz2=idz1+1;     prz2=one-prz1
      do jc=decomp_uxNew%yst(2),decomp_uxNew%yen(2)
        idy1=indyc(jc);  pry1=Ratyc(jc)
        idy2=idy1+1;     pry2=one-pry1
        do ic=decomp_uxNew%yst(1),decomp_uxNew%yen(1)
          idx1=indxp(ic);  prx1=Ratxp(ic)
          idx2=idx1+1;     prx2=one-prx1
          uxNew(ic,jc,kc)= uxOld(idx1,idy1,idz1)*prx1*pry1*prz1 +uxOld(idx2,idy1,idz1)*prx2*pry1*prz1 &
                          +uxOld(idx1,idy2,idz1)*prx1*pry2*prz1 +uxOld(idx2,idy2,idz1)*prx2*pry2*prz1 &
                          +uxOld(idx1,idy1,idz2)*prx1*pry1*prz2 +uxOld(idx2,idy1,idz2)*prx2*pry1*prz2 &
                          +uxOld(idx1,idy2,idz2)*prx1*pry2*prz2 +uxOld(idx2,idy2,idz2)*prx2*pry2*prz2
        enddo
      enddo
    enddo
    if(IsUxConst) then
      MeanValue=zero
      do kc=decomp_uxNew%yst(3),decomp_uxNew%yen(3)
        do jc=decomp_uxNew%yst(2),decomp_uxNew%yen(2)
          do ic=decomp_uxNew%yst(1),decomp_uxNew%yen(1)
            MeanValue=MeanValue+uxNew(ic,jc,kc)*(ypNew(jc+1)-ypNew(jc))
          enddo
        enddo
      enddo
      MeanValue=MeanValue/(real(nxcNew,RK))/(real(nzcNew,RK))/yly
      call MPI_ALLREDUCE(MeanValue,MeanValueR,1,real_type,MPI_SUM,MPI_COMM_WORLD,ierror)
      uxNew= (ubulk/MeanValueR)*uxNew
    endif
    call decomp2d_writeVar(fhNew,dispNew,uxNew,decomp_uxNew)
    end associate asso_ux
    deallocate(ArrNew)

    ! uy grid  =============
    allocate(ArrNew(decomp_uyNew%yst(1):decomp_uyNew%yen(1),decomp_uyNew%yst(2):decomp_uyNew%yen(2),decomp_uyNew%yst(3):decomp_uyNew%yen(3)))
    asso_uy: associate( uyOld=>ArrOld, uyNew=>ArrNew); uyOld=zero
    call decomp2d_readVar(fhOld,dispOld,uyOld(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
    call SetBC_and_UpdateHalo_uy(uyOld,mbOld)
    do kc=decomp_uyNew%yst(3),decomp_uyNew%yen(3)
      idz1=indzc(kc);  prz1=Ratzc(kc)
      idz2=idz1+1;     prz2=one-prz1
      do jc=decomp_uyNew%yst(2),decomp_uyNew%yen(2)
        idy1=indyp(jc);  pry1=Ratyp(jc)
        idy2=idy1+1;     pry2=one-pry1
        do ic=decomp_uyNew%yst(1),decomp_uyNew%yen(1)
          idx1=indxc(ic);  prx1=Ratxc(ic)
          idx2=idx1+1;     prx2=one-prx1
          uyNew(ic,jc,kc)= uyOld(idx1,idy1,idz1)*prx1*pry1*prz1 +uyOld(idx2,idy1,idz1)*prx2*pry1*prz1 &
                          +uyOld(idx1,idy2,idz1)*prx1*pry2*prz1 +uyOld(idx2,idy2,idz1)*prx2*pry2*prz1 &
                          +uyOld(idx1,idy1,idz2)*prx1*pry1*prz2 +uyOld(idx2,idy1,idz2)*prx2*pry1*prz2 &
                          +uyOld(idx1,idy2,idz2)*prx1*pry2*prz2 +uyOld(idx2,idy2,idz2)*prx2*pry2*prz2
        enddo
      enddo
    enddo
    call decomp2d_writeVar(fhNew,dispNew,uyNew,decomp_uyNew)
    end associate asso_uy
    deallocate(ArrNew)

    ! uz grid  =============
    allocate(ArrNew(decomp_uzNew%yst(1):decomp_uzNew%yen(1),decomp_uzNew%yst(2):decomp_uzNew%yen(2),decomp_uzNew%yst(3):decomp_uzNew%yen(3)))
    asso_uz: associate( uzOld=>ArrOld, uzNew=>ArrNew); uzOld=zero
    call decomp2d_readVar(fhOld,dispOld,uzOld(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
    call SetBC_and_UpdateHalo_uz(uzOld,mbOld)
    do kc=decomp_uzNew%yst(3),decomp_uzNew%yen(3)
      idz1=indzp(kc);  prz1=Ratzp(kc)
      idz2=idz1+1;     prz2=one-prz1
      do jc=decomp_uzNew%yst(2),decomp_uzNew%yen(2)
        idy1=indyc(jc);  pry1=Ratyc(jc)
        idy2=idy1+1;     pry2=one-pry1
        do ic=decomp_uzNew%yst(1),decomp_uzNew%yen(1)
          idx1=indxc(ic);  prx1=Ratxc(ic)
          idx2=idx1+1;     prx2=one-prx1
          uzNew(ic,jc,kc)= uzOld(idx1,idy1,idz1)*prx1*pry1*prz1 +uzOld(idx2,idy1,idz1)*prx2*pry1*prz1 &
                          +uzOld(idx1,idy2,idz1)*prx1*pry2*prz1 +uzOld(idx2,idy2,idz1)*prx2*pry2*prz1 &
                          +uzOld(idx1,idy1,idz2)*prx1*pry1*prz2 +uzOld(idx2,idy1,idz2)*prx2*pry1*prz2 &
                          +uzOld(idx1,idy2,idz2)*prx1*pry2*prz2 +uzOld(idx2,idy2,idz2)*prx2*pry2*prz2
        enddo
      enddo
    enddo
    call decomp2d_writeVar(fhNew,dispNew,uzNew,decomp_uzNew)
    end associate asso_uz
    deallocate(ArrNew)

    ! pr grid  =============
    allocate(ArrNew(decomp_prNew%yst(1):decomp_prNew%yen(1),decomp_prNew%yst(2):decomp_prNew%yen(2),decomp_prNew%yst(3):decomp_prNew%yen(3)))
    asso_pr: associate( prOld=>ArrOld, prNew=>ArrNew); prOld=zero
    call decomp2d_readVar(fhOld,dispOld,prOld(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
    call SetBC_and_UpdateHalo_pr(prOld,mbOld)
    do kc=decomp_prNew%yst(3),decomp_prNew%yen(3)
      idz1=indzc(kc);  prz1=Ratzc(kc)
      idz2=idz1+1;     prz2=one-prz1
      do jc=decomp_prNew%yst(2),decomp_prNew%yen(2)
        idy1=indyc(jc);  pry1=Ratyc(jc)
        idy2=idy1+1;     pry2=one-pry1
        do ic=decomp_prNew%yst(1),decomp_prNew%yen(1)
          idx1=indxc(ic);  prx1=Ratxc(ic)
          idx2=idx1+1;     prx2=one-prx1
          prNew(ic,jc,kc)= prOld(idx1,idy1,idz1)*prx1*pry1*prz1 +prOld(idx2,idy1,idz1)*prx2*pry1*prz1 &
                          +prOld(idx1,idy2,idz1)*prx1*pry2*prz1 +prOld(idx2,idy2,idz1)*prx2*pry2*prz1 &
                          +prOld(idx1,idy1,idz2)*prx1*pry1*prz2 +prOld(idx2,idy1,idz2)*prx2*pry1*prz2 &
                          +prOld(idx1,idy2,idz2)*prx1*pry2*prz2 +prOld(idx2,idy2,idz2)*prx2*pry2*prz2
        enddo
      enddo
    enddo
    call decomp2d_writeVar(fhNew,dispNew,prNew,decomp_prNew)
    end associate asso_pr
    deallocate(ArrNew)

#ifdef ScalarFlow
    ! saclar grid  =============
    allocate(ArrNew(decomp_prNew%yst(1):decomp_prNew%yen(1),decomp_prNew%yst(2):decomp_prNew%yen(2),decomp_prNew%yst(3):decomp_prNew%yen(3)))
    asso_scalar: associate( scalarOld=>ArrOld, scalarNew=>ArrNew); scalarOld=zero
    call decomp2d_readVar(fhOld,dispOld,scalarOld(ystart(1):yend(1),ystart(2):yend(2),ystart(3):yend(3)))
    call SetBC_and_UpdateHalo_scalar(scalarOld,mbOld)
    do kc=decomp_prNew%yst(3),decomp_prNew%yen(3)
      idz1=indzc(kc);  prz1=Ratzc(kc)
      idz2=idz1+1;     prz2=one-prz1
      do jc=decomp_prNew%yst(2),decomp_prNew%yen(2)
        idy1=indyc(jc);  pry1=Ratyc(jc)
        idy2=idy1+1;     pry2=one-pry1
        do ic=decomp_prNew%yst(1),decomp_prNew%yen(1)
          idx1=indxc(ic);  prx1=Ratxc(ic)
          idx2=idx1+1;     prx2=one-prx1
          scalarNew(ic,jc,kc)= scalarOld(idx1,idy1,idz1)*prx1*pry1*prz1 +scalarOld(idx2,idy1,idz1)*prx2*pry1*prz1 &
                              +scalarOld(idx1,idy2,idz1)*prx1*pry2*prz1 +scalarOld(idx2,idy2,idz1)*prx2*pry2*prz1 &
                              +scalarOld(idx1,idy1,idz2)*prx1*pry1*prz2 +scalarOld(idx2,idy1,idz2)*prx2*pry1*prz2 &
                              +scalarOld(idx1,idy2,idz2)*prx1*pry2*prz2 +scalarOld(idx2,idy2,idz2)*prx2*pry2*prz2
        enddo
      enddo
    enddo
    if(IsScalarConst) then
      MeanValue=zero
      do kc=decomp_prNew%yst(3),decomp_prNew%yen(3)
        do jc=decomp_prNew%yst(2),decomp_prNew%yen(2)
          do ic=decomp_prNew%yst(1),decomp_prNew%yen(1)
            MeanValue=MeanValue+scalarNew(ic,jc,kc)*(ypNew(jc+1)-ypNew(jc))
          enddo
        enddo
      enddo
      MeanValue=MeanValue/(real(nxcNew,RK))/(real(nzcNew,RK))/yly
      call MPI_ALLREDUCE(MeanValue,MeanValueR,1,real_type,MPI_SUM,MPI_COMM_WORLD,ierror)
      scalarNew= (ScalarMean/MeanValueR)*scalarNew
    endif
    call decomp2d_writeVar(fhNew,dispNew,scalarNew,decomp_prNew)
    end associate asso_scalar
    deallocate(ArrNew)
#endif

    call MPI_FILE_CLOSE(fhOld,ierror)
    call MPI_FILE_CLOSE(fhNew,ierror)
  end subroutine Interp_uvwp

end module m_Interp

max_stars_repo_name: GongZheng-Justin/ParaTC

max_issues_repo_name: GongZheng-Justin/ParaTC

max_forks_repo_name: GongZheng-Justin/ParaTC


text:
      program interarg
      integer s, i
      s = 0
      do i=1, 10
         call addto(s, i)
      enddo
      print *, s
      end
      subroutine addto(x, j)
      integer x, j
      x = x + j
      end
!tps$ echo PROPER REDUCTIONS
!tps$ activate PRINT_CODE_PROPER_REDUCTIONS
!tps$ display PRINTED_FILE($ALL)
!tps$ echo CUMULATED REDUCTIONS
!tps$ activate PRINT_CODE_CUMULATED_REDUCTIONS
!tps$ display PRINTED_FILE($ALL)

ext: f


text:
module rojff_parser_m
    use erloff, only: error_list_t, fatal_t, message_type_t, module_t, procedure_t
    use iso_c_binding, only: c_char, c_double, c_ptr, c_null_char, c_null_ptr
    use rojff_cursor_m, only: cursor_t
    use rojff_fallible_json_value_m, only: &
        fallible_json_value_t, move_into_fallible_json
    use rojff_file_cursor_m, only: file_cursor_t
    use rojff_json_array_m, only: move_into_array
    use rojff_json_bool_m, only: create_json_bool
    use rojff_json_element_m, only: json_element_t
    use rojff_json_integer_m, only: create_json_integer
    use rojff_json_linked_list_m, only: json_linked_list_t
    use rojff_json_member_m, only: json_member_t
    use rojff_json_null_m, only: create_json_null
    use rojff_json_number_m, only: create_json_number
    use rojff_json_object_m, only: move_into_object
    use rojff_json_string_m, only: json_string_t, create_json_string_unsafe
    use rojff_json_value_m, only: json_value_t
    use rojff_member_linked_list_m, only: member_linked_list_t
    use rojff_string_cursor_m, only: string_cursor_t
    use rojff_utils_m, only: to_string

    implicit none
    private
    public :: &
            parse_json_string, &
            parse_json_from_file, &
            parse_json_from_string, &
            INVALID_INPUT

    type(message_type_t), parameter :: INVALID_INPUT = message_type_t( &
            "Invalid Input")
    character(len=*), parameter :: MODULE_NAME = "rojff_parser_m"
contains
    subroutine parse_json(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=*), parameter :: PROCEDURE_NAME = "parse_json"

        call skip_whitespace(cursor)
        if (cursor%finished()) then
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t(PROCEDURE_NAME), &
                    "No json found"))
            return
        end if
        call parse_json_value(cursor, json, errors)
        if (errors%has_any()) then
            errors = error_list_t(errors, module_t(MODULE_NAME), procedure_t(PROCEDURE_NAME))
            return
        end if
        if (.not.cursor%finished()) then
            block
                character(len=:), allocatable :: trailing_content
                allocate(character(len=0) :: trailing_content)
                do while (.not. cursor%finished())
                    trailing_content = trailing_content // cursor%peek()
                    call cursor%next()
                end do
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected trailing content: " // trailing_content))
            end block
        end if
    end subroutine

    recursive subroutine parse_json_value(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=*), parameter :: PROCEDURE_NAME = "parse_json_value"
        character(len=1) :: next_character

        next_character = cursor%peek()

        select case (next_character)
        case ("n")
            call parse_json_null(cursor, json, errors)
        case ("t")
            call parse_json_true(cursor, json, errors)
        case ("f")
            call parse_json_false(cursor, json, errors)
        case ("+", "-", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9")
            call parse_json_number(cursor, json, errors)
        case ('"')
            call cursor%next()
            call parse_json_string(cursor, json, errors)
        case ('[')
            call cursor%next()
            call parse_json_array(cursor, json, errors)
        case ('{')
            call cursor%next()
            call parse_json_object(cursor, json, errors)
        case default
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t(PROCEDURE_NAME), &
                    "At line " // to_string(cursor%current_line()) &
                    // " and column " // to_string(cursor%current_column()) &
                    // " found " // next_character &
                    // ', but expected one of ", {, [, +, -, 0-9, true, false, or null'))
            return
        end select
        if (errors%has_any()) then
            errors = error_list_t(errors, module_t(MODULE_NAME), procedure_t(PROCEDURE_NAME))
        else
            call skip_whitespace(cursor)
        end if
    end subroutine

    subroutine parse_json_null(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=4) :: null_string
        integer :: starting_line, starting_column, i

        starting_line = cursor%current_line()
        starting_column = cursor%current_column()

        do i = 1, 4
            null_string(i:i) = cursor%peek()
            call cursor%next()
            if (cursor%finished()) exit
        end do
        if (null_string == "null") then
            call create_json_null(json)
        else
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t("parse_json_null"), &
                    "At line " // to_string(starting_line) &
                    // " and column " // to_string(starting_column) &
                    // " found " // null_string(1:min(i, 4)) &
                    // ', but expected null'))
        end if
    end subroutine

    subroutine parse_json_true(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=4) :: true_string
        integer :: starting_line, starting_column, i

        starting_line = cursor%current_line()
        starting_column = cursor%current_column()

        do i = 1, 4
            true_string(i:i) = cursor%peek()
            call cursor%next()
            if (cursor%finished()) exit
        end do
        if (true_string == "true") then
            call create_json_bool(json, .true.)
        else
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t("parse_json_true"), &
                    "At line " // to_string(starting_line) &
                    // " and column " // to_string(starting_column) &
                    // " found " // true_string(1:min(i, 4)) &
                    // ', but expected true'))
        end if
    end subroutine

    subroutine parse_json_false(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=5) :: false_string
        integer :: starting_line, starting_column, i

        starting_line = cursor%current_line()
        starting_column = cursor%current_column()

        do i = 1, 5
            false_string(i:i) = cursor%peek()
            call cursor%next()
            if (cursor%finished()) exit
        end do
        if (false_string == "false") then
            call create_json_bool(json, .false.)
        else
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t("parse_json_false"), &
                    "At line " // to_string(starting_line) &
                    // " and column " // to_string(starting_column) &
                    // " found " // false_string(1:min(i, 4)) &
                    // ', but expected false'))
        end if
    end subroutine

    subroutine parse_json_number(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        interface
            function strtod( str, endptr ) result(d) bind(C, name="strtod" )
                ! <stdlib.h> :: double strtod(const char *str, char **endptr)
                import
                character(kind=c_char,len=1),dimension(*),intent(in) :: str
                type(c_ptr), intent(inout) :: endptr
                real(c_double) :: d
            end function strtod
        end interface

        character(len=*), parameter :: NUMBER_SYMBOLS = "0123456789+-.eE"
        character(len=*), parameter :: PROCEDURE_NAME = "parse_json_number"
        character(len=1) :: next_character
        character(len=:), allocatable :: number_string
        integer :: starting_line, starting_column, stat, number_i, precision, exponent_location
        double precision :: number_d
        type(c_ptr) :: endptr

        starting_line = cursor%current_line()
        starting_column = cursor%current_column()

        allocate(character(len=0) :: number_string)
        do while (.not. cursor%finished())
            next_character = cursor%peek()
            if (index(NUMBER_SYMBOLS, next_character) /= 0) then
                number_string = number_string // next_character
                call cursor%next()
            else
                exit
            end if
        end do
        if ( &
                index(number_string, "e") == 0 &
                .and. index(number_string, "E") == 0 &
                .and. index(number_string, ".") == 0) then
            read(number_string, *, iostat=stat) number_i
            if (stat == 0) then
                call create_json_integer(json, number_i)
            else
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(starting_line) &
                        // " and column " // to_string(starting_column) &
                        // " unable to parse " // number_string &
                        // " as an integer"))
            end if
        else
            endptr = c_null_ptr
            number_d = strtod(number_string // c_null_char, endptr)
            if ((.not. number_d < 0.0d0) .and. (.not. number_d > 0.0d0)) then
                read(number_string, *, iostat=stat) number_d  ! not efficient - might really be 0.0
            else
                stat = 0
            end if
            if (stat == 0) then
                exponent_location = index(number_string, "e")
                if (exponent_location == 0) exponent_location = index(number_string, "E")
                if (exponent_location == 0) then
                    precision = len_trim(number_string) - 1
                else
                    precision = exponent_location - merge(1, 2, index(number_string, ".") == 0)
                end if
                call create_json_number(json, number_d, precision)
            else
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(starting_line) &
                        // " and column " // to_string(starting_column) &
                        // " unable to parse " // number_string &
                        // " as a number"))
            end if
        end if
    end subroutine

    subroutine parse_json_string(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=*), parameter :: PROCEDURE_NAME = "parse_json_string"
        character(len=:), allocatable :: the_string
        character(len=1) :: next_character
        integer :: i

        allocate(character(len=0) :: the_string)
        do
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing string " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            next_character = cursor%peek()
            select case (next_character)
            case ('\')
                the_string = the_string // next_character
                call cursor%next()
                if (cursor%finished()) then
                    errors = error_list_t(fatal_t( &
                            INVALID_INPUT, &
                            module_t(MODULE_NAME), &
                            procedure_t(PROCEDURE_NAME), &
                            "Unexpected end of input after escape character " &
                            // "At line " // to_string(cursor%current_line()) &
                            // " and column " // to_string(cursor%current_column())))
                    return
                end if
                next_character = cursor%peek()
                select case (next_character)
                case ('"', '\', '/', 'b', 'f', 'n', 'r', 't')
                    the_string = the_string // next_character
                case ('u')
                    the_string = the_string // next_character
                    do i = 1, 4
                        call cursor%next()
                        if (cursor%finished()) then
                            errors = error_list_t(fatal_t( &
                                    INVALID_INPUT, &
                                    module_t(MODULE_NAME), &
                                    procedure_t(PROCEDURE_NAME), &
                                    "Unexpected end of input while parsing escaped unicode string " &
                                    // "At line " // to_string(cursor%current_line()) &
                                    // " and column " // to_string(cursor%current_column())))
                            return
                        end if
                        next_character = cursor%peek()
                        if (index("0123456789abcdefABCDEF", next_character) == 0) then
                            errors = error_list_t(fatal_t( &
                                    INVALID_INPUT, &
                                    module_t(MODULE_NAME), &
                                    procedure_t(PROCEDURE_NAME), &
                                    "At line " // to_string(cursor%current_line()) &
                                    // " and column " // to_string(cursor%current_column()) &
                                    // " unexpected character in escaped unicode: found " // next_character &
                                    // ', but expected one of 0-9, a-f, or A-F'))
                            return
                        end if
                        the_string = the_string // next_character
                    end do
                case default
                    errors = error_list_t(fatal_t( &
                            INVALID_INPUT, &
                            module_t(MODULE_NAME), &
                            procedure_t(PROCEDURE_NAME), &
                            "At line " // to_string(cursor%current_line()) &
                            // " and column " // to_string(cursor%current_column()) &
                            // " unexpected escaped character in string: found " // next_character &
                            // ', but expected one of ", \, /, b, f, n, r, t, or u'))
                    return
                end select
            case ('"')
                call cursor%next()
                exit
            case default
                the_string = the_string // next_character
            end select
            call cursor%next()
        end do
        call create_json_string_unsafe(json, the_string)
    end subroutine

    recursive subroutine parse_json_array(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=*), parameter :: PROCEDURE_NAME = "parse_json_array"
        type(json_element_t), allocatable :: elements(:)
        type(json_linked_list_t) :: parsed

        call skip_whitespace(cursor)
        if (cursor%finished()) then
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t(PROCEDURE_NAME), &
                    "Unexpected end of input while parsing array " &
                    // "At line " // to_string(cursor%current_line()) &
                    // " and column " // to_string(cursor%current_column())))
            return
        end if
        if (cursor%peek() == ']') then
            allocate(elements(0))
            call move_into_array(json, elements)
            call cursor%next()
            return
        end if
        do
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing array " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            call parse_json_value(cursor, json, errors)
            if (errors%has_any()) then
                errors = error_list_t(errors, module_t(MODULE_NAME), procedure_t(PROCEDURE_NAME))
                return
            end if
            call parsed%append(json)
            call skip_whitespace(cursor)
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing array " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            select case (cursor%peek())
            case (',')
                call cursor%next()
                call skip_whitespace(cursor)
            case (']')
                call cursor%next()
                exit
            case default
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column()) &
                        // " unexpected character in array: found " // cursor%peek() &
                        // ', but expected one of ",", or "]"'))
                return
            end select
        end do
        call parsed%move_into_elements(elements)
        call move_into_array(json, elements)
    end subroutine

    recursive subroutine parse_json_object(cursor, json, errors)
        class(cursor_t), intent(inout) :: cursor
        class(json_value_t), allocatable, intent(out) :: json
        type(error_list_t), intent(out) :: errors

        character(len=*), parameter :: PROCEDURE_NAME = "parse_json_object"
        type(json_member_t), allocatable :: members(:)
        type(member_linked_list_t) :: parsed
        class(json_value_t), allocatable :: key

        call skip_whitespace(cursor)
        if (cursor%finished()) then
            errors = error_list_t(fatal_t( &
                    INVALID_INPUT, &
                    module_t(MODULE_NAME), &
                    procedure_t(PROCEDURE_NAME), &
                    "Unexpected end of input while parsing object " &
                    // "At line " // to_string(cursor%current_line()) &
                    // " and column " // to_string(cursor%current_column())))
            return
        end if
        if (cursor%peek() == '}') then
            allocate(members(0))
            call move_into_object(json, members)
            call cursor%next()
            return
        end if
        do
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing object " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            if (cursor%peek() /= '"') then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column()) &
                        // " unexpected character in object: found " // cursor%peek() &
                        // ", but expected '""'"))
                return
            else
                call cursor%next()
            end if
            call parse_json_string(cursor, key, errors)
            if (errors%has_any()) then
                errors = error_list_t(errors, module_t(MODULE_NAME), procedure_t(PROCEDURE_NAME))
                return
            end if
            call skip_whitespace(cursor)
            if (cursor%peek() /= ":") then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column()) &
                        // " unexpected character in object: found " // cursor%peek() &
                        // ', but expected ":"'))
                return
            else
                call cursor%next()
            end if
            call skip_whitespace(cursor)
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing object " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            call parse_json_value(cursor, json, errors)
            if (errors%has_any()) then
                errors = error_list_t(errors, module_t(MODULE_NAME), procedure_t(PROCEDURE_NAME))
                return
            end if
            select type (key)
            type is (json_string_t)
                call parsed%append(key%string, json)
            class default
                error stop "somehow parsing a string didn't produce a string?"
            end select
            if (cursor%finished()) then
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "Unexpected end of input while parsing object " &
                        // "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column())))
                return
            end if
            call skip_whitespace(cursor)
            select case (cursor%peek())
            case (',')
                call cursor%next()
                call skip_whitespace(cursor)
            case ('}')
                call cursor%next()
                exit
            case default
                errors = error_list_t(fatal_t( &
                        INVALID_INPUT, &
                        module_t(MODULE_NAME), &
                        procedure_t(PROCEDURE_NAME), &
                        "At line " // to_string(cursor%current_line()) &
                        // " and column " // to_string(cursor%current_column()) &
                        // " unexpected character in array: found " // cursor%peek() &
                        // ', but expected one of ",", or "]"'))
                return
            end select
        end do
        call parsed%move_into_members(members)
        call move_into_object(json, members)
    end subroutine

    function parse_json_from_file(filename) result(fallible_json)
        character(len=*), intent(in) :: filename
        type(fallible_json_value_t) :: fallible_json

        type(file_cursor_t) :: cursor
        class(json_value_t), allocatable :: json
        type(error_list_t) :: errors
        integer :: unit

        open( &
                newunit = unit, &
                file = filename, &
                status="OLD", &
                action="READ", &
                access="STREAM", &
                form="FORMATTED")
        cursor = file_cursor_t(unit)
        call parse_json(cursor, json, errors)
        close(unit)
        if (errors%has_any()) then
            fallible_json = fallible_json_value_t(error_list_t( &
                    errors, &
                    module_t(MODULE_NAME), &
                    procedure_t("parse_json_from_string")))
        else
            call move_into_fallible_json(fallible_json, json)
        end if
    end function

    function parse_json_from_string(string) result(fallible_json)
        character(len=*), intent(in) :: string
        type(fallible_json_value_t) :: fallible_json

        type(string_cursor_t) :: cursor
        class(json_value_t), allocatable :: json
        type(error_list_t) :: errors

        cursor = string_cursor_t(string)
        call parse_json(cursor, json, errors)
        if (errors%has_any()) then
            fallible_json = fallible_json_value_t(error_list_t( &
                    errors, &
                    module_t(MODULE_NAME), &
                    procedure_t("parse_json_from_string")))
        else
            call move_into_fallible_json(fallible_json, json)
        end if
    end function

    subroutine skip_whitespace(cursor)
        class(cursor_t), intent(inout) :: cursor

        character(len=1), parameter :: TAB = char(9)
        character(len=1), parameter :: NEWLINE = char(10)
        character(len=1), parameter :: CARRIAGE_RETURN = char(13)
        character(len=1), parameter :: SPACE = char(32)
        character(len=*), parameter :: WHITESPACE = &
                            TAB // NEWLINE // CARRIAGE_RETURN // SPACE

        do while ((.not.cursor%finished()) .and. index(WHITESPACE, cursor%peek()) /= 0)
            call cursor%next()
        end do
    end subroutine
end module

max_stars_repo_name: everythingfunctional/rojff

max_issues_repo_name: everythingfunctional/rojff

max_forks_repo_name: everythingfunctional/rojff


text:
program problem35
  use euler
  implicit none
  integer :: i,answer=1

  do i=3,1000000,2
     if (is_circular_prime(int8(i))) then
        answer=answer+1
     end if
  end do

  print *, answer

contains
  pure function is_circular_prime(n)
    implicit none
    integer*8, intent(in) :: n
    integer*8             :: temp
    integer               :: i
    logical               :: is_circular_prime

    is_circular_prime=.true.

    temp=n
    do i=1,number_length(n)
       if (.not. is_prime(temp)) then
          is_circular_prime=.false.
          exit
       end if
       temp=rotate(temp)
    end do

  end function is_circular_prime

  pure function rotate(n)
    implicit none
    integer*8, intent(in) :: n
    integer*8             :: rotate
    integer               :: p,d

    p=number_length(n)-1
    d=mod(n,10)
    rotate=n/10 + d*(10**p)

  end function rotate

end program problem35

max_stars_repo_name: jamesmcclain/Dioskouroi

max_issues_repo_name: jamesmcclain/Dioskouroi

max_forks_repo_name: jamesmcclain/Dioskouroi


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Computes local number of quads where ABC conditions are imposed.
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0
!> @param[in] nm number of materials
!> @param[in] sd ploynomial degree vector
!> @param[in] tag_mat labels for material vector
!> @param[in] cs_nnz_loc length of cs_loc
!> @param[in] cs_loclocal connectivity vector
!> @param[in] nn_loc number of local nodes
!> @param[in] local_n_num  local node numeration vector
!> @param[in] i4count i4count(i) equal to  1 means that the i-th local node belongs
!                        to an abc surface
!> @param[out] nel_abl  local number of abc surface

      subroutine GET_DIME_ABC(nm, sd, tag_mat, cs_nnz_loc, cs_loc, &
                           nn_loc, local_n_num, &
                           nel_abc, i4count)

     implicit none

     integer*4 :: nm, cs_nnz_loc, nn_loc, nel_abc
     integer*4 :: im, nn, ie, ned
     integer*4 :: ne1, ne2, ne3, ne4, ic1, ic2, ic3, ic4

     integer*4, dimension(nm) :: tag_mat, sd
     integer*4, dimension(0:cs_nnz_loc) :: cs_loc
     integer*4, dimension(nn_loc) :: local_n_num, i4count

     ned = cs_loc(0) - 1

      do im = 1,nm

         nn = sd(im) +1

         do ie = 1,ned

            if (cs_loc(cs_loc(ie -1) + 0) .eq. tag_mat(im)) then

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(1 -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +(nn -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(1 -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +(nn -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

               ne1 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(1 -1) +1)
               ne2 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +(nn -1) +1)
               ne3 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(nn -1) +1)
               ne4 = cs_loc(cs_loc(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +(1 -1) +1)

               if ((i4count(ne1).ne.0) .and. (i4count(ne2).ne.0) .and. (i4count(ne3).ne.0) .and. (i4count(ne4).ne.0)) then
                  nel_abc = nel_abc +1
               endif

            endif
         enddo
      enddo

      return

      end subroutine GET_DIME_ABC


!> Various definitions and tools for running an NGA2 simulation
module simulation
   use precision,     only: WP
   use geometry,      only: cfg1,grp1,isInGrp1
   use geometry,      only: cfg2,grp2,isInGrp2
   use ensight_class, only: ensight
   use coupler_class, only: coupler
   implicit none
   private

   !> Public declarations
   public :: simulation_init,simulation_run,simulation_final

   !> Ensight postprocessing
   type(ensight) :: ens1,ens2

   !> Give ourselves two work arrays
   real(WP), dimension(:,:,:), allocatable :: U1,U2

   !> Coupler
   type(coupler) :: cpl

contains

   !> Initialization of problem solver
   subroutine simulation_init
      implicit none

      ! Both groups prepare the coupler
      if (isInGrp1.or.isInGrp2) then
         ! Create the coupler
         cpl=coupler(src_grp=grp1,dst_grp=grp2,name='test')
         ! Set the grids
         if (isInGrp1) call cpl%set_src(cfg1)
         if (isInGrp2) call cpl%set_dst(cfg2)
         ! Initialize the metrics
         call cpl%initialize()
      end if

      ! Group1 does its initialization work here
      if (isInGrp1) then

         ! Group1 allocates and initializes the field to transfer
         create_field1: block
            use mathtools, only: twoPi
            integer :: i,j,k
            ! Allocate array
            allocate(U1(cfg1%imino_:cfg1%imaxo_,cfg1%jmino_:cfg1%jmaxo_,cfg1%kmino_:cfg1%kmaxo_))
            ! Initialize to solid body rotation
            do k=cfg1%kmino_,cfg1%kmaxo_
               do j=cfg1%jmino_,cfg1%jmaxo_
                  do i=cfg1%imino_,cfg1%imaxo_
                     U1(i,j,k)=-twoPi*cfg1%ym(j)
                  end do
               end do
            end do
         end block create_field1

         ! Group1 also outputs to Ensight
         ens1=ensight(cfg=cfg1,name='grid1')
         call ens1%add_scalar('U',U1)
         call ens1%write_data(0.0_WP)

      end if

      ! Group2 allocates the field to receive
      if (isInGrp2) then
         allocate(U2(cfg2%imino_:cfg2%imaxo_,cfg2%jmino_:cfg2%jmaxo_,cfg2%kmino_:cfg2%kmaxo_))
         U2=-1.0_WP
      end if

      ! Both groups work on the coupling
      coupling_step: block
         if (isInGrp1) call cpl%push(U1)
         call cpl%transfer()
         if (isInGrp2) call cpl%pull(U2)
      end block coupling_step

      ! Group2 outputs its received data
      if (isInGrp2) then
         ens2=ensight(cfg=cfg2,name='grid2')
         call ens2%add_scalar('U',U2)
         call ens2%add_scalar('overlap',cpl%overlap)
         call ens2%write_data(0.0_WP)
      end if

   end subroutine simulation_init

   !> Perform an NGA2 simulation
   subroutine simulation_run
      implicit none
      ! Nothing here
   end subroutine simulation_run

   !> Finalize the NGA2 simulation
   subroutine simulation_final
      implicit none
      ! Nothing here
   end subroutine simulation_final

end module simulation


text:
!     Calculate the total number of people in the US reported by the 1980 census.

      PROGRAM sum1980
!     declare variable types
      INTEGER spop,tpop

      OPEN(UNIT=1,FILE="../data/1980_census.txt",ACTION="READ")

      tpop=0

      DO i=1,51
         READ(UNIT=1, FMT=*) spop
         tpop=tpop+spop
      ENDDO

      PRINT *, tpop

      END PROGRAM sum1980

max_stars_repo_name: IQSS/workshops

max_issues_repo_name: IQSS/workshops

max_forks_repo_name: IQSS/workshops


text:
      subroutine gridal (imajx,iminx,imajy,iminy,ixlab,iylab,iflg,x,y)
      common /sysplt/ mmajx  ,mmajy  ,mminx  ,mminy  ,mxlab  ,mylab  ,
     1                mflg   ,mtype  ,mxa    ,mya    ,mxb    ,myb    ,
     2                mx     ,my     ,mtypex ,mtypey ,xxa    ,yya    ,
     3                xxb    ,yyb    ,xxc    ,yyc    ,xxd    ,yyd    ,
     4                xfactr ,yfactr ,xadd   ,yadd   ,xx     ,yy     ,
     5                mfmtx(3)       ,mfmty(3)       ,mumx   ,mumy   ,
     6                msizx  ,msizy  ,mxdec  ,mydec  ,mxor   ,mop(19),
     7                mname(19)      ,mxold  ,myold  ,mxmax  ,mymax  ,
     8                mxfac  ,myfac  ,modef  ,mf2er  ,mshftx ,mshfty ,
     9                mmgrx  ,mmgry  ,mmnrx  ,mmnry  ,mfrend ,mfrlst ,
     +                mcrout ,mpair1 ,mpair2 ,msblen ,mflcnt ,mjxmin ,
     1                mjymin ,mjxmax ,mjymax ,mnxsto ,mnysto ,mxxsto ,
     2                mxysto ,mprint ,msybuf(360)    ,mncpw  ,minst  ,
     3                mbufa  ,mbuflu ,mfwa(12)       ,mlwa(12)       ,
     4                mipair ,mbprs(16)      ,mbufl  ,munit  ,mbswap ,
     5                small
c
c non-compact version of gridal
c
c     ray bovet ishft changed to ishfta patch
      dimension       nmaj(2),nmin(2),nlab(2),nflg(2),num(2) ,zza(2) ,
     1                zzb(2) ,zzc(2) ,zzd(2) ,ichars(5)      ,
     2                ifmt(3,2)      ,iz(2)  ,iza(2) ,izb(2) ,imz(2) ,
     3                izdec(2)       ,isiz(2),imajl(2)       ,
     4                iminl(2)       ,itype(2)       ,zz(2)  ,
     5                ishfta(2)       ,izaa(2),izbb(2),kz(4)
c     ray bovet patch to avoid small integers being set to 0
      integer x,y,xx,yy
c
c
c     ray bovet ishft changed to ishfta patch
      equivalence     (xxa,zza(1))   ,(xxb,zzb(1))   ,(xxc,zzc(1))   ,
     1                (xxd,zzd(1))   ,(mfmtx(1),ifmt(1,1)),
     2                (mx,iz(1))     ,(mxa,iza(1))   ,(mxb,izb(1))   ,
     3                (majx,nmaj(1)) ,(minx,nmin(1)) ,(mumx,num(1))  ,
     4                (mxdec,izdec(1))       ,(msizx,isiz(1)),
     5                (mmgrx,imajl(1))       ,(mmnrx,iminl(1))       ,
     6                (mtypex,itype(1))      ,(xx,kz(1))     ,
     7                (xx,zz(1))     ,(mshftx,ishfta(1))
c
c set up variables for loop
c
      nmaj(1) = imajx
      nmaj(2) = imajy
      nmin(1) = iminx
      nmin(2) = iminy
      nlab(1) = ixlab
      nlab(2) = iylab
      nflg(1) = ishift(iflg,-2)-1
      nflg(2) = iand(iflg,3)-1
      izaa(1) = iza(1)
      izaa(2) = iza(2)
      izbb(1) = izb(1)
      izbb(2) = izb(2)
      if (nflg(1).le.0 .and. nflg(2).le.0) go to 101
      xx = x
      yy = y
      call trans
      if (nflg(2) .gt. 0) izaa(1) = mx
      if (nflg(1) .gt. 0) izaa(2) = my
      if (nflg(2) .gt. 0) izbb(1) = mx
      if (nflg(1) .gt. 0) izbb(2) = my
  101 continue
      call optn (4hdpat,65535)
      do 121 i=1,2
c
c i=1 for x axis with ticks in y direction
c i=2 for y axis with ticks in x direction
c
         if (nlab(i)) 121,102,102
  102    continue
c
c ior.ne.0 posibility for x only
c
         ixor = (2-i)*90*mxor
         imaj = max0(nmaj(i),1)
         imin = max0(nmin(i),1)
         begin = iza(i)
         biginc = float(izb(i)-iza(i))/float(imaj)
         smlinc = biginc/float(imin)
         start = zzc(i)
         dif = (zzd(i)-zzc(i))/float(imaj)
         iop = 3-i
c
c iop is the opposit axis to i
c
         idec = izdec(iop)
         if (idec .eq. 0) idec = izaa(iop)-izbb(iop)-655
         if (ixor .eq. i-1) go to 103
c
c labels and axis are orthogonal
c
         icent = isign(1,idec-1)
         go to 104
c
c labels and axis are parallel
c
  103    icent = 0
  104    continue
         if (itype(i) .eq. 0) go to 105
         fact = 10.**imaj
         if (zzc(i) .gt. zzd(i)) fact = 1./fact
         val = zzc(i)/fact
         delval = val
         if (imin.le.10 .and. imaj.eq.1) imin = 9
         if (imin .ne. 9) imin = 1
         imaj = abs(alog10(zzd(i)/zzc(i)))+1.0001
  105    imajp1 = imaj+1
         iminm1 = imin-1
         do 119 j=1,imajp1
            part = j-1
c
c draw major line or tick
c
            call optn (4hintn,4hhigh)
            if (itype(i) .ne. 0) go to 106
            iz(i) = begin+part*biginc
            go to 107
  106       val = val*fact
            zz(i) = val
            kz(iop) = 1
            call trans
            delval = delval*fact
            if (iz(i)-10 .gt. izb(i)) go to 120
  107       continue
            iz(iop) = izaa(iop)
            minst = 0
            call put42
            if (nflg(i)) 108,109,109
  108       iz(iop) = izb(iop)
            minst = 1
            call put42
            go to 111
  109       iz(iop) = izaa(iop)+imajl(iop)
            minst = 1
            call put42
            if (nflg(i)) 110,110,111
  110       iz(iop) = izb(iop)
            minst = 0
            call put42
            iz(iop) = izb(iop)-imajl(iop)
            minst = 1
            call put42
  111       continue
c
c form label if needed
c
            if (nlab(i) .le. 0) go to 112
            if (itype(i) .eq. 0) val = start+part*dif
            call encode (num(i),ifmt(1,i),ichars,val)
c     ray bovet ishft changed to ishfta patch
            imz(i) = ishift(iz(i),-ishfta(i))
            imz(iop) = max0(1,ishift(izaa(iop)-idec,-ishfta(iop)))
            njust = num(i)
            if (icent .eq. 0) call justfy (ichars,num(i),njust)
            call pwrit (imz(1),imz(2),ichars,njust,isiz(i),ixor,icent)
c
c put in minor ticks
c
  112       if (iminm1.le.0 .or. j.eq.imajp1) go to 119
            call optn (4hintn,3hlow)
            do 118 k=1,iminm1
               if (itype(i) .ne. 0) go to 113
               iz(i) = begin+part*biginc+float(k)*smlinc
               go to 114
  113          zz(i) = val+float(k)*delval
               if (zzc(i) .gt. zzd(i)) zzi = val-float(k)*delval*.1
               kz(iop) = 1
               call trans
               if (iz(i) .gt. izb(i)) go to 120
               if (iz(i) .lt. iza(i)) go to 118
  114          continue
               iz(iop) = izaa(iop)
               minst = 0
               call put42
               if (nflg(i)) 115,116,116
  115          iz(iop) = izb(iop)
               minst = 1
               call put42
               go to 118
  116          iz(iop) = izaa(iop)+iminl(iop)
               minst = 1
               call put42
               if (nflg(i)) 117,117,118
  117          iz(iop) = izb(iop)
               minst = 0
               call put42
               iz(iop) = izb(iop)-iminl(iop)
               minst = 1
               call put42
  118       continue
  119    continue
         call optn (4hintn,4hhigh)
  120    if (nflg(iop) .lt. 0) go to 121
c
c draw axis line
c
         iz(i) = iza(i)
         iz(iop) = izaa(iop)
         minst = 0
         call put42
         iz(i) = izb(i)
         iz(iop) = izaa(iop)
         minst = 1
         call put42
         if (nflg(i) .gt. 0) go to 121
         iz(i) = iza(i)
         iz(iop) = izb(iop)
         minst = 0
         call put42
         iz(i) = izb(i)
         iz(iop) = izb(iop)
         minst = 1
         call put42
  121 continue
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
!Copyright 2017, UT-Battelle, LLC
!
!SPDX-License-Identifier: BSD-3-Clause
!License-Filename: LICENSE
program test_TpetraImportExport
#include "ForTrilinosTpetra_config.hpp"
#include "FortranTestUtilities.h"
  use iso_fortran_env
  use, intrinsic :: iso_c_binding
  use forteuchos
  use fortpetra

  implicit none
  type(TeuchosComm) :: comm
  integer(global_size_type), parameter :: invalid=-1
  character(len=30), parameter :: FILENAME="test_tpetraimport_export.f90"

  SETUP_TEST()

#ifdef HAVE_MPI
  comm = create_TeuchosComm(MPI_COMM_WORLD); FORTRILINOS_CHECK_IERR()
#else
  comm = create_TeuchosComm()
#endif

  ADD_SUBTEST_AND_RUN(TpetraImportExport_Basic)
  ADD_SUBTEST_AND_RUN(TpetraImportExport_GetNeighborsForward)

  !ADD_SUBTEST_AND_RUN(TpetraImport_isLocallyComplete)
  !ADD_SUBTEST_AND_RUN(TpetraImport_createRemoteOnlyImport)

  call comm%release()

  TEARDOWN_TEST()

contains

  ! -----------------------------------Basic --------------------------------- !
  FORTRILINOS_UNIT_TEST(TpetraImportExport_Basic)
    type(TpetraMap) :: src, tgt
    type(TpetraImport) :: importer
    integer(size_type), parameter :: ten=10, five=5
    integer(size_type) :: same, permute, remote, expected_sum, the_sum

    OUT0("Starting TpetraImportExport_Basic!")

    ! create Maps
    src = create_TpetraMap(invalid, ten, comm)
    tgt = create_TpetraMap(invalid, five, comm)

    ! create Import object
    importer = create_TpetraImport(src, tgt)

    same = importer%getNumSameIDs()
    permute = importer%getNumPermuteIDs()
    remote = importer%getNumRemoteIDs()
    the_sum = same + permute + remote
    expected_sum = tgt%getNodeNumElements()
    TEST_EQUALITY(the_sum, expected_sum)

    OUT0("Finished TpetraImportExport_Basic!")

  END_FORTRILINOS_UNIT_TEST(TpetraImportExport_Basic)

  ! ----------------------------- GetNeighborsForward ------------------------ !

  FORTRILINOS_UNIT_TEST(TpetraImportExport_GetNeighborsForward)
    type(TpetraMap) :: src, tgt
    type(TpetraImport) :: importer
    type(TpetraExport) :: exporter
    type(TpetraMultiVector) :: mv_mine, mv_with_neighbors
    type(TpetraMultiVector) :: mine_parent, neigh_parent
    real(scalar_type), allocatable :: val(:)
    real(scalar_type), pointer :: a(:)
    real(scalar_type), parameter :: zero=0
    integer(size_type), parameter :: ten=10, five=5
    integer(size_type) :: num_images, my_image_id, num_local, num_vecs
    integer(size_type) :: tnum, j, n
    integer(local_ordinal_type) :: lclrow
    integer(global_ordinal_type), allocatable :: neighbors(:), cols(:)

    OUT0("Starting TpetraImportExport_GetNeighborsForward!")

    ! import with the importer to duplicate
    ! export with the exporter to add and reduce
    num_images = comm%getSize()
    my_image_id = comm%getRank()

    if (num_images < 2) return

    ! create a Map
    num_local = 1
    num_vecs = 5

    ! my neighbors: my_image_id-1, me, my_image_id+1
    if (my_image_id == 0 .or. my_image_id == num_images-1) then
      allocate(neighbors(2), a(2), val(2))
    else
      allocate(neighbors(3), a(3), val(3))
    end if

    if (my_image_id == 0) then
      neighbors(:) = [my_image_id+1, my_image_id+2]
    else if (my_image_id == num_images - 1) then
      neighbors(:) = [my_image_id, my_image_id+1]
    else
      neighbors(:) = [my_image_id, my_image_id+1, my_image_id+2]
    end if

    ! two maps: one has one entries per node, the other is the 1-D neighbors
    src = create_TpetraMap(invalid, num_local, comm)
    tgt = create_TpetraMap(invalid, neighbors, comm)

    do tnum = 1, 2

      ! for tnum=1, these are contiguously allocated multivectors
      ! for tnum=2, these are non-contiguous views of multivectors
      if (tnum == 1) then
        mv_mine = create_TpetraMultiVector(src, num_vecs)
        mv_with_neighbors = create_TpetraMultiVector(tgt, num_vecs)
      else
        allocate(cols(5))
        cols(:) = [1,7,4,5,6]
        mine_parent = create_TpetraMultiVector(src, num_vecs+2)
        neigh_parent = create_TpetraMultiVector(tgt, num_vecs+2)
        TEST_ASSERT((num_vecs == 5))
        mv_mine = mine_parent%subViewNonConst(cols)
        mv_with_neighbors = neigh_parent%subViewNonConst(cols)
      end if

      ! mv_mine = [my_image_id  my_image_id+num_images ... my_image_id+4*num_images]
      do j = 1, num_vecs
        lclrow = 1
        val(1) = real(my_image_id + j*num_images, kind=scalar_type)
        call mv_mine%replaceLocalValue(lclrow, j, val(1))
      end do

      ! create Import from src to tgt, Export from tgt to src, test them
      importer = create_TpetraImport(src, tgt)
      exporter = create_TpetraExport(tgt, src)

      ! importer testing
      TEST_ASSERT((src%isSameAs(importer%getSourceMap())))
      TEST_ASSERT((tgt%isSameAs(importer%getTargetMap())))

      n = merge(1, 0, my_image_id == 0)
      TEST_EQUALITY(importer%getNumSameIDs(), n)

      n = merge(0, 1, my_image_id == 0)
      TEST_EQUALITY(importer%getNumPermuteIDs(), n)

      n = merge(1, 2, (my_image_id == 0 .or. my_image_id == num_images - 1))
      TEST_EQUALITY(importer%getNumExportIDs(), n)
      TEST_EQUALITY(importer%getNumRemoteIDs(), n)

      ! exporter testing
      TEST_ASSERT(tgt%isSameAs(exporter%getSourceMap()))
      TEST_ASSERT(src%isSameAs(exporter%getTargetMap()))

      n = merge(1, 0, my_image_id == 0)
      TEST_EQUALITY(importer%getNumSameIDs(), n)

      n = merge(0, 1, my_image_id == 0)
      TEST_EQUALITY(exporter%getNumPermuteIDs(), n)

      ! import neighbors, test their proper arrival
      !                     [ 0    n     2n    3n    4n ]
      ! mv_with_neighbors = [...  ....  ....  ....  ....]
      !                     [n-1  2n-1  3n-1  4n-1  5n-1]
      call mv_with_neighbors%doImport(mv_mine, importer, TpetraREPLACE)
      do j = 1, num_vecs
        a => mv_with_neighbors%getData(j)
        if (my_image_id == 0) then
          val(1) = real(my_image_id+j*num_images,kind=scalar_type)
          val(2) = real(j*num_images+1, kind=scalar_type)
        else if (my_image_id == num_images-1) then
          val(1) = real(my_image_id+j*num_images-1, kind=scalar_type)
          val(2) = real(my_image_id+j*num_images, kind=scalar_type)
        else
          val(1) = real(my_image_id+j*num_images-1, kind=scalar_type)
          val(2) = real(my_image_id+j*num_images, kind=scalar_type)
          val(3) = real(my_image_id+j*num_images+1, kind=scalar_type)
        end if
        TEST_FLOATING_ARRAY_EQUALITY(a, val, epsilon(val(1))) ! me
        TEST_FLOATING_ARRAY_EQUALITY(a, val, epsilon(val(1))) ! neighbor
      end do

      ! export values, test
      call mv_mine%putScalar(zero)
      call mv_mine%doExport(mv_with_neighbors, exporter, TpetraADD)
      do j = 1, num_vecs
        a => mv_mine%getData(j)
        if (my_image_id == 0 .or. my_image_id == num_images-1) then
          ! contribution from me and one neighbor: double original value
          val(1) = real(2.0*(my_image_id+j*num_images), kind=scalar_type)
        else
          val(1) = real(3.0*(my_image_id+j*num_images), kind=scalar_type)
        end if
        TEST_FLOATING_EQUALITY(a(1), val(1), epsilon(val(1)))
      end do

      call mv_mine%release()
      call mv_with_neighbors%release()
      call mine_parent%release()
      call neigh_parent%release()
      call importer%release()
      call exporter%release()

    end do

    deallocate(neighbors, val)
    call src%release()
    call tgt%release()

    OUT0("Finished TpetraImportExport_GetNeighborsForward!")

  END_FORTRILINOS_UNIT_TEST(TpetraImportExport_GetNeighborsForward)

  ! --------------------------------- AbsMax --------------------------------- !
  FORTRILINOS_UNIT_TEST(TpetraImportExport_AbsMax)
    type(TpetraMap) :: smap, dmap
    type(TpetraImport) :: importer
    type(TpetraExport) :: exporter
    type(TpetraMultiVector) :: svec, dvec
    integer(size_type) :: num_images
    real(scalar_type), pointer :: a(:)
    integer(local_ordinal_type) :: lclrow
    integer(global_ordinal_type) :: my_only_gid, cols(2)
    integer(size_type), parameter :: one=1
    integer(global_ordinal_type), parameter :: twog=2

    OUT0("Starting TpetraImportExport_AbsMax!")

    ! test ABSMAX CombineMode
    ! test with local and remote entries, as copyAndPermute() and
    ! unpackAndCombine() both need to be tested
    num_images = comm%getSize()

    if (num_images < 2) return;

    ! create a Map
    smap = create_TpetraMap(invalid, one, comm)

    lclrow = 1
    my_only_gid = smap%getGlobalElement(lclrow)
    cols(1) = my_only_gid
    cols(2) = mod(my_only_gid+1, num_images)
    dmap = create_TpetraMap(twog, cols, comm)

    svec = create_TpetraMultiVector(smap, one)
    call svec%putScalar(-1.0_scalar_type)

    dvec = create_TpetraMultiVector(dmap, one)
    call dvec%putScalar(-3.0_scalar_type)

    ! first item of dvec is local (w.r.t. srcVec), while the second is remote
    ! ergo, during the import:
    ! - the first will be over-written (by 1.0) from the source, while
    ! - the second will be "combined", i.e., abs(max(1.0,3.0)) = 3.0 from the dest
    importer = create_TpetraImport(smap, dmap)
    call dvec%doImport(svec, importer, TpetraABSMAX)

    a => dvec%get1dView()
    TEST_FLOATING_EQUALITY(a(1), -1.0_scalar_type, epsilon(a(1)))
    TEST_FLOATING_EQUALITY(a(2), 3.0_scalar_type, epsilon(a(1)))

    OUT0("Finished TpetraImportExport_AbsMax!")

  END_FORTRILINOS_UNIT_TEST(TpetraImportExport_AbsMax)

  ! ----------------------------isLocallyComplete----------------------------- !
  FORTRILINOS_UNIT_TEST(TpetraImport_isLocallyComplete)
    type(TpetraImport) :: Obj
    OUT0("Starting TpetraImport_isLocallyComplete!")

    success = .false.

    !Obj = create_TpetraImport(); TEST_IERR()
    !fresult = Obj%isLocallyComplete(); TEST_IERR()

    !call Obj%release(); TEST_IERR()

    write(*,*) 'TpetraImport_isLocallyComplete: Test not yet implemented'

    OUT0("Finished TpetraImport_isLocallyComplete!")

  END_FORTRILINOS_UNIT_TEST(TpetraImport_isLocallyComplete)

  ! --------------------------createRemoteOnlyImport-------------------------- !
  FORTRILINOS_UNIT_TEST(TpetraImport_createRemoteOnlyImport)
    type(TpetraImport) :: Obj
    type(TpetraMap) :: remotetarget
    OUT0("Starting TpetraImport_createRemoteOnlyImport!")

    success = .false.

    !remotetarget = create_TpetraMap(); TEST_IERR()
    !Obj = create_TpetraImport(); TEST_IERR()
    !fresult = Obj%createRemoteOnlyImport(remotetarget); TEST_IERR()

    !call remotetarget%release(); TEST_IERR()
    !call Obj%release(); TEST_IERR()

    write(*,*) 'TpetraImport_createRemoteOnlyImport: Test not yet implemented'

    OUT0("Finished TpetraImport_createRemoteOnlyImport!")

  END_FORTRILINOS_UNIT_TEST(TpetraImport_createRemoteOnlyImport)

end program test_TpetraImportExport

max_stars_repo_name: sethrj/ForTrilinos

max_issues_repo_name: sethrj/ForTrilinos

max_forks_repo_name: sethrj/ForTrilinos


text:
Tour de Cluck is an annual bicycle tour of Davis chicken coops organized as a fundraiser for Davis Farm to School, which supports local food in a healthy school environment. The fourth annual event will take place on Saturday, May 24, 2014.

The Tour de Cluck Kickoff begins at 9 a.m. at the Davis Farmers Market, US Bicycling Hall of Fame and Central Park.  A Courage to Cluck exhibition, a Fowl Food Fair and a chicken, bicycle and barnyardthemed Art Exhibit and Silent Auction give visitors good reason to linger awhile!

At 10:00am, those ready to head to the coops will gather around one of four pairs of Bicycle Guides which will disperse groups of riders to one of four Coop Loops (neighborhoodcentered clusters of coops) in a bicycle parade.  Coop owners and chicken experts will be at the coops to share their experiences and answer the questions of the chicken curious.  While bikes are obviously strongly encouraged, car travelers are also welcome.  A wonderful extra resource, Chicken Skool Workshops will address all your fowl questions!

The first annual Tour de Cluck with 16 coop stops took place  on May 22, 2010.

 More details

The coops featured will be presented as a selfguided tour in neighborhood clusters (coop loops!) with most having easy access to Davis’s extensive network of bike paths and bike lanes. In addition to the eggcellent Chicken Coop Crawl, The Tour de Cluck will once again celebrate our best local food purveyors at the tasty Fowl Food Fair at the renowned Davis Farmers Market.

Our Community Education series, Chicken Skool, will feature, among other things, backyard poultry experts to inform your chicken decisions, from breed selection and feeding options to coop design and construction.

And of course, it’s never too early to begin practicing your cluck, bawk and cockledoodledo for the nottobemissed Community Kickoff Clucking Exhibition held on the U.S. Bicycling Hall of Fame “back porch”.  We suggest you begin to form your flock now!

To learn more, become an eggsteemed sponsor, volunteer on the Cluckster event team (its fun!), submit silent auction items and more fowl details, visit our website at http://www.tourdecluck.org or contact us at info@tourdecluck.org.

 Tickets

2014 Tickets are currently on sale now at http://www.brownpapertickets.com/event/581755
 Media Coverage
20100421 http://www.davisvoice.com/2010/04/joinlocalcluckstersincoopsolidarity/ Join Local Clucksters in Coop Solidarity on The Davis Voice
20100513 http://www.sacbee.com/2010/05/13/2745862/tourdavischickencoopsontwo.html The Adventure: Tour Davis Chicken Coops on Two Wheels in The Sacramento Bee
20100525 http://www.davisvoice.com/2010/05/bravotourdecluckbravo/ Bravo, Tour de Cluck, Bravo! on The Davis Voice
20100608 http://www.davisvoice.com/2010/07/tourdecluckanewdavistradition/ Tour de Cluck: A New Davis Tradition (includes awesome video) on The Davis Voice
20110211 http://www.slowfoodyolo.com/index.php/archives/218 Seeking Fowl Frivolity? Look No Further! by Slow Food Yolo
20120723 http://www.nytimes.com/2012/07/25/dining/tourdecluckboomletasurveyofchickencoops.html?pagewanted1&srctwr in the New York Times

 Fowl Resources

http://www.sfgate.com/cgibin/article.cgi?f/c/a/2010/04/25/HOK91CMT9H.DTL Poulet Chalets on The San Francisco Chronicle

 Photo Gallery

 201005  Dagon Jones is sad he did not hear about this event sooner so he could show off his backyard chickens coop he made from 100% recycled materials.

20110308 11:20:33 nbsp Are there any sanitation/safety precautions taken in going from coop to coop?  In all seriousness, its entirely possible that this could spread disease. Users/TomGarberson

   Good quesion. Yes, the have lots of safety precautions.  Check out this document: File(TdeC Chicken Health.pdf) .  Users/KemblePope

20130718 15:19:41 nbsp So awesome! Im putting this on my calendar for next May :) Users/jenb

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*
*     -------------------------------------------------------------
*      O R B I T O R B I T
*     -------------------------------------------------------------
*                                                                  *
*     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
*     OF  ORBIT - ORBIT  INTERACTION                               *
*                                                                  *
*                              2                                   *
*     (n l L S  n l L S ::   -a   Sqrt(2k+1)  (2k-1)/(k+1)         *
*       1 1 1 1  2 2 2 2                                           *
*                                                                  *
*           (k-1) (1)(k)   (k-1) (k)(1)(k)(0)       k-1  k+2       *
*        [[[C   * L  ] * [[C   * L  ] ]   ]        r  / r          *
*              1    1        2    2                 <    >         *
*                                                                  *
*                                     ::n l L S  n l L S )         *
*                                        3 3 3 3  4 4 4 4          *
*                                                                  *
*     Written by G. Gaigalas,                                      *
*     Vanderbilt University,  Nashville            February 1997   *
*
      SUBROUTINE ORBITORBIT(L1,L2,L3,L4,KL,AA)
      IMPLICIT DOUBLEPRECISION (A-H,O-Z)
      COMMON/CONSTS/ZERO,TENTH,HALF,ONE,TWO,THREE,FOUR,SEVEN,ELEVEN,EPS
      AA=ZERO
      IF(ITTK(L1,L3,KL).EQ.0)RETURN
      IF(ITTK(L2,L4,KL).EQ.0)RETURN
      AA=RME (L1,L3,KL)
      IF(DABS(AA).LT.EPS)RETURN
      AA=AA*RME (L2,L4,KL)
      IF(DABS(AA).LT.EPS)RETURN
      K=KL-1
      S=DBLE(2*K+1)*DBLE(L1+L3+K+2)*DBLE(L1+L3-K)*DBLE(L1-L3+K+1)
     :*DBLE(L3-L1+K+1)*DBLE(L2+L4+K+2)*DBLE(L2+L4-K)*DBLE(L2-L4+K+1)
     :*DBLE(L4-L2+K+1)
      AA=AA*TWO*DSQRT(S)/DBLE(K*KL)
      RETURN
      END


text:
!
! This is not a parser or lexer bug but a reporting error.
! action.io_implied_do_control should report a hasStride variable.
!
! from NOTE 9.36
!
! To pass this test should output a line with the following:
!
!  R 919:io-implied-do-control: hasStride=true
!
write(lp, fmt='(10F8.2)') (log(A(i)), i = 1, n+9, k), g

end

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
Surge IV is actually a series of temporary buildings (TB 200  TB 203). They are located across the parking lot from the Silo and across the street from Hutchison Hall. The Engineering Deans Office (Scholars Center) and Vet Med Deans deans office are here.

Surge IV is recognizable from its very 70s style covered courtyard and wooden deck.  The cover is by no means waterproof, which makes it rather interesting when it rains since it has all sorts of random slippery deck locations and localized downpours.

Pictures


! ===================================================================
!
!                      Saturation_Curve
!
! ===================================================================
! This subroutine create a grid in the 2-phase region by discretizing P_sat
! from P_tri to P_crit to obtain a table of the saturation properties:
!
!                 P_sat, V_v, V_l, e_v, e_l, T_sat
!
! The spline coefficients are also computed, so for given values of
! saturation pressure, coefficients of spline are ble to express
!                   u(p_sat), v(p_sat) , T(p_sat)
!
! ===================================================================
MODULE saturation
!
     USE def_constants, ONLY: pr, ord_spline, P_cr,P_tri, NNN_TP, NNN_sat_TP,&
&                             rho_tri_L, rho_tri_R, T_tri,e_cr,NNN_TP2, NNN_sat_TP2,&
&                             rho_cr, T_cr,e_tri_L,e_tri_R
     USE def_variables, ONLY: saturP, saturP_sat, vL_psat_spline, vV_psat_spline,&
&                             uL_psat_spline, uV_psat_spline, Tsat_psat_spline,  &
&                             saturP2, saturP_sat2, vL_psat_spline2, vV_psat_spline2,&
&                             uL_psat_spline2, uV_psat_spline2, Tsat_psat_spline2
     USE non_linear_solvers, ONLY: New_Rap3D
     USE properties, ONLY: inter_energy
     USE grid_functions, ONLY: polyfit

     PRIVATE
     PUBLIC  saturation_curve, sat_curve
     CONTAINS
!========================================================================
!
     SUBROUTINE saturation_curve()
!
!========================================================================
        IMPLICIT NONE
!
        INTEGER  :: i, j, exitflag,Niter
!
        REAL(pr)  :: delta, delta_sat,v_v,v_l,Tsat,u_l,u_v,&
 &                   resnorm, guess1,guess2,guess3,in_2
!
        REAL(pr), DIMENSION (NNN_sat_TP) :: vL, vV, uL, uV,saturT
!
!
!
        delta      = (P_cr - P_tri) / (NNN_TP-1)
        saturP     = P_tri + (/(i*delta, i=0, NNN_TP-1)/)
        delta_sat  = (P_cr - P_tri) / (NNN_sat_TP-1)
        saturP_sat = P_tri + (/(i*delta_sat, i=0,NNN_sat_TP-1)/)
!
! points on saturation curve
!
        guess1 = 1_pr/rho_tri_L
        guess2 = 1_pr/rho_tri_R
        guess3 = T_tri
!
        DO i = 1, NNN_sat_TP-1

         CALL New_Rap3D(3,v_l , v_v, Tsat, &
     &   resnorm, Niter, exitflag, saturP_sat(i),in_2, guess1, guess2,guess3)

        IF (resnorm > 1e-7_pr) THEN
        print*, "saturation curve", resnorm, i
        STOP
        ENDIF
!
!        print*, v_l,v_v
!        print*, u_l,u_v
!        print*," "
         CALL inter_energy(Tsat,v_v,u_v)
         CALL inter_energy(Tsat,v_l,u_l)
           vL(i)     = v_l
           vV(i)     = v_v
           uL(i)     = u_l
           uV(i)     = u_v
           saturT(i) = Tsat
!
           guess1 = v_l
           guess2 = v_v
           guess3 = Tsat
        ENDDO
!
        vL(NNN_sat_TP)     = 1_pr/rho_cr
        vV(NNN_sat_TP)     = 1_pr/rho_cr
        uL(NNN_sat_TP)     = e_cr
        uV(NNN_sat_TP)     = e_cr
        saturT(NNN_sat_TP) = T_cr
! intervals on saturation curve
!
        DO i = 1, NNN_TP-1
           j = 3 * (i-1) + 1
             CALL polyfit(vL_psat_spline(:,i), saturP_sat (j:j+3), vL(j:j+3),ord_spline)
             CALL polyfit(vV_psat_spline(:,i), saturP_sat (j:j+3), vV(j:j+3),ord_spline)
             CALL polyfit(uL_psat_spline(:,i), saturP_sat (j:j+3), uL(j:j+3),ord_spline)
             CALL polyfit(Tsat_psat_spline(:,i), saturP_sat(j:j+3),saturT(j:j+3), ord_spline)
             CALL polyfit(uV_psat_spline(:,i), saturP_sat (j:j+3), uV(j:j+3),ord_spline)
        ENDDO
!print*,'-------------------------------------------------------------------------'
!print*,'           construction for saturation curve with 3e order spline        '
!print*,'-------------------------------------------------------------------------'
!
!
     END SUBROUTINE saturation_curve
!
!===================================================================================
!
     SUBROUTINE sat_curve()
!
!===================================================================================
        IMPLICIT NONE
!
        INTEGER  :: i, j, exitflag,Niter
!
        REAL(pr)  :: delta, delta_sat,v_v,v_l,Tsat,u_l,u_v,&
 &                   resnorm, guess1,guess2,guess3,in_2
!
        REAL(pr), DIMENSION (NNN_sat_TP2) :: vL, vV, uL, uV,saturT
!
!
!
        delta      = (P_cr - P_tri) / (NNN_TP2-1)
        saturP2     = P_tri + (/(i*delta, i=0, NNN_TP2-1)/)
        delta_sat  = (P_cr - P_tri) / (NNN_sat_TP2-1)
        saturP_sat2 = P_tri + (/(i*delta_sat, i=0,NNN_sat_TP2-1)/)
!
! points on saturation curve
!
        guess1 = 1_pr/rho_tri_L
        guess2 = 1_pr/rho_tri_R
        guess3 = T_tri
!
        DO i = 2, NNN_sat_TP2-1

         CALL New_Rap3D(3,v_l , v_v, Tsat, &
     &   resnorm, Niter, exitflag, saturP_sat2(i),in_2, guess1, guess2,guess3)

        IF (resnorm > 1e-12_pr) THEN
        print*, "sat  curve", resnorm, i, Niter
        STOP
        ENDIF
!
!        print*, v_l,v_v
!        print*, u_l,u_v
!        print*," "
         CALL inter_energy(Tsat,v_v,u_v)
         CALL inter_energy(Tsat,v_l,u_l)
           vL(i)     = v_l
           vV(i)     = v_v
           uL(i)     = u_l
           uV(i)     = u_v
           saturT(i) = Tsat
!
           guess1 = v_l
           guess2 = v_v
           guess3 = Tsat
        ENDDO
!
        vL(NNN_sat_TP2)     = 1_pr/rho_cr
        vV(NNN_sat_TP2)     = 1_pr/rho_cr
        uL(NNN_sat_TP2)     = e_cr
        uV(NNN_sat_TP2)     = e_cr
        saturT(NNN_sat_TP2) = T_cr
!
        vL(1)     = 1_pr/rho_tri_L
        vV(1)     = 1_pr/rho_tri_R
        uL(1)     = e_tri_L
        uV(1)     = e_tri_R
        saturT(1) = T_tri

! intervals on saturation curve
!
        DO i = 1, NNN_TP2-1
           j = 5 * (i-1) + 1
             CALL polyfit(vL_psat_spline2(:,i),   saturP_sat2 (j:j+5), vL(j:j+5),     5)
             CALL polyfit(vV_psat_spline2(:,i),   saturP_sat2 (j:j+5), vV(j:j+5),     5)
             CALL polyfit(uL_psat_spline2(:,i),   saturP_sat2 (j:j+5), uL(j:j+5),     5)
             CALL polyfit(Tsat_psat_spline2(:,i), saturP_sat2 (j:j+5), saturT(j:j+5), 5)
             CALL polyfit(uV_psat_spline2(:,i),   saturP_sat2 (j:j+5), uV(j:j+5),     5)
        ENDDO
!print*,'-------------------------------------------------------------------------'
!print*,'construction for saturation curve with 5e order spline                   '
!print*,'-------------------------------------------------------------------------'
!
!
     END SUBROUTINE sat_curve
!
!=================================================================================
!
!     SUBROUTINE satprop(mode, psat, Tsat, vvsat, vlsat, uvsat, ulsat)
!
!=================================================================================
!      IMPLICIT NONE
!!
!      INTEGER  :: i, j, j_sat
!!
!      REAL(pr), INTENT(in)  :: psat
!      INTEGER , INTENT(in)  :: mode
!!
!      REAL(pr), INTENT(out) :: Tsat, vvsat,vlsat,uvsat,ulsat
!!
!      REAL(pr) :: delta, pp, ratio, c_out, a_out, x_out, v_in, du_dp_x, dv_dp_x, temp
!      REAL(pr) :: duL_dp, duV_dp, dvL_dp, dvV_dp
!      REAL(pr) :: vL, uL, vV, uV
!
!
!use 3e order spline
!
!      IF (mode == 3) THEN
!!
!       delta = saturP(2) - saturP(1)
!       j_sat = INT((psat  - saturP(1))/delta) + 1
!!
!!##computing derivatives
!!
!       duL_dp = 0_pr;  duV_dp = 0_pr;  dvL_dp = 0_pr;  dvV_dp = 0_pr
!       DO i = 1, ord_spline
!          pp      = psat**(ord_spline - i)
!          duL_dp  = duL_dp + (ord_spline+1-i) * uL_psat_spline(i,j_sat) *pp
!          duV_dp  = duV_dp + (ord_spline+1-i) * uV_psat_spline(i,j_sat) *pp
!          dvL_dp  = dvL_dp + (ord_spline+1-i) * vL_psat_spline(i,j_sat) *pp
!          dvV_dp  = dvV_dp + (ord_spline+1-i) * vV_psat_spline(i,j_sat) *pp
!       ENDDO
!
!##computing saturation quantities
!
!       vL = 0_pr; uL = 0_pr;  vV = 0_pr;  uV = 0_pr; temp = 0_pr
!       DO i = 1, ord_spline+1
!          pp = psat**(i-1)
!          vL   = vL  + vL_psat_spline  (ord_spline+2-i, j_sat) * pp
!          uL   = uL  + uL_psat_spline  (ord_spline+2-i, j_sat) * pp
!          vV   = vV  + vV_psat_spline  (ord_spline+2-i, j_sat) * pp
!          uV   = uV  + uV_psat_spline  (ord_spline+2-i, j_sat) * pp
!          temp = temp+ Tsat_psat_spline(ord_spline+2-i, j_sat) * pp
!       ENDDO
!
!       ratio = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!       du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!       dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!
!       c_out = SQRT((psat  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!
!       a_out = x_out*vV/v_in
!
!##for output saturation quantities
!
!       Tsat  = temp
!       vvsat = vV
!       vlsat = vL
!       uvsat = uV
!       ulsat = uL
!
!use 5e order spline
!
!      ELSEIF (mode == 5) THEN
!
!       delta = saturP2(2) - saturP2(1)
!       j_sat = INT((psat  - saturP2(1))/delta) + 1
!
!##computing derivatives
!
!       duL_dp = 0_pr;  duV_dp = 0_pr;  dvL_dp = 0_pr;  dvV_dp = 0_pr; temp=0_pr
!       DO i = 1, 5
!          pp      = psat**(5 - i)
!          duL_dp  = duL_dp + (5+1-i) * uL_psat_spline2(i,j_sat) *pp
!          duV_dp  = duV_dp + (5+1-i) * uV_psat_spline2(i,j_sat) *pp
!          dvL_dp  = dvL_dp + (5+1-i) * vL_psat_spline2(i,j_sat) *pp
!          dvV_dp  = dvV_dp + (5+1-i) * vV_psat_spline2(i,j_sat) *pp
!       ENDDO
!
!##computing saturation quantities
!
!       vL = 0_pr; uL = 0_pr;  vV = 0_pr;  uV = 0_pr
!       DO i = 1, 6
!          pp = psat**(i-1)
!          vL   = vL  + vL_psat_spline2  (5+2-i, j_sat) * pp
!          uL   = uL  + uL_psat_spline2  (5+2-i, j_sat) * pp
!          vV   = vV  + vV_psat_spline2  (5+2-i, j_sat) * pp
!          uV   = uV  + uV_psat_spline2  (5+2-i, j_sat) * pp
!          temp = temp+ Tsat_psat_spline2(5+2-i, j_sat) * pp
!       ENDDO
!
!       ratio = ((uV - uL)/(vV - vL))   ! (J/kg)/(m3/kg)
!       du_dp_x = x_out * duV_dp + (1_pr - x_out) * duL_dp
!       dv_dp_x = x_out * dvV_dp + (1_pr - x_out) * dvL_dp
!
!       c_out = SQRT((psat  + ratio)/(du_dp_x - ratio * dv_dp_x)) * v_in ! (m/s)
!
!       a_out = x_out*vV/v_in
!
!##for output saturation quantities
!
!       Tsat  = temp
!       vvsat = vV
!       vlsat = vL
!       uvsat = uV
!       ulsat = uL
!
!      ENDIF
!     END SUBROUTINE satprop
!
!
!
END MODULE saturation


subroutine eqfor(bu, bv, bsq, rmag, zmag)

      use stel_kinds, only: dp
      use name0, only: cp5, c1p5, czero
      use name1, only: nmax, nznt
      use scalars, only: ns, dnorm, isigng, ohs, twopi, voli
      use profs, only: pres, phips, vp, iotas
      use scalefac, only: wint
      use extfld, only: bscale
      use current, only: buco, bvco, jv, ju
      use fsqu, only: equif
      use spectra, only: specw
      use mnarray, only: mscale, nscale

      implicit none

      ! TODO: more elegant definition of these BLAS functions
      real(kind=dp) :: ddot

      real(kind=dp), intent(in) :: bu  (ns,1)
      real(kind=dp), intent(in) :: bv  (ns,1)
      real(kind=dp), intent(in) :: bsq (ns,1)
      real(kind=dp), intent(in) :: rmag(ns,0:nmax)
      real(kind=dp), intent(in) :: zmag(ns,0:nmax)

      integer       :: i, js, n
      real(kind=dp) :: betaxis, es, aiotaf, volf
      real(kind=dp) :: t0, t1, t2, t3
      real(kind=dp) :: phi1, chi1

      write(3,5)
  5   format(/,'     S       EQUIF      PHI       CHI     <JTHETA>     IOTA     <JZETA>   <BSUBV>    P''/V''    V''/PHIP   <M>',/)

      ! compute beta at axis before pressure re-normalization
      betaxis =  c1p5*pres(2)/(dnorm * ddot(nznt, bsq(2,1), ns, wint(2), ns) - pres(2)) &
               -  cp5*pres(3)/(dnorm * ddot(nznt, bsq(3,1), ns, wint(2), ns) - pres(3))

      phi1 = czero
      chi1 = czero
      do i = 2,ns
        buco(i)  = dnorm/bscale * ddot(nznt, bu(i,1), ns, wint(2), ns)
        bvco(i)  = dnorm/bscale * ddot(nznt, bv(i,1), ns, wint(2), ns)

        pres(i)  = pres(i)/bscale**2
        phips(i) = phips(i)/bscale
      end do

      do js = 2, ns-1
        es = (js-1)/real(ns-1)

        t0     = cp5*(   vp(js) +    vp(js+1))
        aiotaf = cp5*(iotas(js) + iotas(js+1))

        jv(js) =  isigng*ohs*(buco(js+1) - buco(js))/t0
        ju(js) = -isigng*ohs*(bvco(js+1) - bvco(js))/t0

        t1 = jv(js)*aiotaf
        t2 = ohs*(pres(js+1) - pres(js))/t0
        t3 = cp5*(vp(js+1)/phips(js+1) + vp(js)/phips(js))

        phi1 = phi1 + phips(js)          /ohs
        chi1 = chi1 + phips(js)*iotas(js)/ohs

        equif(js) = (ju(js) - t1 - t2*t3) / (abs(ju(js)) + abs(t1) + abs(t2*t3))

        write(3,30) es, equif(js), &
                    isigng*twopi*phi1, &
                    isigng*twopi*chi1, &
                    ju(js), aiotaf, jv(js), bvco(js), t2, t3, specw(js)
      end do
 30   FORMAT(1P5E10.2,1p1E11.3,1P3E10.2,1P1E11.3,0PF7.3)

      volf = (twopi**2)*sum(vp(1:ns))/real(ns-1)
      write(3,40) voli, volf, betaxis
 40   format(/,'  INITIAL VOLUME  =  ',1pe20.10,'  FINAL VOLUME  =  ',1pe20.10,/, &
               '  BETA ON AXIS    =  ',2x,1pe12.4,/)

      write(3,50)
 50   format(2x,'MAGNETIC AXIS COEFFICIENTS',/,&
                '    n     rmag       zmag',/)
      do n = 0, nmax
        t1 = mscale(0)*nscale(n)

        ! TODO: why is the z coeff inverted ?
        write(3,70) n, t1*rmag(1,n), -t1*zmag(1,n)
      end do
 70   format(i5,1p2e12.4)

      return
end


text:
! single_site_tmat provides an interface for the single site
! scattering routines.

      subroutine single_site_tmat(nrel_rel,n_spin_cant,is,
     >                       n_spin_pola,
     >                       mtasa,rws,
     >                       nrelv,clight,lmax,kkrsz,
     >                       energy,prel,pnrel,
     >                       vr,h,jmt,jws,r_mesh,
     >                       tmat_l,tmat_g,matom,
     >                       zlr,jlr,
     >                       gz,fz,gj,fj,nuz,indz,
     >                       ubr,ubrd,dmat,dmatp,
     >                       r_sph,iprint,istop)

! tmat_l returns the spin-diagonal blocks of the tmatrix in the local frame
!        in the non-relativistic case only. Its value in the
!        relativistic case is undefined.
! tmat_g returns the tmatrix in the global frame of reference.

      implicit none

      include    'atom_param.h'

      character  sname*32
      character  istop*32
      character  idpot*10

      integer    nuzp
      parameter  (nuzp=2)

      integer    nrel_rel     ! nrel_rel .ne. 0 : relativistic calc.
      integer    n_spin_cant
      integer    n_spin_pola
      integer    is
      integer    nrelv
      integer    mtasa
      integer    lmax
      integer    kkrsz
      integer    jmt
      integer    jws
      integer    iprint
      integer    nuz(2*kkrsz)
      integer    indz(nuzp,2*kkrsz)
      integer    kmymax
      integer    j
      integer    info
      integer    ipvt(2*kkrsz)
      integer    isq
      integer    isp

      real*8     vr(iprpts,n_spin_pola),v0
      real*8     clight
      real*8     h
      real*8     r_sph
      real*8     rws
      real*8     r_mesh(iprpts)
      real*8     vrr(iprpts)
      real*8     brr(iprpts)
      real*8     boprr(iprpts,2)
      real*8     soscal

      complex*16 matom(lmax+1,n_spin_pola)
      complex*16 tmat_l(kkrsz*kkrsz,n_spin_cant)
      complex*16 tmat_g(kkrsz*n_spin_cant,kkrsz*n_spin_cant)
      complex*16 energy
      complex*16 prel
      complex*16 psq
      complex*16 pnrel
      complex*16 zlr(iprpts*(iplmax+1),ip_spin_cant)
      complex*16 jlr(iprpts*(iplmax+1),ip_spin_cant)
      complex*16 gz(iprpts,nuzp,2*kkrsz)
      complex*16 fz(iprpts,nuzp,2*kkrsz)
      complex*16 gj(iprpts,nuzp,2*kkrsz)
      complex*16 fj(iprpts,nuzp,2*kkrsz)
      complex*16 wbig(4*kkrsz*kkrsz)
      complex*16 ubr(4)
      complex*16 ubrd(4)
      complex*16 dmat(2*kkrsz,2*kkrsz)
      complex*16 dmatp(2*kkrsz,2*kkrsz)
      complex*16 detl

      parameter (sname='single_site_tmat')
! soscal: spin-orbit scaling:
! soscal=1 : Dirac eq with full spin-orbit coupling
!        0 : S.O. completly scaled out
      parameter (soscal=1.d0)
      idpot='x123456789'
c     if(vr(1,1).gt.-1.d-3) idpot='Vacuum    '
c     if(vr(1,1).gt.-1.d-3) v0=vr(1,1)/r_mesh(1)
      v0=0.d0
c     if(vr(1,1).gt.-1.d-3) write(*,*)idpot
!     write(6,*) 'entering ',sname
      kmymax = kkrsz*n_spin_cant
      if(nrel_rel.eq.0) then
!       -------------------------------------------------
!        non relativistic
!       -------------------------------------------------
        do isp=1,n_spin_cant
          isq=max(is,isp)
c         -------------------------------------------------------------
          if(iprint.ge.1) then
            write(6,*) 'before SINGLE_SCATTERER_NONREL'
            call flush(6)
          end if
          call single_scatterer_nonrel(nrelv,clight,lmax,kkrsz,
     >                     energy,prel,pnrel,
     >                     vr(1,isq),r_mesh,h,jmt,jws,
     >                     tmat_l(1,isp),matom(1,isp),
     >                      zlr(1,isp),jlr(1,isp),
     >                     r_sph,iprint,istop)
c         -------------------------------------------------------------
          if(iprint.ge.1) then
c           ===========================================================
c           write out Single-site t-matrix (local frame)...............
            write(6,'(/,''GETTAU:: t-matrixfor isp='',i5)') isp
c           -----------------------------------------------------------
            call wrtmtx(tmat_l(1,isp),kkrsz,istop)
c           -----------------------------------------------------------
          endif
        enddo

        if(n_spin_cant.eq.2) then
          call trltog(kkrsz,kkrsz,ubr,ubrd,
     >               tmat_l(1,1),tmat_l(1,2),tmat_g)
        else
          call zcopy(kkrsz*kkrsz,tmat_l(1,1),1,tmat_g,1)
        end if

      else
!       -------------------------------------------------
!       relativistic
!       -------------------------------------------------
!       first build the Clebsch - Gordan coefs (should be done better
!       and not here) !!!!
!       write(6,*) sname,1
!       call clebsch  !moved to beginning of zplanint
!       currently the relativistic part works only for ASA:
        if(mtasa.ne.1) then
          write(6,*) '  relativistic calcultions work for ASA only!'
          call fstop(sname)
        end if
!       first we transform the potential from the up-down form
!        to el-mag form
!       write(6,*) sname,2
        do j=1,jws
!         do we need a factor of 1/2 here or not?
          vrr(j) = 0.5d0*(vr(j,1)+vr(j,2))
          if(vr(1,1).gt.-1.d-4)vrr(j)=vrr(j)-1.d-4
          brr(j) = 0.5d0*(vr(j,1)-vr(j,2))
!         vrr(j) =       (vr(j,1)+vr(j,2))
!         brr(j) =       (vr(j,1)-vr(j,2))
          boprr(j,1) = 0.d0
          boprr(j,2) = 0.d0
        end do
!         write(6,*) h,jws,rws
!         write(6,200) (vrr(j),j=1,jws)
!         write(6,200) (brr(j),j=1,jws)
  200   format(4d20.12)
        psq = energy  +energy*energy/(clight*clight)
!       h = 1.200000000000000d-002
!       rws = 2.60957430788315d0
        call single_scatterer_rel(energy,psq,lmax,kmymax,
     >            'x123456789',v0,
     >            vrr,brr,boprr,h,jws,rws,
     >            tmat_g,gz,fz,gj,fj,nuz,indz,1,soscal,
     >            iprint,istop)

!       now we have tinv but we need t:
!       write(6,*) sname,4
!       call zgetrf(kmymax,kmymax,tmat_g,kkrsz,ipvt,info)
!       call zgetri(kmymax,tmat_g,kmymax,ipvt,wbig,kkrsz*kkrsz,info)
        call gjinv(tmat_g,kmymax,kmymax,detl)
!       t now is in the local frame of refference, now we have to
!       transform it into the global frame:
!       write(6,*) sname,5
        call tripmt(dmat,tmat_g,dmatp,kmymax,kmymax,kmymax)

      end if

      if(istop.eq.sname) then
        write(6,*) 'tmat_g:'
        call outmat1(tmat_g,kmymax,kmymax,kmymax,1.d-15,6)
        call gjinv(tmat_g,kmymax,kmymax,detl)
        call replms(wbig,tmat_g,lmax,kmymax)
        write(6,*) 'diag(tmat_lms):'
        do j=1,2*(lmax+1)**2
          write(6,*) j,j,wbig(j)
        end do
        write(6,*) 't^(-1) :'
        call outmat1(tmat_g,kmymax,kmymax,kmymax,1.d-15,6)
        call fstop(sname)
      end if

      end

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST