text:
      subroutine test_mpi3rma_win_attr( comm )

      use ,intrinsic :: iso_c_binding
      use utilities, only: check_fompi_status

#if INCLUDE_TYPE==1
      implicit none

      include 'fompif.h'
#endif
#if INCLUDE_TYPE==2
      use fompi

      implicit none
#endif

      integer, intent(in) :: comm

      character(20) :: create_flavor
      integer :: win
      integer :: rank
      integer :: j
      integer :: mpierr

      character, dimension(6) :: space

      integer(kind=MPI_ADDRESS_KIND) :: win_size

      integer(kind=MPI_ADDRESS_KIND) :: argument
      type(c_ptr) :: ptr
      character, dimension(:), pointer :: base
      integer, dimension(1) :: shape

      logical :: flag

      character(50) :: filename

      write (filename,*) "win_attr.f90"

      call MPI_Comm_rank( comm, rank, mpierr )

      do j=0,2
        select case (j)
          case (0)
            win_size = 6
            call foMPI_Win_create( space(1), win_size, 3, MPI_INFO_NULL, comm, win, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 30, rank)
            write (create_flavor,*) "create"

            call foMPI_Win_get_attr( win, foMPI_WIN_BASE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 34, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Base address retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

!              if ( return_ptr != &space[0] ) {
!                fprintf(stderr, "Error: Base address retrieved by MPI_Win_get_attr is %p instead of %p for MPI_Win_%s.\n", return_ptr, &space[0], create_flavor);
!                exit(EXIT_FAILURE);
!              }

            call foMPI_Win_get_attr( win, foMPI_WIN_SIZE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 47, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win size retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 6 ) then
              write (*,*) "Error: Window size retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", 6, &
                " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_DISP_UNIT, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 61, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win disp retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 3 ) then
              write (*,*) "Error: Window displacement unit retrieved by", &
                " MPI_Win_get_attr is ", argument, " instead of 3 for ", &
                "MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_CREATE_FLAVOR, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 76, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win create flavor retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_FLAVOR_CREATE ) then
              write (*,*) "Error: Window create flavor retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_FLAVOR_CREATE, " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_MODEL, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 91, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win model retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_UNIFIED ) then
              write (*,*) "Error: Window memory model retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_UNIFIED, " for MPI_Win_", create_flavor
              STOP
            endif

          case (1)
            win_size = 15
            call foMPI_Win_allocate( win_size, 5, MPI_INFO_NULL, comm, ptr, win, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 107, rank)
            shape = 15
            call c_f_pointer( ptr, base, shape )
            write (create_flavor,*) "allocate"

            call foMPI_Win_get_attr( win, foMPI_WIN_BASE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 114, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win base retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

!              if ( return_ptr != ptr ) {
!                fprintf(stderr, "Error: Base address retrieved by MPI_Win_get_attr is %p instead of %p for MPI_Win_%s.\n", return_ptr, ptr, create_flavor);
!                exit(EXIT_FAILURE);
!              }

            call foMPI_Win_get_attr( win, foMPI_WIN_SIZE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 128, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win size retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 15 ) then
              write (*,*) "Error: Window size retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", 15, &
                " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_DISP_UNIT, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 143, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win disp retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 5 ) then
              write (*,*) "Error: Window displacement unit retrieved ", &
                "by MPI_Win_get_attr is ", argument, &
                " instead of 5 for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_CREATE_FLAVOR, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 157, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win create flavor retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_FLAVOR_ALLOCATE ) then
              write (*,*) "Error: Window create flavor retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_FLAVOR_ALLOCATE, " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_MODEL, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 171, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win model retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_UNIFIED ) then
              write (*,*) "Error: Window memory model retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_UNIFIED, " for MPI_Win_", create_flavor
              STOP
            endif

          case (2)
            call foMPI_Win_create_dynamic( MPI_INFO_NULL, comm, win, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 186, rank)
            write (create_flavor, *) "create_dynamic"

            call foMPI_Win_get_attr( win, foMPI_WIN_BASE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 190, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win base address retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. MPI_BOTTOM ) then
              write (*,*) "Error: Base address retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                MPI_BOTTOM, " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_SIZE, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 204, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win size retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 0 ) then
              write (*,*) "Error: Window size retrieved by ", &
                "MPI_Win_get_attr is ", argument, &
                " instead of 0 for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_DISP_UNIT, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 218, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win disp retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. 1 ) then
              write (*,*) "Error: Window displacement unit retrieved ", &
                "by MPI_Win_get_attr is ", argument, &
                " instead of 1 for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_CREATE_FLAVOR, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 232, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win flavor retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_FLAVOR_DYNAMIC ) then
              write (*,*) "Error: Window create flavor retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_FLAVOR_DYNAMIC, " for MPI_Win_", create_flavor
              STOP
            endif

            call foMPI_Win_get_attr( win, foMPI_WIN_MODEL, argument, flag, mpierr )
            call check_fompi_status(mpierr, MPI_SUCCESS, filename, 246, rank)
            if (flag .neqv. .true.) then
              write (*,*) "Error: Win model retrieved by ", &
                "MPI_Win_get_attr ended with a false flag.", create_flavor
              STOP
            endif

            if ( argument .ne. foMPI_WIN_UNIFIED ) then
              write (*,*) "Error: Window memory model retrieved by ", &
                "MPI_Win_get_attr is ", argument, " instead of ", &
                foMPI_WIN_UNIFIED, " for MPI_Win_", create_flavor
              STOP
            endif

        end select

        call foMPI_Win_free( win, mpierr )
        call check_fompi_status(mpierr, MPI_SUCCESS, filename, 261, rank)
      enddo

      write (*,*) "No Error"

      end subroutine test_mpi3rma_win_attr

max_stars_repo_name: spcl/foMPI-NA

max_issues_repo_name: spcl/foMPI-NA

max_forks_repo_name: spcl/foMPI-NA


text:
!-------------------------------------------------------------------!
! SUBROUTINE RESTORE                                                !
!-------------------------------------------------------------------!
!
! Users should modify this subroutine in order to have his/her own
! restoration method.
!
! In its current version, the subroutine uses the general nonlinear
! Augmented Lagrangian solver ALGENCAN to solve the orthogonal
! projection problem
!
!
!                min  ||x - x_k||_2^2
!                s.t. g(x) <= gamma_k
!                     l <= x <= u
!
! All the parameters of the subroutine are carefully explained in the
! interface file 'rinterface.f90' and in the document.

subroutine restore(n,x,l,u,me,mi,uevalc,uevaljac,epsfeas,verbose, &
     infeas,flag)

  use userinterface, only: evalc, evaljac
  use algencan_restoration

  implicit none

  ! SCALAR ARGUMENTS
  integer :: flag,me,mi,n
  real(8) :: epsfeas,infeas
  logical :: verbose

  ! ARRAY ARGUMENTS
  real(8) :: l(n),u(n),x(n)

  intent(in   ) :: epsfeas,l,me,mi,n,u,verbose
  intent(out  ) :: flag,infeas
  intent(inout) :: x

  ! EXTERNAL SUBROUTINES
  procedure(evalc)   :: uevalc
  procedure(evaljac) :: uevaljac

  call algrestore(n,x,l,u,me,mi,uevalc,uevaljac,epsfeas,verbose, &
       infeas,flag)

end subroutine restore

max_stars_repo_name: fsobral/fkss

max_issues_repo_name: fsobral/fkss

max_forks_repo_name: fsobral/fkss


text:
program Nqueens
  implicit none

  integer, parameter :: n = 8  ! size of board
  integer :: file = 1, rank = 1, queens = 0
  integer :: i
  logical :: board(n,n) = .false.

  do while (queens < n)
    board(file, rank) = .true.
    if(is_safe(board, file, rank)) then
      queens = queens + 1
      file = 1
      rank = rank + 1
    else
      board(file, rank) = .false.
      file = file + 1
      do while(file > n)
         rank = rank - 1
         if (rank < 1) then
           write(*, "(a,i0)") "No solution for n = ", n
           stop
         end if
         do i = 1, n
           if (board(i, rank)) then
             file = i
             board(file, rank) = .false.
             queens = queens - 1
             file = i + 1
             exit
           end if
         end do
       end do
    end if
  end do

  call Printboard(board)

contains

function is_safe(board, file, rank)
  logical :: is_safe
  logical, intent(in) :: board(:,:)
  integer, intent(in) :: file, rank
  integer :: i, f, r

  is_safe = .true.
  do i = rank-1, 1, -1
    if(board(file, i)) then
      is_safe = .false.
      return
    end if
  end do

  f = file - 1
  r = rank - 1
  do while(f > 0 .and. r > 0)
    if(board(f, r)) then
      is_safe = .false.
      return
    end if
    f = f - 1
    r = r - 1
  end do

  f = file + 1
  r = rank - 1
  do while(f <= n .and. r > 0)
    if(board(f, r)) then
      is_safe = .false.
      return
    end if
    f = f + 1
    r = r - 1
  end do
end function

subroutine Printboard(board)
  logical, intent(in) :: board(:,:)
  character(n*4+1) :: line
  integer :: f, r

  write(*, "(a, i0)") "n = ", n
  line = repeat("+---", n) // "+"
  do r = 1, n
    write(*, "(a)") line
    do f = 1, n
      write(*, "(a)", advance="no") "|"
      if(board(f, r)) then
        write(*, "(a)", advance="no") " Q "
      else if(mod(f+r, 2) == 0) then
        write(*, "(a)", advance="no") "   "
      else
        write(*, "(a)", advance="no") "###"
      end if
    end do
    write(*, "(a)") "|"
  end do
  write(*, "(a)") line
end subroutine
end program

ext: f

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
!-----------------------------------------------------------------------
! Interface file for pplib2.f
      module pplib2_h
      implicit none
!
      interface
         subroutine PPINIT2(idproc,nvp)
         implicit none
         integer :: idproc, nvp
         end subroutine
      end interface
!
      interface
         subroutine PPFNDGRP(locl,kstrt,nvp,idev,ndev)
         implicit none
         integer :: kstrt, nvp, idev, ndev
         integer, dimension(nvp) :: locl
         end subroutine
      end interface
!
      interface
         subroutine PPEXIT
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine PPABORT
         implicit none
         end subroutine
      end interface
!
      interface
         subroutine PWTIMERA(icntrl,time,dtime)
         implicit none
         integer :: icntrl
         real :: time
         double precision :: dtime
         end subroutine
      end interface
!
      interface
         subroutine PPSUM(f,g,nxp)
         implicit none
         integer :: nxp
         real, dimension(nxp) :: f, g
         end subroutine
      end interface
!
      interface
         subroutine PPDSUM(f,g,nxp)
         implicit none
         integer :: nxp
         double precision, dimension(nxp) :: f, g
         end subroutine
      end interface
!
      interface
         subroutine PPIMAX(if,ig,nxp)
         implicit none
         integer :: nxp
         integer, dimension(nxp) :: if, ig
         end subroutine
      end interface
!
      interface
         subroutine PPPNCGUARD2L(scs,scr,kstrt,nvp,nxv)
         implicit none
         integer :: kstrt, nvp, nxv
         real, dimension(nxv) :: scs, scr
         end subroutine
      end interface
!
      interface
         subroutine PPPNAGUARD2L(scs,scr,kstrt,nvp,nxv)
         implicit none
         integer :: kstrt, nvp, nxv
         real, dimension(nxv) :: scs, scr
         end subroutine
      end interface
!
      interface
         subroutine PPPTPOSE(sm,tm,nx,ny,kxp,kyp,kstrt,nvp)
         implicit none
         integer :: nx, ny, kxp, kyp, kstrt, nvp
         complex, dimension(kxp*kyp,nvp) :: sm, tm
         end subroutine
      end interface
!
      interface
         subroutine PPPTPOSEN(sm,tm,nx,ny,kxp,kyp,kstrt,nvp,ndim)
         implicit none
         integer :: nx, ny, kxp, kyp, kstrt, nvp, ndim
         complex, dimension(kxp*ndim*kyp,nvp) :: sm, tm
         end subroutine
      end interface
!
      interface
         subroutine ACSNDREC(stm,idproc,nsize,ntag,mode)
         implicit none
         integer :: idproc, nsize, ntag, mode
         complex, dimension(nsize) :: stm
         end subroutine
      end interface
!
      interface
         subroutine PPPMOVE2(sbufr,sbufl,rbufr,rbufl,ncll,nclr,mcll,mclr&
     &,kstrt,nvp,idimp,nbmax,mx1)
         implicit none
         integer :: kstrt, nvp, idimp, nbmax, mx1
         real, dimension(idimp,nbmax) :: sbufr, sbufl, rbufr, rbufl
         integer, dimension(3,mx1) :: ncll, nclr, mcll, mclr
         end subroutine
      end interface
!
      end module

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


text:
COLONY is a supercomputer cluster run by CSE. COLONY is the largest computer ever built by wiki:WikiPedia:Sun_Microsystems Sun. It will help professor Jim Crutchfield study wiki:WikiPedia:Cellular_automaton Cellular automata.  It has 14 motherboards, and each motherboard has 64 chips.  Each chip contains 64 32bit processors, for a total of 57,344 processors.  Wowzers!

http://cse.ucdavis.edu/ CSE

http://www.news.ucdavis.edu/search/news_detail.lasso?id7440 Press Release

http://research.sun.com/spotlight/20050503_COLONY.html Suns press release

It is kind of funny that it has 146464 processors, not 166464... Users/DanielBrown

True, but its built on a Sun Ultra Enterprise 10k chassis, which has 14 motherboards.  I believe there are actually 2 of these available, only one is currently had room, power, and network available at the moment.Users/BillBroadley

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine spin(xxm,nclu)

c Randomly spin each collider

      implicit none

      integer mnat
      parameter(mnat=100) ! max number of atoms for each collider

      double precision xxm(3,mnat),ppm(3,mnat),tmp1,tmp2,tmp3
C
      logical iangle
      double precision rot(3,3),q0,q1,q2,q3
C
      integer i,j,nclu
      double precision s1,s2,s3,cs1,ss1,cs2,ss2,cs3,ss3,ranm(4)
C
      iangle=.false.
      if (.not.iangle) then
        s1=2.0d0
        s2=2.0d0
C     generate two random numbers(-1,1) with sum of square < 1
        do while(s1.ge.1.0d0)
          ranm(1) = -1.0+2.0*dble(rand())
          ranm(2) = -1.0+2.0*dble(rand())
          s1 = ranm(1)**2 + ranm(2)**2
        enddo
C     generate another two random numbers with sum of square < 1
        do while(s2.ge.1.0d0 .or. s2.eq.0.d0)
          ranm(3) = -1.0+2.0*dble(rand())
          ranm(4) = -1.0+2.0*dble(rand())
          s2 = ranm(3)**2 + ranm(4)**2
        enddo
        q0 = ranm(1)
        q1 = ranm(2)
        q2 = ranm(3)*dsqrt( (1-s1)/s2 )
        q3 = ranm(4)*dsqrt( (1-s1)/s2 )
C     Construction the transformation matrix
        rot(1,1)=q0**2 + q1**2 - q2**2 - q3**2
        rot(1,2)=2.d0*(q1*q2+q0*q3)
        rot(1,3)=2.d0*(q1*q3-q0*q2)
C
        rot(2,1)=2.d0*(q1*q2-q0*q3)
        rot(2,2)=q0**2 - q1**2 + q2**2 - q3**2
        rot(2,3)=2.d0*(q2*q3+q0*q1)
C
        rot(3,1)=2.d0*(q1*q3+q0*q2)
        rot(3,2)=2.d0*(q2*q3-q0*q1)
        rot(3,3)=q0**2 - q1**2 - q2**2 + q3**2
C     Change on angles s1(theta): 0-pi (angle to +z), s2(phi),
C     s3(posi): 0-2pi
      else
      endif

      do i=1,nclu
        tmp1 = rot(1,1)*xxm(1,i)+rot(1,2)*xxm(2,i)+rot(1,3)*xxm(3,i)
        tmp2 = rot(2,1)*xxm(1,i)+rot(2,2)*xxm(2,i)+rot(2,3)*xxm(3,i)
        tmp3 = rot(3,1)*xxm(1,i)+rot(3,2)*xxm(2,i)+rot(3,3)*xxm(3,i)
        xxm(1,i) = tmp1
        xxm(2,i) = tmp2
        xxm(3,i) = tmp3
c        tmp1 = rot(1,1)*ppm(1,i)+rot(1,2)*ppm(2,i)+rot(1,3)*ppm(3,i)
c        tmp2 = rot(2,1)*ppm(1,i)+rot(2,2)*ppm(2,i)+rot(2,3)*ppm(3,i)
c        tmp3 = rot(3,1)*ppm(1,i)+rot(3,2)*ppm(2,i)+rot(3,3)*ppm(3,i)
c        ppm(1,i) = tmp1
c        ppm(2,i) = tmp2
c        ppm(3,i) = tmp3
      enddo

      return
      end

max_stars_repo_name: mobergd/OneDMin

max_issues_repo_name: mobergd/OneDMin

max_forks_repo_name: mobergd/OneDMin


text:
!-------------------------------------------------------------------------------!
!                                                                               !

subroutine season(y,w,f,s,nlocal,seasonpar,amplitudecutoff,minimum)

!  Detrend the time-series by a running average                                 !
!                                                                               !
!  Fit sine and cosine functions to determine the number of annual seasons      !
!  The current version of the code only allows one or two annual seasons        !
!                                                                               !
!  Determine points that define the fitting intervals for local functions       !
!                                                                               !
!  Authors:                                                                     !
!  Per J\"onsson, Malm\"o University, Sweden                                    !
!  e-mail per.jonsson@ts.mah.se                                                 !
!                                                                               !
!  Lars Eklundh, Lund University, Sweden                                        !
!  e-mail lars.eklundh@nateko.lu.se                                             !
!                                                                               !
!  Updated May 2006                                                             !
!                                                                               !
!-------------------------------------------------------------------------------!

implicit none

integer :: i,j,nyear,nptperyear,npt,nseason,nenvi,printflag,debugflag,win(3),m1,m2
integer :: s(4*nyear+2),nlocal,yextpos(4),posmin1(1),posmin2(1),posmax1(1)
integer :: posmax2(1),yminpos,minimum
double precision :: y(npt),w(npt),ydt(npt),yav(npt),yfit(npt),f(npt,5),c(5)
double precision :: wfact,amp,yext(4), seasonpar, amplitudecutoff

common /parameters/ nyear,nptperyear,npt,nenvi,printflag,debugflag,wfact,win

!---- Use running average with large window to remove trend in time-series -----

do i = 1,npt

!---- Set window ---------------------------------------------------------------

  m1 = max(1,i - nptperyear/2)
  m2 = min(npt,i + nptperyear/2)

!---- Compute average in the window --------------------------------------------

  yav(i) = sum(w(m1:m2)*y(m1:m2))/sum(w(m1:m2))

end do

!---- Remove the trend ---------------------------------------------------------

ydt = y - yav

!---- Fit sine and cosine functions to detrended time-series -------------------

call linlsq(npt,5,f,ydt,w,c)

!---- Evaluate the fitted function ---------------------------------------------

yfit = c(1)*f(:,1)+c(2)*f(:,2)+c(3)*f(:,3)+c(4)*f(:,4)+c(5)*f(:,5)

!---- Check if the amplitude is larger than the criteria -----------------------

if ( (maxval(yfit) - minval(yfit)) <= amplitudecutoff ) then
  nlocal = 0
  s = 0
  return
end if

!---- Determine points that define the fitting intervals for local functions ---

j = 0
if ( (yfit(1)-yfit(nptperyear))*(yfit(2)-yfit(1)) < 0.d0 ) then
  j = j + 1
  yext(j) = yfit(1)
  yextpos(j) = 1
end if

do i = 2,nptperyear-1
  if ( (yfit(i)-yfit(i-1))*(yfit(i+1)-yfit(i)) < 0.d0 ) then
    j = j + 1
    yext(j) = yfit(i)
    yextpos(j) = i
  end if
end do

if ( (yfit(nptperyear)-yfit(nptperyear-1))*(yfit(1)-yfit(nptperyear)) < 0.d0 ) then
  j = j + 1
  yext(j) = yfit(nptperyear)
  yextpos(j) = nptperyear
end if

if (j == 2) then
  s(1:2*nyear) = (/ ((i-1)*nptperyear+yextpos(1:2), i = 1,nyear) /)
  s(1:2*nyear+2) = (/ 1,s(1:2*nyear),npt /)
  nlocal = 2*nyear
else
  posmax1(1:1) = maxloc(yext)
  posmax2(1:1) = maxloc(yext,mask=yext /= yext(posmax1(1)))
  posmin1(1:1) = minloc(yext)
  posmin2(1:1) = minloc(yext,mask=yext /= yext(posmin1(1)))

  if ( (yext(posmax2(1)) - yext(posmin2(1)))/(yext(posmax1(1)) - yext(posmin1(1))) > seasonpar ) then

!---- Two seasons ------------------------------------------------------------

    s(1:4*nyear) = (/ ((i-1)*nptperyear+yextpos(1:4), i = 1,nyear) /)
    s(1:4*nyear+2) = (/ 1,s(1:4*nyear),npt /)
    nlocal = 4*nyear

  else

!---- One season -------------------------------------------------------------

!---- Fit single sine and cosine functions to detrended time-series ----------

!   write(*,*) sum(ydt),sum(w)
    call linlsq(npt,3,f(:,1:3),ydt,w,c)

!---- Evaluate the fitted function -------------------------------------------

    yfit = c(1)*f(:,1)+c(2)*f(:,2)+c(3)*f(:,3)

    j = 0
    if ( (yfit(1)-yfit(nptperyear))*(yfit(2)-yfit(1)) < 0.d0 ) then
      j = j + 1
      yext(j) = yfit(1)
      yextpos(j) = 1
    end if

    do i = 2,nptperyear-1
      if ( (yfit(i)-yfit(i-1))*(yfit(i+1)-yfit(i)) < 0.d0 ) then
        j = j + 1
        yext(j) = yfit(i)
        yextpos(j) = i
      end if
    end do

    if ( (yfit(nptperyear)-yfit(nptperyear-1))*(yfit(1)-yfit(nptperyear)) < 0.d0 ) then
      j = j + 1
      yext(j) = yfit(nptperyear)
      yextpos(j) = nptperyear
    end if

	s(1:2*nyear) = (/ ((i-1)*nptperyear+yextpos(1:2), i = 1,nyear) /)
    s(1:2*nyear+2) = (/ 1,s(1:2*nyear),npt /)
    nlocal = 2*nyear

  end if
end if

!---- Find out if position corresponds to a minimum or a maximum. This information is ---
!     used in PHENOLOGY

if ( yfit(s(2)) < yfit(s(3)) ) then
  minimum = 1
else
  minimum = 0
end if

end subroutine season

max_stars_repo_name: bucricket/projectMASpreprocess

max_issues_repo_name: bucricket/projectMASpreprocess

max_forks_repo_name: bucricket/projectMASpreprocess


subroutine save_game(auxiliary_empty_cells,auxiliary_empty_cells_history,current_empty_cells_history,empty_cells,empty_cells_history,matrix,matrix_history,score_history)
	implicit none!all variables must be declared
	logical::file_exists=.false.!if true, the temporary file to store the matrix is present
	character(1000)::filename!the name of the file to save to
	character(1)::auxiliary_input='-'!used in case of overwriting
	integer(4),intent(inout)::auxiliary_empty_cells(0:255)!this vector keeps track of where the index of a certain cell is located in the empty_cells array, for example: empty_cells(0:3)=(/1,3,0,2/) would mean that the values of this array are the following: auxiliary_empty_cells(0:3)=(/2,0,1,3/), thus, only the index of a cell is needed when its interal data are required
	integer(4),intent(inout)::auxiliary_empty_cells_history(0:255,0:99)!same as matrix_history, but the player doesn't interact with this variable, it's only kept so it's faster to move back in the history of movements
	integer(4)::cells!number of cells per side of the square to be used in the game
	integer(4)::current_empty_cells!this is the amount of empty cells in a certain moment
	integer(4),intent(inout)::current_empty_cells_history(0:99)!same as matrix_history, but the player doesn't interact with this variable, it's only kept so it's faster to move back in the history of movements
	integer(4),intent(inout)::empty_cells(0:255)!the indeces of empty cells are sequentially kept here with no occupied cells in between, all of those are at the end of the array
	integer(4),intent(inout)::empty_cells_history(0:255,0:99)!same as matrix_history, but the player doesn't interact with this variable, it's only kept so it's faster to move back in the history of movements
	integer(4)::history_position!this variable keeps track of where should the program write in the history arrays
	integer(4)::i=0!first auxiliary index for the matrix
	integer(4)::j=0!second auxiliary index for the matrix
	integer(4)::k=0!third auxiliary index for the matrix_history array
	integer(4),intent(inout)::matrix(0:15,0:15)!all the number cells are here
	integer(4),intent(inout)::matrix_history(0:15,0:15,0:99)!history of player movements, when stored movements reach 99, it saves the next snapshot of the matrix in 0, overwriting the previous value, and starts overwriting the next ones as more snapshots are generated by the player
	integer(4)::score!this keeps the puntuation of the game
	integer(4),intent(inout)::score_history(0:99)!same as matrix_history, but the player doesn't interact with this variable, it's only kept so it's faster to move back in the history of movements
	integer(4)::snapshots!this is the amount of movements stored in the third index of matrix_history, when it reaches 100 it doesn't increase anymore
	common /basic/ cells,current_empty_cells,history_position,snapshots!sharing the variables needed for this subroutine except score
	common /score/ score!sharing the score variable too

	write(*,*) 'set title "Which file do you want to save to? Introduce its name or path."'!asking the player for the name of the file
	write(*,*) 'replot'!updating the output
10	read(*,*) filename!reading the input
	inquire(file=filename,exist=file_exists)!checking if the file is there
	if(file_exists)then
		write(*,*) 'set title "File already exists, do you want to overwrite it? (Just introduce y to overwrite or anything else to try a different name)"'!making sure nothing is unnecesarily overwritten
		write(*,*) 'replot'!updating the output
		read(*,*) auxiliary_input!reading the input
		if(auxiliary_input=='y')then
			open(20,file=filename,action='write',form='unformatted',access='stream',status='replace')!opening the file and overwriting it
		else
			write(*,*) 'set title "Try a different name."'!making sure nothing is unnecesarily overwritten
			write(*,*) 'replot'!updating the output
			go to 10!cycling through this if loop till a valid input is obtained
		end if
	else
		open(20,file=filename,action='write',form='unformatted',access='stream',status='new')!opening the new file
	end if
	write(20) cells,current_empty_cells,history_position,score,snapshots!reading the scalar variables
	do j=0,15
		do i=0,15
			write(20) matrix(i,j)!reading the main matrix
		end do
	end do
	do i=0,255
		write(20) auxiliary_empty_cells(i),empty_cells(i)!reading both empty cells arrays
	end do
	do k=0,99
		do j=0,15
			do i=0,15
				write(20) matrix_history(i,j,k)!reading the matrix_history array
			end do
		end do
	end do
	do j=0,99
		do i=0,255
			write(20) auxiliary_empty_cells_history(i,j),empty_cells_history(i,j)!reading both empty cells history arrays
		end do
	end do
	do i=0,99
		write(20) current_empty_cells_history(i),score_history(i)!reading the current empty cells and score history arrays
	end do
	close(20)!closing the file to save changes
	return
end subroutine save_game


text:
      subroutine dqk15(f,a,b,result,abserr,resabs,resasc)
c***begin prologue  dqk15
c***date written   800101   (yymmdd)
c***revision date  830518   (yymmdd)
c***category no.  h2a1a2
c***keywords  15-point gauss-kronrod rules
c***author  piessens,robert,appl. math. & progr. div. - k.u.leuven
c           de doncker,elise,appl. math. & progr. div - k.u.leuven
c***purpose  to compute i = integral of f over (a,b), with error
c                           estimate
c                       j = integral of abs(f) over (a,b)
c***description
c
c           integration rules
c           standard fortran subroutine
c           double precision version
c
c           parameters
c            on entry
c              f      - double precision
c                       function subprogram defining the integrand
c                       function f(x). the actual name for f needs to be
c                       declared e x t e r n a l in the calling program.
c
c              a      - double precision
c                       lower limit of integration
c
c              b      - double precision
c                       upper limit of integration
c
c            on return
c              result - double precision
c                       approximation to the integral i
c                       result is computed by applying the 15-point
c                       kronrod rule (resk) obtained by optimal addition
c                       of abscissae to the7-point gauss rule(resg).
c
c              abserr - double precision
c                       estimate of the modulus of the absolute error,
c                       which should not exceed abs(i-result)
c
c              resabs - double precision
c                       approximation to the integral j
c
c              resasc - double precision
c                       approximation to the integral of abs(f-i/(b-a))
c                       over (a,b)
c
c***references  (none)
c***routines called  d1mach
c***end prologue  dqk15
c
      double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1,
     *  d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc,
     *  resg,resk,reskh,result,uflow,wg,wgk,xgk
      integer j,jtw,jtwm1
      external f
c
      dimension fv1(7),fv2(7),wg(4),wgk(8),xgk(8)
c
c           the abscissae and weights are given for the interval (-1,1).
c           because of symmetry only the positive abscissae and their
c           corresponding weights are given.
c
c           xgk    - abscissae of the 15-point kronrod rule
c                    xgk(2), xgk(4), ...  abscissae of the 7-point
c                    gauss rule
c                    xgk(1), xgk(3), ...  abscissae which are optimally
c                    added to the 7-point gauss rule
c
c           wgk    - weights of the 15-point kronrod rule
c
c           wg     - weights of the 7-point gauss rule
c
c
c gauss quadrature weights and kronron quadrature abscissae and weights
c as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
c bell labs, nov. 1981.
c
      data wg  (  1) / 0.1294849661 6886969327 0611432679 082 d0 /
      data wg  (  2) / 0.2797053914 8927666790 1467771423 780 d0 /
      data wg  (  3) / 0.3818300505 0511894495 0369775488 975 d0 /
      data wg  (  4) / 0.4179591836 7346938775 5102040816 327 d0 /
c
      data xgk (  1) / 0.9914553711 2081263920 6854697526 329 d0 /
      data xgk (  2) / 0.9491079123 4275852452 6189684047 851 d0 /
      data xgk (  3) / 0.8648644233 5976907278 9712788640 926 d0 /
      data xgk (  4) / 0.7415311855 9939443986 3864773280 788 d0 /
      data xgk (  5) / 0.5860872354 6769113029 4144838258 730 d0 /
      data xgk (  6) / 0.4058451513 7739716690 6606412076 961 d0 /
      data xgk (  7) / 0.2077849550 0789846760 0689403773 245 d0 /
      data xgk (  8) / 0.0000000000 0000000000 0000000000 000 d0 /
c
      data wgk (  1) / 0.0229353220 1052922496 3732008058 970 d0 /
      data wgk (  2) / 0.0630920926 2997855329 0700663189 204 d0 /
      data wgk (  3) / 0.1047900103 2225018383 9876322541 518 d0 /
      data wgk (  4) / 0.1406532597 1552591874 5189590510 238 d0 /
      data wgk (  5) / 0.1690047266 3926790282 6583426598 550 d0 /
      data wgk (  6) / 0.1903505780 6478540991 3256402421 014 d0 /
      data wgk (  7) / 0.2044329400 7529889241 4161999234 649 d0 /
      data wgk (  8) / 0.2094821410 8472782801 2999174891 714 d0 /
c
c
c           list of major variables
c           -----------------------
c
c           centr  - mid point of the interval
c           hlgth  - half-length of the interval
c           absc   - abscissa
c           fval*  - function value
c           resg   - result of the 7-point gauss formula
c           resk   - result of the 15-point kronrod formula
c           reskh  - approximation to the mean value of f over (a,b),
c                    i.e. to i/(b-a)
c
c           machine dependent constants
c           ---------------------------
c
c           epmach is the largest relative spacing.
c           uflow is the smallest positive magnitude.
c
c***first executable statement  dqk15
      epmach = d1mach(4)
      uflow = d1mach(1)
c
      centr = 0.5d+00*(a+b)
      hlgth = 0.5d+00*(b-a)
      dhlgth = dabs(hlgth)
c
c           compute the 15-point kronrod approximation to
c           the integral, and estimate the absolute error.
c
      fc = f(centr)
      resg = fc*wg(4)
      resk = fc*wgk(8)
      resabs = dabs(resk)
      do 10 j=1,3
        jtw = j*2
        absc = hlgth*xgk(jtw)
        fval1 = f(centr-absc)
        fval2 = f(centr+absc)
        fv1(jtw) = fval1
        fv2(jtw) = fval2
        fsum = fval1+fval2
        resg = resg+wg(j)*fsum
        resk = resk+wgk(jtw)*fsum
        resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2))
   10 continue
      do 15 j = 1,4
        jtwm1 = j*2-1
        absc = hlgth*xgk(jtwm1)
        fval1 = f(centr-absc)
        fval2 = f(centr+absc)
        fv1(jtwm1) = fval1
        fv2(jtwm1) = fval2
        fsum = fval1+fval2
        resk = resk+wgk(jtwm1)*fsum
        resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2))
   15 continue
      reskh = resk*0.5d+00
      resasc = wgk(8)*dabs(fc-reskh)
      do 20 j=1,7
        resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh))
   20 continue
      result = resk*hlgth
      resabs = resabs*dhlgth
      resasc = resasc*dhlgth
      abserr = dabs((resk-resg)*hlgth)
      if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00)
     *  abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00)
      if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1
     *  ((epmach*0.5d+02)*resabs,abserr)
      return
      end

max_stars_repo_name: Ennosigaeon/scipy

max_issues_repo_name: Ennosigaeon/scipy

max_forks_repo_name: Ennosigaeon/scipy


text:
      subroutine magi(fractionOfDB,
     &                table,
     &                u_ubvrij,
     &                b_ubvrij,
     &                v_ubvrij,
     &                r_ubvrij,
     &                i_ubvrij,
     &                j_ubvrij)
      use external_types
C     This subroutine calculates ltc,cbv,cvi,cvr,cuv visual absolute
C     and apparent magnitude of the WDs.
      implicit real (a-h,m,o-z)

      integer numberOfStars,iseed,i,numberOfWDs,in
      real lum,teff,xlog,c1,c2,c3,c4,c5,c6,n1,n2,n3,n4,n5
      real fractionOfDB
      real mone

      parameter (numberOfStars=6000000)
      parameter (mone=1.14)

      real :: luminosityOfWD(numberOfStars),
     &        massOfWD(numberOfStars),
     &        metallicityOfWD(numberOfStars),
     &        effTempOfWD(numberOfStars)
      integer :: flagOfWD(numberOfStars)
      real :: coolingTime(numberOfStars)
      integer typeOfWD(numberOfStars)
      integer disk_belonging(numberOfStars)

      real :: u_ubvrij(numberOfStars),
     &        b_ubvrij(numberOfStars),
     &        v_ubvrij(numberOfStars),
     &        r_ubvrij(numberOfStars),
     &        i_ubvrij(numberOfStars),
     &        j_ubvrij(numberOfStars)

      real :: rgac(numberOfStars)
      double precision :: lgac(numberOfStars),
     &                    bgac(numberOfStars)

      TYPE(FileGroupInfo),DIMENSION(11) :: table

      common /enanas/ luminosityOfWD,massOfWD,metallicityOfWD,
     &                effTempOfWD
      common /index/ flagOfWD,numberOfWDs,disk_belonging
      common /paral/ rgac
      common /cool/ coolingTime
      common /indexdb/ typeOfWD
      common /lb/ lgac,bgac

      n1=0
      n2=0
      n3=0
      n4=0
      n5=0

C     ---  Interpolating Mv, luminosity, colors and other variables
C          from coolingTime and the mwd  ---
      do i = 1, numberOfWDs
C       ---  ATENTION! choosing only if .lt.1.1!!!  ---
        if (massOfWD(i) .le. 1.4) then
C         ---  IF CO core ---
          if (massOfWD(i) .lt. mone) then
C           --- Atention We put the "old" ones to 0.6Msol ---
C           --- Distribucion DA/DB ---
            call dbd_fid(iseed,fractionOfDB,in)
C           --- End of distribution ---
C           ---  IF DA ---
            if(in.eq.0) then
              typeOfWD(i)=0
              n1=n1+1
              call interlumda(coolingTime(i),massOfWD(i),
     &             metallicityOfWD(i),lum,teff,xlog,c1,c2,c3,c4,c5,c6,
     &             table)
C           ---  ELSE DB  ---
            else
              n3=n3+1
              typeOfWD(i)=1
              call interlumdb(coolingTime(i),massOfWD(i),
     &             metallicityOfWD(i),lum,c1,c2,c3,c4,c5,c6,teff,xlog,
     &             table)
              if(teff.lt.6000) n5=n5+1
            end if
C           ---  END IF DB/NON-DB
C         ---  ELSE ONe ---
          else
            n2=n2+1
            typeOfWD(i)=2
            call interlumone(coolingTime(i),massOfWD(i),lum,c1,c2,c3,c4,
     &           c5,c6,teff,xlog)
          end if
C         ---  END IF CO/ONe ---
          if (teff < 6000) then
            n4 = n4 + 1
          end if

          luminosityOfWD(i) = -lum
          effTempOfWD(i) = teff

          u_ubvrij(i) = c1
          b_ubvrij(i) = c2
          v_ubvrij(i) = c3
          r_ubvrij(i) = c4
          i_ubvrij(i) = c5
          j_ubvrij(i) = c6

C       ---  ELSE mass >= 1.4  --- EXPLOTA, exceeding Chandrasekar limit
        else
          typeOfWD(i) = 5
        end if
      end do

      write(*,*) "DA CO ",n1
      write(*,*) "DA ONe ",n2
      write(*,*) "DB",n3
      write(*,*) "<6000 DA", n4
      write(*,*) "<6000 DB", n5

      return
      end

max_stars_repo_name: wolvespack/alcor

max_issues_repo_name: wolvespack/alcor

max_forks_repo_name: wolvespack/alcor


text:
! MIT License
!
! Copyright (c) 2010-present David A. Kopriva and other contributors: AUTHORS.md
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!
! FTObjectLibrary contains code that, to the best of our knowledge, has been released as
! public domain software:
! * `b3hs_hash_key_jenkins`: originally by Rich Townsend,
! https://groups.google.com/forum/#!topic/comp.lang.fortran/RWoHZFt39ng, 2005
!
! --- End License

!
!////////////////////////////////////////////////////////////////////////
!
!      SparseMatrixTests.f90
!      Created: July 29, 2013 2:17 PM
!      By: David Kopriva
!
!////////////////////////////////////////////////////////////////////////
!
      SUBROUTINE SparseMatrixTests
         USE FTValueClass
         USE FTAssertions
         USE FTSparseMatrixClass
         IMPLICIT NONE
!
!        --------------------------------
!        We will store values in a matrix
!        --------------------------------
!
         TYPE (FTValue) , POINTER :: v
         CLASS(FTValue) , POINTER :: vTest
         CLASS(FTObject), POINTER :: obj

         TYPE(FTSparseMatrix), POINTER :: SparseMatrix

         INTEGER :: i, j, N
         INTEGER :: h1, h2
         LOGICAL :: t
!
!        -------------------------------------
!        initialize a table with four elements
!        -------------------------------------
!
         N = 4
         ALLOCATE(SparseMatrix)
         CALL SparseMatrix % initWithSize(4)
         CALL FTAssertEqual(N,SparseMatrix % SparseMatrixSize(),"Table size size")

         obj => SparseMatrix % objectForKeys(2,3)
         CALL FTAssertEqual(.FALSE.,ASSOCIATED(obj),"Empty table test")
!
!        ---------------------------------------------------------
!        Add an object to the table, retrieve it and then destroy
!        the table
!        ---------------------------------------------------------
!
         ALLOCATE(v)
         CALL v % initWithValue(42)

         obj => v
         CALL SparseMatrix % addObjectForKeys(obj,2,3)
         CALL FTAssertEqual(2,v % refcount(),"Add object to table reference count")

         t = SparseMatrix % containsKeys(2,3)
         CALL FTAssertEqual(.TRUE.,t,"contains value for key")

         obj   => SparseMatrix % objectforKeys(2,3)
         vTest => valueFromObject(obj)
         CALL FTAssertEqual(42,vTest % integerValue(),"Table entry retrieval")

         CALL releaseFTSparseMatrix(SparseMatrix)

         CALL FTAssertEqual(1, v % refCount(),"Table release object refCount")
         CALL releaseFTValue(v)
         CALL FTAssertEqual(.FALSE.,ASSOCIATED(v),"Released object should have been deallocated")
!
!        ---------------------------------------------
!        Now create a hash table with a lot of entries
!        ---------------------------------------------
!
         ALLOCATE(SparseMatrix)
         CALL SparseMatrix % initWithSize(4)
         DO j = 1,N
            DO i = 1,N
               h1 = Hash1([i,j])
               h2 = Hash2([i,j])

               ALLOCATE(v)
               CALL v % initWithValue(i+j)
               obj => v
               CALL SparseMatrix % addObjectForKeys(obj,h1,h2)
               CALL releaseFTValue(v)

            END DO
         END DO
!
!        ------------
!        Get them out
!        ------------
!
         DO j = 1,N
            DO i = 1,N
               h1    =  Hash1([i,j])
               h2    =  Hash2([i,j])
               obj   => SparseMatrix % objectForKeys(h1,h2)
               vTest => valueFromObject(obj)

               CALL FTAssertEqual(1  , vTest % refCount(),"Table add refCount")
               CALL FTAssertEqual(i+j, vTest % integerValue(), "Table entry value")

            END DO
         END DO
         CALL vTest % retain()
         CALL FTAssertEqual(2  , vTest % refCount(),"Retain refCount")

         CALL releaseFTSparseMatrix(SparseMatrix)

         CALL FTAssertEqual(1, vTest % refCount(),"Table release object refCount")
         CALL releaseFTValue(v)
         IF ( ASSOCIATED(v) )     THEN
            CALL FTAssert(.FALSE.,"Release object count")
         END IF

      END SUBROUTINE SparseMatrixTests

max_stars_repo_name: HigherOrderMethods/FTObjectLibrary

max_issues_repo_name: HigherOrderMethods/FTObjectLibrary

max_forks_repo_name: HigherOrderMethods/FTObjectLibrary


text:
      SUBROUTINE FRMLTA (IFILE,Z,Y,ZM)
C
C     LOWER TRIANGULAR TRANSPOSE WITH OFF-DIAGONAL SWITCH
C     SINGLE PRECISION VERSION
C
C     LAST REVISED  11/91, BY G.CHAN/UNISYS
C     ADDITIONAL OF A NEW METHODS WHICH IS MORE EFFICIENT, AND IS
C     ALREADY GOOD FOR VECTORIZATION
C
      REAL             Z(1)    ,Y(1)    ,ZM(1)
      DIMENSION        IFILE(7),NAM(2)
      COMMON  /UNPAKX/ ITYP    ,IP      ,NP      ,INCR
      COMMON  /FEERXX/ DM18(18),NZM
      COMMON  /ZZZZZZ/ IZ(1)
      COMMON  /SYSTEM/ IBUF    ,NOUT
      EQUIVALENCE      (DP,IDP)
      DATA     NAM   / 4HFRML  ,4HTA  /
C
      N    = IFILE(2)
      IFL  = IFILE(1)
      IF (IFILE(7) .LT. 0) IFL = -IFILE(7)
      CALL REWIND (IFL)
      IF (IFILE(7) .LT. 0) GO TO 30
      CALL SKPREC (IFL,1)
      ITYP = IFILE(5)
C
C     NASTRAN ORIGINAL METHOD
C
      INCR = 1
      DO 20 I = 1,N
      Y(I) = 0.0
      IP   = 0
      CALL UNPACK (*30,IFL,ZM(1))
      IF (IP .EQ. I) ZM(1) = -ZM(1)
      SUM  = 0.0
      II   = 0
      DO 10 J = IP,NP
      II   = II + 1
   10 SUM  = SUM - ZM(II)*Z(J)
      Y(I) = SUM
   20 CONTINUE
      GO TO 150
C
C     NEW METHOD
C
C     UNLIKE FRMLT, IFL WAS UNPACKED BACKWARD FIRST, THEN FORWARD BY
C     UNPSCR/FEER3. SO WE SKIP BACKWARD PASS BEFORE READING DATA
C
   30 NREC = IFILE(4)/10
      CALL SKPREC (IFL,NREC+1)
      NREC = 0
      LL2  = 0
      NTMS = 1
      DO 70 I = 1,N
      IF (NTMS .LT. LL2) GO TO 50
      NREC = NREC + 1
      CALL READ (*100,*40,IFL,ZM,NZM,1,LL)
      CALL MESAGE (-8,0,NAM)
   40 LL2  = LL
      NTMS = 1
   50 DP   = ZM(NTMS)
      II   = IDP
      IF (II .NE. I) GO TO 120
      DP   = ZM(NTMS+1)
      JJ   = IDP
      ZM(NTMS+2) = -ZM(NTMS+2)
      SUM  = 0.0
      LL   = NTMS + 1
      DO 60 J = II,JJ
      LL   = LL + 1
   60 SUM  = SUM - ZM(LL)*Z(J)
      Y(I) = SUM
   70 NTMS = NTMS + JJ - II + 3
      GO TO 150
C
  100 J = IFILE(4)/10
      WRITE  (NOUT,110) NREC,I,N,J
  110 FORMAT ('0*** TRY TO READ RECORD',I5,'.  I,N,IFILE(4) =',2I7,I5)
      CALL MESAGE (-2,ILF,NAM)
  120 WRITE  (NOUT,130) II,I
  130 FORMAT ('0*** II AND I MISMATCH =',2I8)
      CALL MESAGE (-37,0,NAM)
C
  150 RETURN
      END


text:
C
C	$Id: gwilod.f,v 1.4 2008-07-27 00:21:07 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE GWILOD(LIST,BITS,COUNT,GKSERR)
C
C  Load the bit string contained in the low order part of each word.
C
C  INPUT
C    LIST   -  A list of words which have bit strings right justified.
C    BITS   -  Number of bits per word to move.
C    COUNT  -  The repitition count (number of words in the list).
C
C  OUTPUT
C    GKSERR -  The error status.
C
C  All data is type INTEGER unless otherwise indicated.
C
      IMPLICIT INTEGER (A-Z)
C
      DIMENSION LIST(*)
C
C  Common for metafile buffer.
C
      include 'gwiio.h'
C
      SAVE
C
      DATA  ALLOK/0/
C
      IF (COUNT.LE.0)  RETURN
      CTEMP = COUNT
      STRT = 1
C
C  Determine the number of packets of size bits left in buffer.
C
 10   CONTINUE
      BLEFT = (WXBITS-WBFPOS)
      REPLFT = BLEFT/BITS
C
C  Compute how may packets to move into the buffer.
C
      IF (CTEMP .LE. REPLFT) THEN
C
C  Room for all.
C
        CMOVE = CTEMP
        CTEMP = 0
      ELSE
C
C  Not enough room for all bits.
C
        CMOVE = REPLFT
        CTEMP = CTEMP - REPLFT
      END IF
C
C  Move the current bit packets.
C
      CALL SBYTES(WOUTBF,LIST(STRT),WBFPOS,BITS,0,CMOVE)
      WBFPOS = WBFPOS + (BITS * CMOVE)
C
C  Check if more bit packets remaining.
C
      IF (CTEMP.NE.0) THEN
C
C  More remaining, flush the current buffer.
C
        CALL GWIFLB(GKSERR)
        IF (GKSERR .NE. ALLOK) RETURN
        STRT = STRT + CMOVE
        GO TO 10
      END IF
C
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
-makelib xcelium_lib/xilinx_vip -sv \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi4stream_vip_axi4streampc.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi_vip_axi4pc.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/xil_common_vip_pkg.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi4stream_vip_pkg.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi_vip_pkg.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi4stream_vip_if.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/axi_vip_if.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/clk_vip_if.sv" \
  "C:/vivado/Vivado/2018.3/data/xilinx_vip/hdl/rst_vip_if.sv" \
-endlib
-makelib xcelium_lib/axi_infrastructure_v1_1_0 \
  "../../../../hello_word.srcs/sources_1/bd/hello_word/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib xcelium_lib/axi_vip_v1_1_4 -sv \
  "../../../../hello_word.srcs/sources_1/bd/hello_word/ipshared/98af/hdl/axi_vip_v1_1_vl_rfs.sv" \
-endlib
-makelib xcelium_lib/processing_system7_vip_v1_0_6 -sv \
  "../../../../hello_word.srcs/sources_1/bd/hello_word/ipshared/70cf/hdl/processing_system7_vip_v1_0_vl_rfs.sv" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  "../../../bd/hello_word/ip/hello_word_processing_system7_0_0/sim/hello_word_processing_system7_0_0.v" \
  "../../../bd/hello_word/sim/hello_word.v" \
-endlib
-makelib xcelium_lib/xil_defaultlib \
  glbl.v
-endlib


text:
program global_work
  call configure_arrays     ! Perform the appropriate allocations
  call compute              ! Use the arrays in computations
end program global_work
module work_arrays          ! An example set of work arrays
  integer n
  real, allocatable, save :: a(:), b(:,:), c(:,:,:)
end module work_arrays
subroutine configure_arrays ! Process to set up work arrays
  use work_arrays
  read (*,*) N
  allocate (a(n), b(n,n), c(n,n,2*n))
end subroutine configure_arrays
subroutine compute
  use work_arrays
  ! ... Computations involving arrays A, B, and C
end subroutine compute

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
Jump to Timeline #Navigation Navigation

  Ben A Madson born


text:
        SUBROUTINE PUTACT (ITEM)
C       PRINT THE ACTION OF AN ITEM
C       GF 11.07.1980
C
        INCLUDE 'PARS.f'
        INCLUDE 'ITES.f'
        INCLUDE 'PROS.f'
        INTEGER*2
     =  ACT     ! THE ACTION OF THE PARSER
     = ,ITEM    ! THE ITEM THE ACTION OF WHICH IS TO BE PRINTED
     = ,PROD    ! THE PRODUCTION FOR 'ACT = REDUCE'
     = ,MNEMO(4) ! THE MNEMONICS FOR THE ACTIONS
        DATA MNEMO/'->', '=:', '=.', '??'/
C
        ACT = ITEACT(ITEM)
        CALL ZZWC (MNEMO(ACT),1,2,3)
        PROD = ITESUC (ITEM)
        CALL ZZWI (PROD,0)
        IF (ACT .NE. REDUCE) GOTO 2
          CALL ZZWC (': ',1,2,0)
          CALL OUTMAP (PROD,0)
2       CONTINUE
        RETURN ! PUTACT
        END

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE ADDDSP (COORDS, DSP)
      DIMENSION COORDS (NUMNP,*), DSP(NUMNP,*)
CC
      include 'nu_numg.blk'
      IF (NDIM .EQ. 2) THEN
         DO 10 J=1,NUMNP
            DSP (J, 1) = COORDS (J, 1) + DSP (J, 1)
            DSP (J, 2) = COORDS (J, 2) + DSP (J, 2)
   10    CONTINUE
      ELSE
         DO 20 J=1,NUMNP
            DSP (J, 1) = COORDS (J, 1) + DSP (J, 1)
            DSP (J, 2) = COORDS (J, 2) + DSP (J, 2)
            DSP (J, 3) = COORDS (J, 3) + DSP (J, 3)
   20    CONTINUE
      END IF
CC
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
      SUBROUTINE KSLIST(IPAIR)
*
*
*       KS perturber selection.
*       -----------------------
*
      INCLUDE 'common6.h'
*
*
*       Set component & c.m. index and form semi-major axis & eccentricity.
      I1 = 2*IPAIR - 1
      I = N + IPAIR
      SEMI = -0.5D0*BODY(I)/H(IPAIR)
      EB = -0.5*BODY(I1)*BODY(I1+1)/SEMI
      ECC2 = (1.0 - R(IPAIR)/SEMI)**2 + TDOT2(IPAIR)**2/(BODY(I)*SEMI)
      ECC = SQRT(ECC2)
*     IT = 0
*
*       Use semi-major axis and/or RMIN for perturber selection.
      IF (EB.LT.EBH) THEN
          RAP = SEMI*(1.0 + ECC)
      ELSE
*       Include a tapered criterion depending on energy for soft binaries.
          IF (EB.LT.0.0) THEN
              ZFAC = 1.0 + ABS(EB - EBH)/ABS(EBH)
          ELSE
              ZFAC = 1.0
          END IF
*       Adopt actual apocentre for perturber selection if R > SEMI.
          IF (SEMI.GT.0.0D0.AND.R(IPAIR).GT.SEMI) THEN
              RAP = SEMI*(1.0 + ECC)
          ELSE
*       Accept separation limit for hyperbolic case.
              RAP = MAX(ZFAC*SEMI,R(IPAIR))
          END IF
*       Ensure extra perturbers at new regularizations (R may be small).
          IF (IPHASE.GT.0.AND.SEMI.GT.0.0) THEN
              RAP = SEMI*(1.0 + ECC)
          END IF
      END IF
*
*       Restrict perturber selection for massive eccentric binary.
      IF (BODY(I).GT.100.0*BODYM.AND.
     &   (R(IPAIR).LT.0.1*SEMI.OR.SEMI.LT.0)) THEN
          RAP = MIN(RMIN,0.1*ABS(SEMI))
      END IF
*
      RCRIT2 = 2.0*RAP**2/BODY(I)
      RCRIT3 = RCRIT2*RAP/GMIN
*       Base fast search on maximum single perturber mass (BODY1).
      RCRIT2 = RCRIT2*BODY1*CMSEP2
      RCRIT6 = RCRIT3**2
*
*       Select new perturbers from the neighbour list.
*   6 NNB1 = 1
      NNB1 = 1
      NNB2 = LIST(1,I) + 1
      DO 10 L = 2,NNB2
          J = LIST(L,I)
          W1 = X(1,J) - X(1,I)
          W2 = X(2,J) - X(2,I)
          W3 = X(3,J) - X(3,I)
          RSEP2 = W1*W1 + W2*W2 + W3*W3
*       Include any merged c.m. or chain c.m. bodies in the fast test.
          IF (RSEP2.LT.RCRIT2.OR.NAME(J).LE.0) THEN
*       Estimate unperturbed distance from tidal limit approximation.
              IF (RSEP2**3.LT.BODY(J)**2*RCRIT6) THEN
                  NNB1 = NNB1 + 1
                  LIST(NNB1,I1) = J
              ELSE IF (J.GT.N) THEN
*       Employ a more generous criterion for possible wide binary.
                  RJ = MIN(10.0*ABS(SEMI),-BODY(J)/H(J-N))
                  IF (RSEP2.LT.CMSEP2*RJ**2) THEN
                      NNB1 = NNB1 + 1
                      LIST(NNB1,I1) = J
                  END IF
              END IF
          END IF
   10 CONTINUE
*
*       Ensure at least one perturber first time (max 5 tries except CHAOS).
*     IF (NNB1.EQ.1.AND.IPHASE.EQ.0.AND.NNB2.GT.1.AND.TIME.GT.0.0) THEN
*         RCRIT2 = 2.0*RCRIT2  ! try to avoid IPHASE for parallel version!
*         RCRIT6 = 2.0*RCRIT6
*         IT = IT + 1
*       Skip repeat for small size (next KSLIST requires many periods).
*         IF ((SEMI*SU.GT.10.0.AND.IT.LE.5).OR.KSTAR(I).EQ.-1) GO TO 6
*     END IF
*
*       Check case of no perturbers (dual purpose).
      IF (NNB1.EQ.1) THEN
*       Add distant perturber for hyperbolic orbit.
          IF (SEMI.LT.0.0) THEN
              NNB1 = 2
              LIST(2,I1) = LIST(2,I)
              GO TO 20
          END IF
*
*       Restrict look-up time to one period for active PN binary (< 1000*RZ).
          IF (KZ(11).GT.0) THEN
              RP = SEMI*(1.0 - ECC)
              IF (RP.LT.1000.0*RZ) THEN
                  IGR = 1
              ELSE
                  IGR = 0
              END IF
*        Note usual transition from perturbed to unperturbed state in KSINT.
              IF (IGR.GT.0) THEN
                  STEP(I1) = TWOPI*SEMI*SQRT(SEMI/BODY(I))
                  DT = 2.0*STEP(I1)
                  CALL STEPK(DT,DTN)
                  STEP(I1) = DTN
                  GO TO 20
              END IF
          END IF
*
          IF (KZ(27).LE.0) THEN
*       Specify one unperturbed period at apocentre (NB! check STEP(I)).
              STEP(I1) = TWOPI*SEMI*SQRT(SEMI/BODY(I))
              STEP(I1) = MIN(STEP(I1),STEP(I))
*       Include quantization (interval shortened a little so compensate).
              DT = 2.0*STEP(I1)
              CALL STEPK(DT,DTN)
              STEP(I1) = DTN
          ELSE
*       Maintain perturbed motion during Chaos event (not ROCHE/SPIRAL).
              IF (KSTAR(I).EQ.-1) THEN
                  IF (LIST(1,I1).GT.0) THEN
                      NNB1 = 2
                      LIST(2,I1) = N
                  END IF
              ELSE
*                 STEP(I1) = TWOPI*SEMI*SQRT(SEMI/BODY(I))
*                 STEP(I1) = MIN(STEP(I1),STEP(I))
*       Avoid possible small period for standard binary (hence use c.m.).
                  STEP(I1) = STEP(I)
                  DT = 2.0*STEP(I1)
                  CALL STEPK(DT,DTN)
                  STEP(I1) = DTN
*       Note reduction in BRAKE4 in case of active PN phase.
              END IF
          END IF
      END IF
*
*       Save perturber membership.
   20 LIST(1,I1) = NNB1 - 1
*
      RETURN
*
      END


text:
      SUBROUTINE slPLTE (DATE, ELONG, PHI, JFORM, EPOCH,
     :                       ORBINC, ANODE, PERIH, AORQ, E,
     :                       AORL, DM, RA, DEC, R, JSTAT)
*+
*     - - - - - - -
*      P L T E
*     - - - - - - -
*
*  Topocentric apparent RA,Dec of a Solar-System object whose
*  heliocentric orbital elements are known.
*
*  Given:
*     DATE     d     MJD of observation (JD - 2400000.5, Notes 1,5)
*     ELONG    d     observer's east longitude (radians, Note 2)
*     PHI      d     observer's geodetic latitude (radians, Note 2)
*     JFORM    i     choice of element set (1-3; Notes 3-6)
*     EPOCH    d     epoch of elements (TT MJD, Note 5)
*     ORBINC   d     inclination (radians)
*     ANODE    d     longitude of the ascending node (radians)
*     PERIH    d     longitude or argument of perihelion (radians)
*     AORQ     d     mean distance or perihelion distance (AU)
*     E        d     eccentricity
*     AORL     d     mean anomaly or longitude (radians, JFORM=1,2 only)
*     DM       d     daily motion (radians, JFORM=1 only )
*
*  Returned:
*     RA,DEC   d     RA, Dec (topocentric apparent, radians)
*     R        d     distance from observer (AU)
*     JSTAT    i     status:  0 = OK
*                            -1 = illegal JFORM
*                            -2 = illegal E
*                            -3 = illegal AORQ
*                            -4 = illegal DM
*                            -5 = numerical error
*
*  Called: slELUE, slPLTU
*
*  Notes:
*
*  1  DATE is the instant for which the prediction is required.  It is
*     in the TT timescale (formerly Ephemeris Time, ET) and is a
*     Modified Julian Date (JD-2400000.5).
*
*  2  The longitude and latitude allow correction for geocentric
*     parallax.  This is usually a small effect, but can become
*     important for near-Earth asteroids.  Geocentric positions can be
*     generated by appropriate use of routines slEVP (or slEPV) and
*     slPLNE.
*
*  3  The elements are with respect to the J2000 ecliptic and equinox.
*
*  4  A choice of three different element-set options is available:
*
*     Option JFORM = 1, suitable for the major planets:
*
*       EPOCH  = epoch of elements (TT MJD)
*       ORBINC = inclination i (radians)
*       ANODE  = longitude of the ascending node, big omega (radians)
*       PERIH  = longitude of perihelion, curly pi (radians)
*       AORQ   = mean distance, a (AU)
*       E      = eccentricity, e (range 0 to <1)
*       AORL   = mean longitude L (radians)
*       DM     = daily motion (radians)
*
*     Option JFORM = 2, suitable for minor planets:
*
*       EPOCH  = epoch of elements (TT MJD)
*       ORBINC = inclination i (radians)
*       ANODE  = longitude of the ascending node, big omega (radians)
*       PERIH  = argument of perihelion, little omega (radians)
*       AORQ   = mean distance, a (AU)
*       E      = eccentricity, e (range 0 to <1)
*       AORL   = mean anomaly M (radians)
*
*     Option JFORM = 3, suitable for comets:
*
*       EPOCH  = epoch of elements and perihelion (TT MJD)
*       ORBINC = inclination i (radians)
*       ANODE  = longitude of the ascending node, big omega (radians)
*       PERIH  = argument of perihelion, little omega (radians)
*       AORQ   = perihelion distance, q (AU)
*       E      = eccentricity, e (range 0 to 10)
*
*     Unused arguments (DM for JFORM=2, AORL and DM for JFORM=3) are not
*     accessed.
*
*  5  Each of the three element sets defines an unperturbed heliocentric
*     orbit.  For a given epoch of observation, the position of the body
*     in its orbit can be predicted from these elements, which are
*     called "osculating elements", using standard two-body analytical
*     solutions.  However, due to planetary perturbations, a given set
*     of osculating elements remains usable for only as long as the
*     unperturbed orbit that it describes is an adequate approximation
*     to reality.  Attached to such a set of elements is a date called
*     the "osculating epoch", at which the elements are, momentarily,
*     a perfect representation of the instantaneous position and
*     velocity of the body.
*
*     Therefore, for any given problem there are up to three different
*     epochs in play, and it is vital to distinguish clearly between
*     them:
*
*     . The epoch of observation:  the moment in time for which the
*       position of the body is to be predicted.
*
*     . The epoch defining the position of the body:  the moment in time
*       at which, in the absence of purturbations, the specified
*       position (mean longitude, mean anomaly, or perihelion) is
*       reached.
*
*     . The osculating epoch:  the moment in time at which the given
*       elements are correct.
*
*     For the major-planet and minor-planet cases it is usual to make
*     the epoch that defines the position of the body the same as the
*     epoch of osculation.  Thus, only two different epochs are
*     involved:  the epoch of the elements and the epoch of observation.
*
*     For comets, the epoch of perihelion fixes the position in the
*     orbit and in general a different epoch of osculation will be
*     chosen.  Thus, all three types of epoch are involved.
*
*     For the present routine:
*
*     . The epoch of observation is the argument DATE.
*
*     . The epoch defining the position of the body is the argument
*       EPOCH.
*
*     . The osculating epoch is not used and is assumed to be close
*       enough to the epoch of observation to deliver adequate accuracy.
*       If not, a preliminary call to slPRTL may be used to update
*       the element-set (and its associated osculating epoch) by
*       applying planetary perturbations.
*
*  6  Two important sources for orbital elements are Horizons, operated
*     by the Jet Propulsion Laboratory, Pasadena, and the Minor Planet
*     Center, operated by the Center for Astrophysics, Harvard.
*
*     The JPL Horizons elements (heliocentric, J2000 ecliptic and
*     equinox) correspond to SLALIB arguments as follows.
*
*       Major planets:
*
*         JFORM  = 1
*         EPOCH  = JDCT-2400000.5D0
*         ORBINC = IN (in radians)
*         ANODE  = OM (in radians)
*         PERIH  = OM+W (in radians)
*         AORQ   = A
*         E      = EC
*         AORL   = MA+OM+W (in radians)
*         DM     = N (in radians)
*
*         Epoch of osculation = JDCT-2400000.5D0
*
*       Minor planets:
*
*         JFORM  = 2
*         EPOCH  = JDCT-2400000.5D0
*         ORBINC = IN (in radians)
*         ANODE  = OM (in radians)
*         PERIH  = W (in radians)
*         AORQ   = A
*         E      = EC
*         AORL   = MA (in radians)
*
*         Epoch of osculation = JDCT-2400000.5D0
*
*       Comets:
*
*         JFORM  = 3
*         EPOCH  = Tp-2400000.5D0
*         ORBINC = IN (in radians)
*         ANODE  = OM (in radians)
*         PERIH  = W (in radians)
*         AORQ   = QR
*         E      = EC
*
*         Epoch of osculation = JDCT-2400000.5D0
*
*     The MPC elements correspond to SLALIB arguments as follows.
*
*       Minor planets:
*
*         JFORM  = 2
*         EPOCH  = Epoch-2400000.5D0
*         ORBINC = Incl. (in radians)
*         ANODE  = Node (in radians)
*         PERIH  = Perih. (in radians)
*         AORQ   = a
*         E      = e
*         AORL   = M (in radians)
*
*         Epoch of osculation = Epoch-2400000.5D0
*
*       Comets:
*
*         JFORM  = 3
*         EPOCH  = T-2400000.5D0
*         ORBINC = Incl. (in radians)
*         ANODE  = Node. (in radians)
*         PERIH  = Perih. (in radians)
*         AORQ   = q
*         E      = e
*
*         Epoch of osculation = Epoch-2400000.5D0
*
*  This revision:  19 June 2004
*
*  Copyright (C) 2004 P.T.Wallace.  All rights reserved.
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program 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 General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION DATE,ELONG,PHI
      INTEGER JFORM
      DOUBLE PRECISION EPOCH,ORBINC,ANODE,PERIH,AORQ,E,
     :                 AORL,DM,RA,DEC,R
      INTEGER JSTAT

      DOUBLE PRECISION U(13)

*  Transform conventional elements to universal elements.
      CALL slELUE(DATE,
     :               JFORM,EPOCH,ORBINC,ANODE,PERIH,AORQ,E,AORL,DM,
     :               U,JSTAT)

*  If successful, make the prediction.
      IF (JSTAT.EQ.0) CALL slPLTU(DATE,ELONG,PHI,U,RA,DEC,R,JSTAT)

      END

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
module da_gpseph

   use module_domain, only : domain, x_type, xb_type
   use module_dm, only : wrf_dm_sum_real

   use da_control, only : obs_qc_pointer,max_ob_levels,missing_r, &
      v_interp_p, v_interp_h, check_max_iv_print, pi, radian, interpolate_level, &
      missing, max_error_uv, max_error_t, rootproc,fails_error_max, &
      max_error_p,max_error_q, check_max_iv_unit,check_max_iv, qcstat_conv_unit, &
      max_stheight_diff,missing_data,max_error_bq,max_error_slp, ob_vars, &
      max_error_bt, max_error_buv, gpseph,max_error_thickness, &
      pseudo_var, num_pseudo, kms,kme,kts,kte, trace_use_dull, &
      anal_type_verify,fails_error_max,pseudo_err,pseudo_x, pseudo_y, stdout, &
      use_gpsephobs, pseudo_z,pseudo_val,max_error_eph, &
      pseudo, jts, jte,its,ite,ids,ide,jds,jde,kds,kde, gpseph_opt, &
      gps_ray_path_step
   use da_define_structures, only : maxmin_type, iv_type, y_type, jo_type, &
      bad_data_type, number_type, bad_data_type, &
      maxmin_type, da_allocate_observations, gpseph_type
   use da_interpolation, only : da_interp_lin_3d,da_interp_lin_3d_adj, &
      da_to_zk, da_splinx, da_splinx_lin, da_splinx_adj
   use da_par_util, only : da_proc_stats_combine, da_patch_to_global
   use da_par_util1, only : da_proc_sum_int
   use da_statistics, only : da_stats_calculate
   use da_tools, only : da_max_error_qc, da_residual, da_convert_zk,da_get_print_lvl
   use da_tools_serial, only : da_free_unit, da_get_unit
   use da_tracing, only : da_trace_entry, da_trace_exit
   use da_reporting, only : da_error
   use da_wrf_interfaces, only : wrf_dm_bcast_real

   ! The "stats_gpseph_type" is ONLY used locally in da_gpseph:

   type residual_gpseph1_type
      real :: eph                        ! GPS linear excess phase
   end type residual_gpseph1_type

   type maxmin_gpseph_stats_type
      type (maxmin_type)         :: eph  ! GPS excess phase
   end type maxmin_gpseph_stats_type

   type stats_gpseph_type
      type (maxmin_gpseph_stats_type)  :: maximum, minimum
      type (residual_gpseph1_type)     :: average, rms_err
   end type stats_gpseph_type

   TYPE ob_in_mean_h    ! for storing calculated pseudo data in the mean altitude of model
      REAL             , pointer :: lat(:)   ! Latitude
      REAL             , pointer :: lon(:)   ! Longitude
      REAL             , pointer :: ref(:)   ! Refractivity
      REAL             , pointer :: azim(:)  ! Azimuth
   END TYPE ob_in_mean_h

   real, dimension(:,:),   allocatable :: global_lat, global_lon, global_terr
   real, dimension(:,:,:), allocatable :: global_h, global_ref, global_xa_ref, global_adj_ref
   real, dimension(:),     allocatable :: global_h_mean
   real, dimension(:,:,:), allocatable :: global_ref_mean_h

   type gpsnode_type
     integer, dimension(:),  allocatable :: i1,i2,i3
     real, dimension(:),     allocatable :: h
     real, dimension(:,:),   allocatable :: w1
     real, dimension(:,:,:), allocatable :: w2
   end type gpsnode_type

   type gpsrays_type
     integer :: nbot
     integer :: ntop
     integer,dimension(:),allocatable :: je2
     integer,dimension(:),allocatable :: ilocal
     type(gpsnode_type),dimension(:),allocatable :: ip123
   end type gpsrays_type

   type(gpsrays_type), dimension(:), allocatable  :: gps_rays

contains

#include "da_ao_stats_gpseph.inc"
#include "da_calculate_grady_gpseph.inc"
#include "da_jo_and_grady_gpseph.inc"
#include "da_oi_stats_gpseph.inc"
#include "da_print_stats_gpseph.inc"
#include "da_transform_xtoy_gpseph.inc"
#include "da_transform_xtoy_gpseph_adj.inc"
#include "da_check_max_iv_gpseph.inc"
#include "da_get_innov_vector_gpseph.inc"
#include "da_obs_ref_to_eph.inc"
#include "da_mdl_ref_to_eph.inc"
#include "da_residual_gpseph.inc"
#include "da_gpseph_create_ob.inc"
#include "da_gpseph_rays.inc"
#include "da_gpseph_init.inc"
#include "da_gpseph_final.inc"

end module da_gpseph

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
	SUBROUTINE GDPDTA  ( iflno, gdatim, gvcord, gfunc, gpoint,
     +			     time, ivcord, ystrt, ystop,
     +   		     rgx, rgy, rlat,
     +			     rlon, npts, x, y, parm, iret )
C************************************************************************
C* GDPDTA								*
C*									*
C* This subroutine gets the data to plot for a profile.			*
C*									*
C* GDPDTA  ( IFLNO, GDATIM, GVCORD, GFUNC, GPOINT, TIME, IVCORD,	*
C*           YSTRT, YSTOP, RGX, RGY, RLAT, RLON, NPTS, X, Y, PARM, 	*
C*	     IRET )							*
C*									*
C* Input parameters:							*
C*      IFLNO           INTEGER         Grid file number                *
C*	GDATIM		CHAR*		User input date/time		*
C*	GVCORD		CHAR*		User input vert coord		*
C*	GFUNC		CHAR*		User input function		*
C*	GPOINT		CHAR*		User input point to plot	*
C*	TIME  (2)	CHAR*		Time to search for levels	*
C*	IVCORD		INTEGER		Vertical coordinate for search	*
C*      YSTRT	 	REAL            Starting vert coord value	*
C*      YSTOP           REAL            Stopping vert coord value       *
C*									*
C* Output parameters:							*
C*	RGX		REAL		X grid coordinate		*
C*	RGY		REAL		Y grid coordinate		*
C*	RLAT		REAL		Latitude			*
C*	RLON		REAL		Longitude			*
C*	NPTS		INTEGER		Number of points		*
C*	X    (NPTS)	REAL		X coordinates			*
C*	Y    (NPTS)	REAL		Y coordinates			*
C*	PARM		CHAR*		Parameter name			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = invalid grid point	*
C*					 -9 = no valid points		*
C*					-10 = no levels at this time	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	11/85						*
C* M. desJardins/GSFC	 9/88	GEMPAK4					*
C* G. Huffman/GSC	 1/89	Treatment of missing sfc vert. coord.	*
C* M. desJardins/GSFC	 4/89	Added call to LV_CCRD			*
C* M. desJardins/GSFC	 5/89	Call to GR_PLOC				*
C* M. desJardins/GSFC	 5/89	Call to GR_INTP				*
C* M. desJardins/GSFC	 7/89	Reversed x and y in call to GR_INTP	*
C* M. desJardins/GSFC	 7/89	Changed error messages for no points	*
C* K. Brill/GSC		11/89   Added subsetting and layer average parms*
C* S. Schotz/GSC	 6/90	Removed respnd flag			*
C* K. Brill/NMC          9/90   Fix so that surface value is used	*
C* R. Tian/SAIC         10/02   Change call to DG_CXGP                  *
C* K. Brill/HPC		12/02	Pass blank IJSKIP into DG_CXGP		*
C* R. Tian/SAIC		 4/04	Changes for new DG FT mngnmnt		*
C* S. Chiswell/Unidata	 3/06	Added DGCONAV call			*
C* S. Chiswell/Unidata	 2/07	restored output navigation		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	gdatim, gvcord, gfunc, gpoint, time (2), parm
	REAL		x (*), y (*)
	LOGICAL		havsfc, lavflg, havgfs
C*
	CHARACTER	dattim (2)*20, glevel*20, pfunc*80, vcparm*4,
     +			gparm*13
	CHARACTER	ggddtt*48, pparm*12, cbuf*8, dparm*6, cfile
	INTEGER		lev (2)
	REAL		grid ( LLMXGD ), rlvl ( LLMXLV )
	REAL		rnvblk ( LLNNAV ), savblk ( LLNNAV )
	LOGICAL		gintp
C*
	INTEGER		level ( 2, LLMXLV )
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
	npts = 0
C
C*	Find plotting location.
C
	CALL DG_CXGP  ( gpoint, ' ', 1, np, rgx, rgy, rlat, rlon, iret )
	IF ( iret .ne. 0 ) RETURN
C
C*	Check to see if the function involves a layer average.
C
	CALL GDPCLA ( gfunc, lavflg, iret )
	IF ( iret .ne. 0 ) RETURN
C
C*	Get levels which might have data.
C*	First translate date/time and vertical coordinate.
C
	dattim ( 1 ) = time ( 1 )
	dattim ( 2 ) = ' '
	CALL DG_GLEV  ( iflno, dattim, ivcord,
     +	                LLMXLV, level, nlev, ier )
	IF  ( nlev .eq. 0 )  THEN
	    iret = -10
	    CALL ER_WMSG  ( 'GDPROF', iret, ' ', ier )
	    RETURN
	END IF
C
C*	Set a new output reference navigation using gpoint as the center
C*	of a 9 point grid. This allows vector components to be interpolated
C*	to the desired grid point first within DG_GRID before the quantity is
C*	calculated, rather than interpolated from surrounding points after
C*	the computation is done. This will improve calculations involving
C*	vector components.
C
	CALL GR_MNAV ( 'CED', 3, 3, rlat - .01, rlon - .01,
     +		rlat + .01, rlon + .01, 0., 0., 0., .true.,
     +		rnvblk, iret )
	IF ( iret .eq. 0 ) THEN
C
C*	    Save the current navigation block, set new reference navigation
C
	    CALL GD_GNAV  ( iflno, savblk, isavsz, iret )
	    CALL DGCONAV ( rnvblk, iret )
	END IF
	IF ( iret .ne. 0 ) THEN
	    iret = 3
	    CALL ER_WMSG  ( 'GDPROF', iret, ' ', ier )
	    gintp = .true.
	  ELSE
	    gintp = .false.
	END IF
C
C*	Float the levels for sorting and look for surface.
C
	havsfc = .false.
	DO i = 1, nlev
	  rlvl ( i ) = FLOAT ( level ( 1, i ) )
	  IF ( level ( 1, i ) .eq. 0 ) havsfc = .true.
	END DO
C*
	CALL LV_SORT ( ivcord, nlev, rlvl, iret )
C
C*      Assign all levels, temporarily.
C
	DO i = 1, nlev
	  y ( i ) = rlvl ( i )
	END DO
C
C*      If it exists, retrieve vert. coord. value at the surface.
C
	ysfc = RMISSD
	IF ( havsfc ) THEN
		        glevel = '0'
C
C*			Set time, and vertical coordinate to use.  If
C*			retrieve is successful, interpolate value;
C*			otherwise set sfc vert. coord. value to missing.
C
			CALL LV_CCRD ( ivcord, vcparm, ier )
			ggddtt = time (1)
C
C*			Convert file number to character.
C
			CALL ST_INCH ( iflno, cfile, ier )
		        dparm = vcparm // '+' // cfile
			CALL DG_GRID ( ggddtt, glevel, vcparm, dparm,
     +				       pfunc, grid, kx, ky,
     +				       dattim, lev, jvcord, pparm, ier )
			IF  ( ier .ge. 0 )  THEN
			    IF ( gintp ) THEN
			        CALL GR_INTP  ( 1, rgx, rgy, 1, kx, ky,
     +						grid, ysfc, ier )
			      ELSE
				ysfc = grid (5)
			    END IF
			  ELSE
			    ysfc = RMISSD
			END IF
C*
	END IF
C
C*	Set surface value.
C
	IF ( havsfc .and. .not. ERMISS ( ysfc ) )
     +          y ( 1 ) = ysfc
	IF ( havsfc .and. ERMISS ( ysfc ) ) THEN
	  ii = 0
	  DO i = 2, nlev
	    ii = ii + 1
            y ( ii ) = rlvl ( i )
	    rlvl ( ii ) = y ( ii )
       	  END DO
	  nlev = nlev - 1
	  havsfc = .false.
	END IF
C
C*	If level 0 (surface) is present, then HAVSFC is T;
C* 	otherwise HAVSFC is F.  If GFUNC can be computed at
C*	level 0, then HAVGFS is T.
C
	havgfs = .false.
	IF ( havsfc .and. ( .not. lavflg ) ) havgfs = .true.
C*
C
C*	Do subset in the vertical--eliminate unneeded levels.
C
	IF ( havsfc ) THEN
	  IF ( ivcord .eq. 1 .and. ystrt .gt. ysfc ) THEN
	    ys1 = ysfc
	  ELSE IF ( ivcord .ne. 1 .and. ystrt .lt. ysfc ) THEN
	    ys1 = ysfc
          ELSE
	    ys1 = ystrt
          END IF
	ELSE
	    ys1 = ystrt
	END IF
	IF ( ys1 .eq. ystrt ) havgfs = .false.
	i = 1
	istrt = 0
	istop = 0
	DO WHILE ( i .lt. nlev .and.
     +           (istrt .eq. 0 .or. istop .eq. 0 ) )
	  i = i + 1
	  IF ( (ys1 .ge. y ( i-1 ) .and. ys1 .lt. y ( i ) )
     +                               .or.
     +         (ys1 .le. y ( i-1 ) .and. ys1 .gt. y ( i ) ) ) THEN
              istrt = i - 1
	      IF ( havgfs ) rlvl ( istrt ) = rlvl ( 1 )
	  END IF
	  IF ( (ystop .gt. y ( i-1 ) .and. ystop .le. y ( i ) )
     +                               .or.
     +         (ystop .lt. y ( i-1 ) .and. ystop .ge. y ( i ) ) )
     +      istop = i
	END DO
	IF ( istrt .eq. 0 ) istrt = 1
	IF ( istop .eq. 0 ) istop = nlev
	IF ( lavflg ) THEN
	  istrt = istrt - 1
	  istop = istop + 1
	  IF ( istrt .lt. 1 ) istrt = 1
	  IF ( istop .gt. nlev ) istop = nlev
	END IF
C
C*	Loop through single levels finding data.
C
	DO  i = istrt, istop
	  IF ( i .eq. 1 .and. lavflg ) THEN
C
C*	     Do nothing.
C
          ELSE IF ( i .eq. 2 .and. lavflg .and. havsfc ) THEN
C
C*	     Do nothing.
C
	  ELSE
C
C*		Encode level and compute function.
C
	    	IF ( lavflg ) THEN
  		  ilvl1 = int ( rlvl ( i ) )
	          ilvl2 = int ( rlvl ( i - 1 ) )
		  CALL ST_INLN ( ilvl1, glevel, lnth, ier )
		  CALL ST_INCH ( ilvl2, cbuf, ier )
		  glevel = glevel(1:lnth) // ':' // cbuf
	 	ELSE
	          intlvl = int ( rlvl ( i ) )
	 	  CALL ST_INCH ( intlvl, glevel, ier )
		END IF
C*
		CALL DG_GRID  ( gdatim, glevel, gvcord, gfunc,
     +				pfunc,  grid,   kx, ky, dattim, lev,
     +				jvcord, gparm, ierdg )
C*
		IF  ( ierdg .eq. 0 )  THEN
		    parm = gparm
C
		    IF ( gintp ) THEN
C
C*		        Check that grid includes point to be found.
C
		        rkx = FLOAT ( kx )
		        rky = FLOAT ( ky )
		        IF  ( ( rgx .gt. rkx ) .or. ( rgy .gt. rky ) )
     +								 THEN
			    iret = -4
			    RETURN
		        END IF
C
C*		        Interpolate to correct point.
C
		        CALL GR_INTP  ( 1, rgx, rgy, 1, kx, ky, grid,
     +					xx, ier )
		      ELSE
			xx = grid(5)
		    END IF
	            IF ( .not. ERMISS ( xx ) ) THEN
                      npts = npts + 1
	              x (npts) = xx
		      IF ( lavflg ) THEN
		        IF ( ivcord .eq. 1 ) THEN
			  y ( npts ) = SQRT ( rlvl ( i ) *
     +			       		    rlvl ( i - 1 ) )
		        ELSE
		          y ( npts ) = .5 * ( rlvl ( i ) +
     +					    rlvl ( i - 1 ) )
		        END IF
		      ELSE
	                y (npts) = rlvl ( i )
		      END IF
		    ELSE IF ( i .eq. istrt .and. havgfs ) THEN
		      havgfs = .false.
	            END IF
		ELSE IF ( i .eq. istrt .and. havgfs ) THEN
                  havgfs = .false.
		END IF
	  END IF
	END DO
C
C*	Reset the navigation block
C
	CALL DGCONAV ( savblk, ier )
C
C*	Check that there are some points.
C
	IF  ( npts .le. 0 )  THEN
	    CALL ER_WMSG  ( 'DG', ierdg , pfunc, ier )
	    iret = -9
	    CALL ER_WMSG  ( 'GDPROF', iret, ' ', ier )
	    RETURN
	END IF
C
C*	Set vertical coordinate for GFUNC value at surface.
C
	IF ( havgfs ) y ( 1 ) = ysfc
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak