text:
Subroutine hbtout(nnew, nt, ntmax)
  Parameter (maxstr=150001, maxr=1)
  Parameter (oneminus=0.99999, oneplus=1.00001)
  Dimension lastkp(maxstr), newkp(maxstr), xnew(3)
  Common /para7/ioscar, nsmbbbar, nsmmeson
  Common /hbt/lblast(maxstr), xlast(4, maxstr), plast(4, maxstr), nlast
  Common /input1/masspr, massta, iseed, iavoid, dt
  Common /aa/r(3, maxstr)
  Common /bb/p(3, maxstr)
  Common /cc/e(maxstr)
  Common /ee/id(maxstr), lb(maxstr)
  Common /lastt/itimeh, bimp
  Common /tdecay/tfdcy(maxstr), tfdpi(maxstr, maxr), tft(maxstr)
  Common /arevt/iaevt, iarun, miss
  Common /snn/efrm, npart1, npart2, epsipz, epsipt, pzproj, pztarg
  Common /hjglbr/nelt, ninthj, nelp, ninp
  Common /ftmax/ftsv(maxstr), ftsvt(maxstr, maxr)
  Common /dpert/dpertt(maxstr, maxr), dpertp(maxstr), dplast(maxstr), dpdcy(maxstr), dpdpi(maxstr, maxr), dpt(maxstr, maxr), dpp1(maxstr, maxr), dppion(maxstr, maxr)
  Common /hparnt/hipr1(100), ihpr2(50), hint1(100), ihnt2(50)
  External iarflv, invflv
  Common /para8/idpert, npertd, idxsec
  Common /phihj/iphirp, phirp
  Save
  Do i = 1, max0(nlast, nnew)
     lastkp(i) = 0
  End Do
  Do i = 1, nnew
     newkp(i) = 0
  End Do
  Do ip = 1, nnew
     Do iplast = 1, nlast
        If (p(1,ip)==plast(1,iplast) .And. p(2,ip)==plast(2,iplast) .And. p(3,ip)==plast(3,iplast) .And. e(ip)==plast(4,iplast) .And. lb(ip)==lblast(iplast) .And. dpertp(ip)==dplast(iplast) .And. lastkp(iplast)==0) Then
           deltat = nt*dt - xlast(4, iplast)
           ene = sqrt(plast(1,iplast)**2+plast(2,iplast)**2+plast(3,iplast)**2+plast(4,iplast)**2)
           Do ii = 1, 3
              xnew(ii) = xlast(ii, iplast) + plast(ii, iplast)/ene*deltat
           End Do
           dr = sqrt((r(1,ip)-xnew(1))**2+(r(2,ip)-xnew(2))**2+(r(3,ip)-xnew(3))**2)
           If (dr<=0.01) Then
              lastkp(iplast) = 1
              newkp(ip) = 1
              If (nt==ntmax .And. ftsv(ip)>((ntmax-1)*dt)) xlast(4, iplast) = ftsv(ip)
              Goto 100
           End If
        End If
     End Do
100 End Do
  Do ip = 1, nnew
     If (newkp(ip)==0) Then
        Do iplast = 1, nnew
           If (lastkp(iplast)==0) Then
              xlast(1, iplast) = r(1, ip)
              xlast(2, iplast) = r(2, ip)
              xlast(3, iplast) = r(3, ip)
              xlast(4, iplast) = nt*dt
              If (nt==ntmax) Then
                 If (tfdcy(ip)>(ntmax*dt+0.001)) Then
                    xlast(4, iplast) = tfdcy(ip)
                 Else If (ftsv(ip)>((ntmax-1)*dt)) Then
                    xlast(4, iplast) = ftsv(ip)
                 End If
              End If
              plast(1, iplast) = p(1, ip)
              plast(2, iplast) = p(2, ip)
              plast(3, iplast) = p(3, ip)
              plast(4, iplast) = e(ip)
              lblast(iplast) = lb(ip)
              lastkp(iplast) = 1
              dplast(iplast) = dpertp(ip)
              Goto 150
           End If
        End Do
     End If
150 End Do
  If (nnew<nlast) Then
     Do iplast = 1, nlast
        If (lastkp(iplast)==0) Then
           Do ip2 = iplast + 1, nlast
              If (lastkp(ip2)==1) Then
                 xlast(1, iplast) = xlast(1, ip2)
                 xlast(2, iplast) = xlast(2, ip2)
                 xlast(3, iplast) = xlast(3, ip2)
                 xlast(4, iplast) = xlast(4, ip2)
                 plast(1, iplast) = plast(1, ip2)
                 plast(2, iplast) = plast(2, ip2)
                 plast(3, iplast) = plast(3, ip2)
                 plast(4, iplast) = plast(4, ip2)
                 lblast(iplast) = lblast(ip2)
                 lastkp(iplast) = 1
                 dplast(iplast) = dplast(ip2)
                 Goto 170
              End If
           End Do
        End If
170  End Do
  End If
  nlast = nnew
  If (nt==ntmax) Then
     ndpert = 0
     Do ip = 1, nlast
        If (dplast(ip)>oneminus .And. dplast(ip)<oneplus) Then
        Else
           ndpert = ndpert + 1
        End If
     End Do
     Write (16, 191) iaevt, iarun, nlast - ndpert, bimp, npart1, npart2, nelp, ninp, nelt, ninthj, phirp
     If (idpert==1 .Or. idpert==2) Write (90, 190) iaevt, iarun, ndpert, bimp, npart1, npart2, nelp, ninp, nelt, ninthj
     Do ip = 1, nlast
        If (abs(plast(1,ip))<=epsipt .And. abs(plast(2,ip))<=epsipt .And. (plast(3,ip)>amax1(0.,pzproj-epsipz) .Or. plast(3,ip)<(-pztarg+epsipz)) .And. (lblast(ip)==1 .Or. lblast(ip)==2)) Then
           If (dplast(ip)>oneminus .And. dplast(ip)<oneplus) Then
              Write (16, 200) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), plast(4, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip)
           Else
              If (idpert==1 .Or. idpert==2) Then
                 Write (90, 250) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip)
              Else
                 Write (99, *) 'Unexpected perturbative particles'
              End If
           End If
        Else If (amax1(abs(xlast(1,ip)),abs(xlast(2,ip)),abs(xlast(3,ip)),abs(xlast(4,ip)))<9999) Then
           If (dplast(ip)>oneminus .And. dplast(ip)<oneplus) Then
              Write (16, 200) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), plast(4, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip)
           Else
              If (idpert==1 .Or. idpert==2) Then
                 Write (90, 250) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip), dplast(ip)
              Else
                 Write (99, *) 'Unexpected perturbative particles'
              End If
           End If
        Else
           If (dplast(ip)>oneminus .And. dplast(ip)<oneplus) Then
              Write (16, 201) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), plast(4, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip)
           Else
              If (idpert==1 .Or. idpert==2) Then
                 Write (90, 251) invflv(lblast(ip)), plast(1, ip), plast(2, ip), plast(3, ip), xlast(1, ip), xlast(2, ip), xlast(3, ip), xlast(4, ip), dplast(ip)
              Else
                 Write (99, *) 'Unexpected perturbative particles'
              End If
           End If
        End If
     End Do
     If (ioscar==1) Call hoscar
  End If
  Return
190 Format (3(I7), F10.4, 5X, 6(I4))
191 Format (3(I7), F10.4, 5X, 6(I4), 5X, F7.4)
200 Format (I6, 2(1X,F8.3), 1X, F11.4, 1X, F6.3, 4(1X,F8.2))
201 Format (I6, 2(1X,F8.3), 1X, F11.4, 1X, F6.3, 4(1X,E8.2))
250 Format (I5, 2(1X,F8.3), 1X, F10.3, 2(1X,F7.1), 1X, F8.2, 1X, F7.2, 1X, E10.4)
251 Format (I5, 2(1X,F8.3), 1X, F10.3, 4(1X,E8.2), 1X, E10.4)
End Subroutine hbtout

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
#:include 'mpifx.fypp'
#:set TYPES = ALL_TYPES
#:set RANKS = range(1, MAX_RANK + 1)

#! ************************************************************************
#! *** mpifx_allgatherv
#! ************************************************************************

#:def mpifx_allgatherv_dr0_template(SUFFIX, TYPE, RANK, MPI_TYPE)

  !> Gathers results of variable length on all processes (type ${SUFFIX}$).
  !!
  !! \param mycomm  MPI communicator.
  !! \param send  Quantity to be sent for gathering.
  !! \param recv  Received data
  !! \param recvcounts Counts of received data from each process
  !! \param displs Entry i specifies where to place data from process rank i-1
  !!               (default: computed from recvcounts assuming order with rank)
  !! \param error  Error code on exit.
  !!
  subroutine mpifx_allgatherv_${SUFFIX}$(mycomm, send, recv, recvcounts, displs, error)
    type(mpifx_comm), intent(in) :: mycomm
    ${TYPE}$, intent(in) :: send${RANKSUFFIX(RANK)}$
    ${TYPE}$, intent(out) :: recv${RANKSUFFIX(RANK)}$
    integer, intent(in) :: recvcounts(:)
    integer, intent(in), optional :: displs(:)
    integer, intent(out), optional :: error

    integer :: error0, ii
    integer, allocatable :: displs0(:)

    @:ASSERT(size(recv) == sum(recvcounts))
    allocate(displs0(mycomm%size))
    if (present(displs)) then
      @:ASSERT(size(displs) == mycomm%size)
      displs0 = displs
    else
      displs0(1) = 0
      do ii = 2, mycomm%size
        displs0(ii) = displs0(ii-1) + recvcounts(ii-1)
      end do
    end if

    call mpi_allgatherv(send, size(send), ${MPI_TYPE}$, recv, recvcounts, displs0, &
        & ${MPI_TYPE}$, mycomm%id, error0)

    call handle_errorflag(error0, "MPI_ALLGATHERV in mpifx_allgatherv_${SUFFIX}$", error)

  end subroutine mpifx_allgatherv_${SUFFIX}$

#:enddef

#:def mpifx_allgatherv_dr1_template(SUFFIX, TYPE, SEND_RANK, SEND_BUFFER_SIZE, RECV_RANK, MPI_TYPE)
  #!
  #!
  #! ${BUFFER_SIZE}$: send buffer size (1 or size(send))
  #! ${MPI_TYPE}$: corresponding MPI type
  #!
  !> Gathers results on one process (type ${SUFFIX}$).
  !!
  !! \param mycomm  MPI communicator.
  !! \param send  Quantity to be sent for gathering.
  !! \param recv  Received data on receive node (indefined on other nodes)
  !! \param recvcounts Counts of received data from each process
  !! \param displs Entry i specifies where to place data from process rank i-1
  !!               (default: computed from recvcounts assuming order with rank)
  !! \param error  Error code on exit.
  !!
  subroutine mpifx_allgatherv_${SUFFIX}$(mycomm, send, recv, recvcounts, displs, error)
    type(mpifx_comm), intent(in) :: mycomm
    ${TYPE}$, intent(in) :: send${RANKSUFFIX(SEND_RANK)}$
    ${TYPE}$, intent(out) :: recv${RANKSUFFIX(RECV_RANK)}$
    integer, intent(in) :: recvcounts(:)
    integer, intent(in), optional :: displs(:)
    integer, intent(out), optional :: error

    integer :: ii, error0
    integer, allocatable :: displs0(:)

    @:ASSERT(size(recv) == sum(recvcounts))
    @:ASSERT(size(recv, dim=${RECV_RANK}$) == mycomm%size)
    allocate(displs0(mycomm%size))
    if (present(displs)) then
      @:ASSERT(size(displs) == mycomm%size)
      displs0 = displs
    else
      displs0(1) = 0
      do ii = 2, mycomm%size
        displs0(ii) = displs0(ii-1) + recvcounts(ii-1)
      end do
    end if

    call mpi_allgatherv(send, ${SEND_BUFFER_SIZE}$, ${MPI_TYPE}$, recv, recvcounts, displs0, &
         & ${MPI_TYPE}$,  mycomm%id, error0)

    call handle_errorflag(error0, "MPI_ALLGATHERV in mpifx_allgatherv_${SUFFIX}$", error)

  end subroutine mpifx_allgatherv_${SUFFIX}$

#:enddef

!> Contains wrapper for \c MPI_allgatherv
module mpifx_allgatherv_module
  use mpifx_common_module
  implicit none
  private

  public :: mpifx_allgatherv

  !> Gathers scalars/arrays of different lengths on all nodes.
  !!
  !! \details All functions have the same argument list only differing in the
  !! type and rank of the second and third arguments. The second and third
  !! arguments can be of type integer (i), real (s), double precision (d),
  !! complex (c), double complex (z) and logical (l). Their rank can vary from
  !! zero (scalars) up to the maximum rank. Both arguments must be of same
  !! type. The third argument must have the size of the second times the number
  !! of processes taking part in the gathering. The fourth argument must be
  !! an array of integers corresponding to the array sizes received from each
  !! processor. The displacements at which to place the incoming data can be
  !! given as an optional argument. By default they are computed from recvcounts,
  !! assuming ordering with processor rank.
  !!
  !! \see MPI documentation (\c MPI_allgatherv)
  !!
  !! Example:
  !!
  !!     program test_allgatherv
  !!       use libmpifx_module
  !!       implicit none
  !!
  !!       type(mpifx_comm) :: mycomm
  !!       real, allocatable :: send1(:)
  !!       real, allocatable :: recv1(:)
  !!       integer, allocatable :: recvcounts(:)
  !!       integer :: ii, nrecv
  !!       character(100) :: formstr
  !!       character(*), parameter :: label = "(I2.2,'-',I3.3,'|',1X"
  !!
  !!       call mpifx_init()
  !!       call mycomm%init()
  !!
  !!       ! I1 -> I1
  !!       allocate(send1(mycomm%rank+1))
  !!       send1 = 1.0*mycomm%rank
  !!       ! recv1 size is 1+2+3+...+mycomm%size
  !!       nrecv = mycomm%size*(mycomm%size+1)/2
  !!       allocate(recv1(nrecv))
  !!       recv1(:) = 0
  !!       allocate(recvcounts(mycomm%size))
  !!       do ii = 1, mycomm%size
  !!         recvcounts(ii) = ii
  !!       end do
  !!
  !!       write(*, *) mycomm%rank, "Send1 buffer:", send1(:)
  !!       call mpifx_allgatherv(mycomm, send1, recv1, recvcounts)
  !!       if (mycomm%lead) then
  !!         write(*, *) mycomm%rank, "Recv1 buffer:", recv1
  !!       end if
  !!
  !!       call mpifx_finalize()
  !!
  !!     end program test_allgatherv
  !!
  interface mpifx_allgatherv
    #:for TYPE in TYPES
      #:for RANK in RANKS
        #:set TYPEABBREV = TYPE_ABBREVS[TYPE]
        module procedure mpifx_allgatherv_${TYPEABBREV}$${RANK}$${TYPEABBREV}$${RANK}$
      #:endfor
      module procedure mpifx_allgatherv_${TYPEABBREV}$0${TYPEABBREV}$1
    #:endfor
  end interface mpifx_allgatherv

contains

  #:for TYPE in TYPES

    #:set FTYPE = FORTRAN_TYPES[TYPE]
    #:set MPITYPE = MPI_TYPES[TYPE]

    #:for RANK in RANKS
      #:set SUFFIX = TYPE_ABBREVS[TYPE] + str(RANK) + TYPE_ABBREVS[TYPE] + str(RANK)
      $:mpifx_allgatherv_dr0_template(SUFFIX, FTYPE, RANK, MPITYPE)
    #:endfor

    #:set SUFFIX = TYPE_ABBREVS[TYPE] + str(0) + TYPE_ABBREVS[TYPE] + str(1)
    $:mpifx_allgatherv_dr1_template(SUFFIX, FTYPE, 0, 1, 1, MPITYPE)

  #:endfor

end module mpifx_allgatherv_module

ext: fpp

max_stars_repo_name: terminationshock/mpifx

max_issues_repo_name: terminationshock/mpifx

max_forks_repo_name: terminationshock/mpifx


program pi_computation

implicit none
include 'mpif.h'

      integer :: ierr, size, rank
      integer :: N, iff, ill, i
      integer, parameter :: iroot = 0
      real :: dx, x, summ, integrand, t_sum

      N = 20

      call MPI_INIT(ierr)

      call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
      call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)

      if (rank.gt.iroot) then

                dx = 1.d0/N
                call decompose(N, size, rank, iff, ill)
                summ = 0
                do i = iff, ill
                        x = (real(i) - 0.5)*dx
                        summ = summ + integrand(x)*dx
                end do
      else
                summ = 0
      end if

      call MPI_REDUCE(summ,t_sum, 1, MPI_REAL, MPI_SUM, iroot, MPI_COMM_WORLD, ierr)

      if(rank.eq.iroot) then
                print *, 'The computed pi value:', t_sum
      end if

      call MPI_FINALIZE(ierr)

end program

!-------------Subroutine decompose ----------------------

subroutine decompose(n, size, rank, iff, ill)

        implicit none
        integer :: n, size, rank, iff, ill, neven

        neven = n/(size-1)
        iff = (rank-1)*neven + 1
        ill = rank*neven

        if (rank.eq.size-1) then
                ill = rank*neven + mod(n, size-1)
        end if

end subroutine
!--------------------------------------------------------
function integrand(x)

        implicit none
        real :: x, integrand

        integrand = 4.d0/(1.d0 + x**2)

end function


text:
      subroutine grdfcn(n,x,g,nprob)
      integer n,nprob
      double precision x(n),g(n)
c     **********
c
c     subroutine grdfcn
c
c     this subroutine defines the gradient vectors of eighteen
c     nonlinear unconstrained minimization problems. the problem
c     dimensions are as described in the prologue comments of objfcn.
c
c     the subroutine statement is
c
c       subroutine grdfcn(n,x,g,nprob)
c
c     where
c
c       n is a positive integer input variable.
c
c       x is an input array of length n.
c
c       g is an output array of length n which contains the components
c         of the gradient vector of the nprob objective function
c         evaluated at x.
c
c       nprob is a positive integer input variable which defines the
c         number of the problem. nprob must not exceed 18.
c
c     subprograms called
c
c       fortran-supplied ... dabs,datan,dcos,dexp,dlog,dsign,dsin,
c                            dsqrt
c
c     argonne national laboratory. minpack project. march 1980.
c     burton s. garbow, kenneth e. hillstrom, jorge j. more
c
c     **********
      integer i,iev,ivar,j
      double precision ap,arg,c2pdm6,cp0001,cp1,cp2,cp25,cp5,c1p5,
     *                 c2p25,c2p625,c3p5,c19p8,c20p2,c25,c29,c100,
     *                 c180,c200,c10000,c1pd6,d1,d2,eight,fifty,five,
     *                 four,one,r,s1,s2,s3,t,t1,t2,t3,ten,th,three,
     *                 tpi,twenty,two,zero
      double precision fvec(50),y(15)
      double precision dfloat
      data zero,one,two,three,four,five,eight,ten,twenty,fifty
     *     /0.0d0,1.0d0,2.0d0,3.0d0,4.0d0,5.0d0,8.0d0,1.0d1,2.0d1,
     *      5.0d1/
      data c2pdm6,cp0001,cp1,cp2,cp25,cp5,c1p5,c2p25,c2p625,c3p5,
     *     c19p8,c20p2,c25,c29,c100,c180,c200,c10000,c1pd6
     *     /2.0d-6,1.0d-4,1.0d-1,2.0d-1,2.5d-1,5.0d-1,1.5d0,2.25d0,
     *      2.625d0,3.5d0,1.98d1,2.02d1,2.5d1,2.9d1,1.0d2,1.8d2,2.0d2,
     *      1.0d4,1.0d6/
      data ap /1.0d-5/
      data y(1),y(2),y(3),y(4),y(5),y(6),y(7),y(8),y(9),y(10),y(11),
     *     y(12),y(13),y(14),y(15)
     *     /9.0d-4,4.4d-3,1.75d-2,5.4d-2,1.295d-1,2.42d-1,3.521d-1,
     *      3.989d-1,3.521d-1,2.42d-1,1.295d-1,5.4d-2,1.75d-2,4.4d-3,
     *      9.0d-4/
      dfloat(ivar) = ivar
c
c     gradient routine selector.
c
      go to (10,20,50,70,80,100,130,190,220,260,270,290,310,350,370,
     *       390,400,410), nprob
c
c     helical valley function.
c
   10 continue
      tpi = eight*datan(one)
      th = dsign(cp25,x(2))
      if (x(1) .gt. zero) th = datan(x(2)/x(1))/tpi
      if (x(1) .lt. zero) th = datan(x(2)/x(1))/tpi + cp5
      arg = x(1)**2 + x(2)**2
      r = dsqrt(arg)
      t = x(3) - ten*th
      s1 = ten*t/(tpi*arg)
      g(1) = c200*(x(1) - x(1)/r + x(2)*s1)
      g(2) = c200*(x(2) - x(2)/r - x(1)*s1)
      g(3) = two*(c100*t + x(3))
      go to 490
c
c     biggs exp6 function.
c
   20 continue
      do 30 j = 1, 6
         g(j) = zero
   30    continue
      do 40 i = 1, 13
         d1 = dfloat(i)/ten
         d2 = dexp(-d1) - five*dexp(-ten*d1) + three*dexp(-four*d1)
         s1 = dexp(-d1*x(1))
         s2 = dexp(-d1*x(2))
         s3 = dexp(-d1*x(5))
         t = x(3)*s1 - x(4)*s2 + x(6)*s3 - d2
         th = d1*t
         g(1) = g(1) - s1*th
         g(2) = g(2) + s2*th
         g(3) = g(3) + s1*t
         g(4) = g(4) - s2*t
         g(5) = g(5) - s3*th
         g(6) = g(6) + s3*t
   40    continue
      g(1) = two*x(3)*g(1)
      g(2) = two*x(4)*g(2)
      g(3) = two*g(3)
      g(4) = two*g(4)
      g(5) = two*x(6)*g(5)
      g(6) = two*g(6)
      go to 490
c
c     gaussian function.
c
   50 continue
      g(1) = zero
      g(2) = zero
      g(3) = zero
      do 60 i = 1, 15
         d1 = cp5*dfloat(i-1)
         d2 = c3p5 - d1 - x(3)
         arg = -cp5*x(2)*d2**2
         r = dexp(arg)
         t = x(1)*r - y(i)
         s1 = r*t
         s2 = d2*s1
         g(1) = g(1) + s1
         g(2) = g(2) - d2*s2
         g(3) = g(3) + s2
   60    continue
      g(1) = two*g(1)
      g(2) = x(1)*g(2)
      g(3) = two*x(1)*x(2)*g(3)
      go to 490
c
c     powell badly scaled function.
c
   70 continue
      t1 = c10000*x(1)*x(2) - one
      s1 = dexp(-x(1))
      s2 = dexp(-x(2))
      t2 = s1 + s2 - one - cp0001
      g(1) = two*(c10000*x(2)*t1 - s1*t2)
      g(2) = two*(c10000*x(1)*t1 - s2*t2)
      go to 490
c
c     box 3-dimensional function.
c
   80 continue
      g(1) = zero
      g(2) = zero
      g(3) = zero
      do 90 i = 1, 10
         d1 = dfloat(i)
         d2 = d1/ten
         s1 = dexp(-d2*x(1))
         s2 = dexp(-d2*x(2))
         s3 = dexp(-d2) - dexp(-d1)
         t = s1 - s2 - s3*x(3)
         th = d2*t
         g(1) = g(1) - s1*th
         g(2) = g(2) + s2*th
         g(3) = g(3) - s3*t
   90    continue
      g(1) = two*g(1)
      g(2) = two*g(2)
      g(3) = two*g(3)
      go to 490
c
c     variably dimensioned function.
c
  100 continue
      t1 = zero
      do 110 j = 1, n
         t1 = t1 + dfloat(j)*(x(j) - one)
  110    continue
      t = t1*(one + two*t1**2)
      do 120 j = 1, n
         g(j) = two*(x(j) - one + dfloat(j)*t)
  120    continue
      go to 490
c
c     watson function.
c
  130 continue
      do 140 j = 1, n
         g(j) = zero
  140    continue
      do 180 i = 1, 29
         d1 = dfloat(i)/c29
         s1 = zero
         d2 = one
         do 150 j = 2, n
            s1 = s1 + dfloat(j-1)*d2*x(j)
            d2 = d1*d2
  150       continue
         s2 = zero
         d2 = one
         do 160 j = 1, n
            s2 = s2 + d2*x(j)
            d2 = d1*d2
  160       continue
         t = s1 - s2**2 - one
         s3 = two*d1*s2
         d2 = two/d1
         do 170 j = 1, n
            g(j) = g(j) + d2*(dfloat(j-1) - s3)*t
            d2 = d1*d2
  170       continue
  180    continue
      t1 = x(2) - x(1)**2 - one
      g(1) = g(1) + x(1)*(two - four*t1)
      g(2) = g(2) + two*t1
      go to 490
c
c     penalty function i.
c
  190 continue
      t1 = -cp25
      do 200 j = 1, n
         t1 = t1 + x(j)**2
  200    continue
      d1 = two*ap
      th = four*t1
      do 210 j = 1, n
         g(j) = d1*(x(j) - one) + x(j)*th
  210    continue
      go to 490
c
c     penalty function ii.
c
  220 continue
      t1 = -one
      do 230 j = 1, n
         t1 = t1 + dfloat(n-j+1)*x(j)**2
  230    continue
      d1 = dexp(cp1)
      d2 = one
      th = four*t1
      do 250 j = 1, n
         g(j) = dfloat(n-j+1)*x(j)*th
         s1 = dexp(x(j)/ten)
         if (j .eq. 1) go to 240
         s3 = s1 + s2 - d2*(d1 + one)
         g(j) = g(j) + ap*s1*(s3 + s1 - one/d1)/five
         g(j-1) = g(j-1) + ap*s2*s3/five
  240    continue
         s2 = s1
         d2 = d1*d2
  250    continue
      g(1) = g(1) + two*(x(1) - cp2)
      go to 490
c
c     brown badly scaled function.
c
  260 continue
      t1 = x(1) - c1pd6
      t2 = x(2) - c2pdm6
      t3 = x(1)*x(2) - two
      g(1) = two*(t1 + x(2)*t3)
      g(2) = two*(t2 + x(1)*t3)
      go to 490
c
c     brown and dennis function.
c
  270 continue
      g(1) = zero
      g(2) = zero
      g(3) = zero
      g(4) = zero
      do 280 i = 1, 20
         d1 = dfloat(i)/five
         d2 = dsin(d1)
         t1 = x(1) + d1*x(2) - dexp(d1)
         t2 = x(3) + d2*x(4) - dcos(d1)
         t = t1**2 + t2**2
         s1 = t1*t
         s2 = t2*t
         g(1) = g(1) + s1
         g(2) = g(2) + d1*s1
         g(3) = g(3) + s2
         g(4) = g(4) + d2*s2
  280    continue
      g(1) = four*g(1)
      g(2) = four*g(2)
      g(3) = four*g(3)
      g(4) = four*g(4)
      go to 490
c
c     gulf research and development function.
c
  290 continue
      g(1) = zero
      g(2) = zero
      g(3) = zero
      d1 = two/three
      do 300 i = 1, 99
         arg = dfloat(i)/c100
         r = (-fifty*dlog(arg))**d1 + c25 - x(2)
         t1 = dabs(r)**x(3)/x(1)
         t2 = dexp(-t1)
         t = t2 - arg
         s1 = t1*t2*t
         g(1) = g(1) + s1
         g(2) = g(2) + s1/r
         g(3) = g(3) - s1*dlog(dabs(r))
  300    continue
      g(1) = two*g(1)/x(1)
      g(2) = two*x(3)*g(2)
      g(3) = two*g(3)
      go to 490
c
c     trigonometric function.
c
  310 continue
      s1 = zero
      do 320 j = 1, n
         g(j) = dcos(x(j))
         s1 = s1 + g(j)
  320    continue
      s2 = zero
      do 330 j = 1, n
         th = dsin(x(j))
         t = dfloat(n+j) - th - s1 - dfloat(j)*g(j)
         s2 = s2 + t
         g(j) = (dfloat(j)*th - g(j))*t
  330    continue
      do 340 j = 1, n
         g(j) = two*(g(j) + dsin(x(j))*s2)
  340    continue
      go to 490
c
c     extended rosenbrock function.
c
  350 continue
      do 360 j = 1, n, 2
         t1 = one - x(j)
         g(j+1) = c200*(x(j+1) - x(j)**2)
         g(j) = -two*(x(j)*g(j+1) + t1)
  360    continue
      go to 490
c
c     extended powell function.
c
  370 continue
      do 380 j = 1, n, 4
         t = x(j) + ten*x(j+1)
         t1 = x(j+2) - x(j+3)
         s1 = five*t1
         t2 = x(j+1) - two*x(j+2)
         s2 = four*t2**3
         t3 = x(j) - x(j+3)
         s3 = twenty*t3**3
         g(j) = two*(t + s3)
         g(j+1) = twenty*t + s2
         g(j+2) = two*(s1 - s2)
         g(j+3) = -two*(s1 + s3)
  380    continue
      go to 490
c
c     beale function.
c
  390 continue
      s1 = one - x(2)
      t1 = c1p5 - x(1)*s1
      s2 = one - x(2)**2
      t2 = c2p25 - x(1)*s2
      s3 = one - x(2)**3
      t3 = c2p625 - x(1)*s3
      g(1) = -two*(s1*t1 + s2*t2 + s3*t3)
      g(2) = two*x(1)*(t1 + x(2)*(two*t2 + three*x(2)*t3))
      go to 490
c
c     wood function.
c
  400 continue
      s1 = x(2) - x(1)**2
      s2 = one - x(1)
      s3 = x(2) - one
      t1 = x(4) - x(3)**2
      t2 = one - x(3)
      t3 = x(4) - one
      g(1) = -two*(c200*x(1)*s1 + s2)
      g(2) = c200*s1 + c20p2*s3 + c19p8*t3
      g(3) = -two*(c180*x(3)*t1 + t2)
      g(4) = c180*t1 + c20p2*t3 + c19p8*s3
      go to 490
c
c     chebyquad function.
c
  410 continue
      do 420 i = 1, n
         fvec(i) = zero
  420    continue
      do 440 j = 1, n
         t1 = one
         t2 = two*x(j) - one
         t = two*t2
         do 430 i = 1, n
            fvec(i) = fvec(i) + t2
            th = t*t2 - t1
            t1 = t2
            t2 = th
  430       continue
  440    continue
      d1 = one/dfloat(n)
      iev = -1
      do 450 i = 1, n
         fvec(i) = d1*fvec(i)
         if (iev .gt. 0) fvec(i) = fvec(i) + one/(dfloat(i)**2 - one)
         iev = -iev
  450    continue
      do 470 j = 1, n
         g(j) = zero
         t1 = one
         t2 = two*x(j) - one
         t = two*t2
         s1 = zero
         s2 = two
         do 460 i = 1, n
            g(j) = g(j) + fvec(i)*s2
            th = four*t2 + t*s2 - s1
            s1 = s2
            s2 = th
            th = t*t2 - t1
            t1 = t2
            t2 = th
  460       continue
  470    continue
      d2 = two*d1
      do 480 j = 1, n
         g(j) = d2*g(j)
  480    continue
  490 continue
      return
c
c     last card of subroutine grdfcn.
c
      end

ext: f

max_stars_repo_name: ianormy/Pyto

max_issues_repo_name: ProjectZeroDays/Pyto

max_forks_repo_name: ProjectZeroDays/Pyto


text:
! ALOHA global configuration parameters
module aloha_config
    USE aloha_constants, ONLY : wp

    ! Vacuum wavenumber
    real(kind=wp)       :: k0

    ! waveguide parameters
    character, dimension(1) :: bcte
    real(kind=wp)       :: a
    real(kind=wp)       :: b(100), z(100)

    ! ?
    real(kind=wp)       :: epsabs
    ! plasma ??
    real(kind=wp)       :: X0, D0
    real(kind=wp)       :: d_vide !V5
    real(kind=wp)       :: X1, D1, d_couche, pertes !V6
    !modmn/
    integer             :: m, n
    ! guid
    integer             :: i,j

    ! ########## numerical integration configuration
    integer :: nlimit
    real(kind=wp) ::  max_nz

    ! Mode number
    integer :: Nmhm, Nmem

    integer :: Nmax=10 ! nombre de mode max
    integer :: Gmax=60 ! nombre de guide max

    ! Airy
    integer :: app

    ! telnum
    integer :: knout

end module aloha_config

max_stars_repo_name: IRFM/ALOHA

max_issues_repo_name: IRFM/ALOHA

max_forks_repo_name: IRFM/ALOHA


text:
! { dg-do compile }
! { dg-options "-Wall" }

! This long comment line should not trigger a line-truncation warning with -Wall even for free-form 132 character line limit (blah blah)

      PROGRAM foo
      WRITE (*,*) "Test"           ! Neither this comment which exceeds the 132 character limit with some random words, too (blah blah)
      WRITE (*,*) "This exactly 132 character long soruce line not, too.  How can people fill 132 characters without sensless stuff"
      END

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


      subroutine ch_read_elements

      use input_file_module
      use maximum_data_module
      use calibration_data_module
      use hydrograph_module
      use sd_channel_module

      implicit none

      character (len=80) :: titldum    !          |title of file
      character (len=80) :: header     !          |header of file
      integer :: eof                   !          |end of file
      integer :: imax                  !          |determine max number for array (imax) and total number in file
      integer :: mcal                  !units     |description
      logical :: i_exist               !          |check to determine if file exists
      integer :: mreg                  !units     |description
      integer :: i                     !none      |counter
      integer :: k                     !units     |description
      integer :: nspu                  !units     |description
      integer :: isp                   !none      |counter
      integer :: ielem1                !units     |description
      integer :: ii                    !none      |counter
      integer :: icha                  !none      |counter
      integer :: ireg                  !none      |counter
      integer :: ires                  !none      |counter

      imax = 0
      mcal = 0
      mreg = 0

    inquire (file=in_regs%def_cha, exist=i_exist)
    if (i_exist .or. in_regs%def_cha /= "null") then
      do
        open (107,file=in_regs%def_cha)
        read (107,*,iostat=eof) titldum
        if (eof < 0) exit
        read (107,*,iostat=eof) mreg
        if (eof < 0) exit
        read (107,*,iostat=eof) header

        !! allocate channel outputs for regions
        allocate (schsd_d(0:mreg)); allocate (schsd_m(0:mreg)); allocate (schsd_y(0:mreg)); allocate (schsd_a(0:mreg))

      do i = 1, mreg

        read (107,*,iostat=eof) k, ccu_out(i)%name, ccu_out(i)%area_ha, nspu
        if (eof < 0) exit
        if (nspu > 0) then
          allocate (elem_cnt(nspu))
          backspace (107)
          read (107,*,iostat=eof) k, ccu_out(i)%name, ccu_out(i)%area_ha, nspu, (elem_cnt(isp), isp = 1, nspu)
          if (eof < 0) exit

          call define_unit_elements (nspu, ielem1)

          allocate (ccu_out(i)%num(ielem1))
          ccu_out(i)%num = defunit_num
          ccu_out(i)%num_tot = ielem1
          deallocate (defunit_num)
        else
          !!all hrus are in region
          allocate (ccu_out(i)%num(sp_ob%hru))
          ccu_out(i)%num_tot = sp_ob%res
          do icha = 1, sp_ob%chan
            ccu_out(i)%num(ires) = icha
          end do
        end if

      end do    ! i = 1, mreg
      exit

      db_mx%cha_reg = mreg
      end do
      end if

    !! setting up regions for channel soft cal and/or output by order
    inquire (file=in_regs%def_cha_reg, exist=i_exist)
    if (i_exist .or. in_regs%def_cha_reg /= "null") then
      do
        open (107,file=in_regs%def_cha_reg)
        read (107,*,iostat=eof) titldum
        if (eof < 0) exit
        read (107,*,iostat=eof) mreg
        if (eof < 0) exit
        read (107,*,iostat=eof) header
        if (eof < 0) exit

      do i = 1, mreg

        read (107,*,iostat=eof) k, ccu_reg(i)%name, ccu_reg(i)%area_ha, nspu
        if (eof < 0) exit
        if (nspu > 0) then
          allocate (elem_cnt(nspu))
          backspace (107)
          read (107,*,iostat=eof) k, ccu_reg(i)%name, ccu_reg(i)%area_ha, nspu, (elem_cnt(isp), isp = 1, nspu)
          if (eof < 0) exit

          call define_unit_elements (nspu, ielem1)

          allocate (ccu_reg(i)%num(ielem1))
          ccu_reg(i)%num = defunit_num
          ccu_reg(i)%num_tot = ielem1

        else
          !!all channels are in region
          allocate (ccu_reg(i)%num(sp_ob%hru))
          ccu_reg(i)%num_tot = sp_ob%chan
          do icha = 1, sp_ob%chan
            ccu_reg(i)%num(icha) = icha
          end do
        end if

      end do    ! i = 1, mreg
      exit

      db_mx%cha_reg = mreg
      end do
      end if

      !! if no regions are input, don"t need elements
      if (mreg > 0) then

      do ireg = 1, mreg
        ccu_cal(ireg)%lum_ha_tot = 0.
        ccu_cal(ireg)%lum_num_tot = 0
        ccu_cal(ireg)%lum_ha_tot = 0.
        !allocate (region(ireg)%lum_ha_tot(db_mx%landuse))
        !allocate (region(ireg)%lum_num_tot(db_mx%landuse))
        !allocate (rwb_a(ireg)%lum(db_mx%landuse))
        !allocate (rnb_a(ireg)%lum(db_mx%landuse))
        !allocate (rls_a(ireg)%lum(db_mx%landuse))
        !allocate (rpw_a(ireg)%lum(db_mx%landuse))
      end do
      end if    ! mreg > 0

      !!read data for each element in all landscape cataloging units
      inquire (file="element.ccu", exist=i_exist)
      if (i_exist ) then
      do
        open (107,file="element.ccu")
        read (107,*,iostat=eof) titldum
        if (eof < 0) exit
        read (107,*,iostat=eof) header
        if (eof < 0) exit
        imax = 0
          do while (eof == 0)
              read (107,*,iostat=eof) i
              if (eof < 0) exit
              imax = Max(i,imax)
          end do

        allocate (ccu_elem(imax))

        rewind (107)
        read (107,*,iostat=eof) titldum
        if (eof < 0) exit
        read (107,*,iostat=eof) header
        if (eof < 0) exit

        do isp = 1, imax
          read (107,*,iostat=eof) i
          if (eof < 0) exit
          backspace (107)
          read (107,*,iostat=eof) k, ccu_elem(i)%name, ccu_elem(i)%obtyp, ccu_elem(i)%obtypno,      &
                                    ccu_elem(i)%bsn_frac, ccu_elem(i)%ru_frac, ccu_elem(i)%reg_frac
          if (eof < 0) exit
        end do
        exit
      end do
      end if

      ! set hru number from element number and set hru areas in the region
      do ireg = 1, mreg
        do icha = 1, ccu_reg(ireg)%num_tot      !elements have to be hru or hru_lte
          ielem1 = ccu_reg(ireg)%num(icha)
          !switch %num from element number to hru number
          ccu_cal(ireg)%num(icha) = ccu_elem(ielem1)%obtypno
          ccu_cal(ireg)%hru_ha(icha) = ccu_elem(ielem1)%ru_frac * ccu_cal(ireg)%area_ha
        end do
      end do

      close (107)

      return
      end subroutine ch_read_elements


text:
*DECK DPOFA
      SUBROUTINE DPOFA (A, LDA, N, INFO)
C***BEGIN PROLOGUE  DPOFA
C***PURPOSE  Factor a real symmetric positive definite matrix.
C***LIBRARY   SLATEC (LINPACK)
C***CATEGORY  D2B1B
C***TYPE      DOUBLE PRECISION (SPOFA-S, DPOFA-D, CPOFA-C)
C***KEYWORDS  LINEAR ALGEBRA, LINPACK, MATRIX FACTORIZATION,
C             POSITIVE DEFINITE
C***AUTHOR  Moler, C. B., (U. of New Mexico)
C***DESCRIPTION
C
C     DPOFA factors a double precision symmetric positive definite
C     matrix.
C
C     DPOFA is usually called by DPOCO, but it can be called
C     directly with a saving in time if  RCOND  is not needed.
C     (time for DPOCO) = (1 + 18/N)*(time for DPOFA) .
C
C     On Entry
C
C        A       DOUBLE PRECISION(LDA, N)
C                the symmetric matrix to be factored.  Only the
C                diagonal and upper triangle are used.
C
C        LDA     INTEGER
C                the leading dimension of the array  A .
C
C        N       INTEGER
C                the order of the matrix  A .
C
C     On Return
C
C        A       an upper triangular matrix  R  so that  A = TRANS(R)*R
C                where  TRANS(R)  is the transpose.
C                The strict lower triangle is unaltered.
C                If  INFO .NE. 0 , the factorization is not complete.
C
C        INFO    INTEGER
C                = 0  for normal return.
C                = K  signals an error condition.  The leading minor
C                     of order  K  is not positive definite.
C
C***REFERENCES  J. J. Dongarra, J. R. Bunch, C. B. Moler, and G. W.
C                 Stewart, LINPACK Users' Guide, SIAM, 1979.
C***ROUTINES CALLED  DDOT
C***REVISION HISTORY  (YYMMDD)
C   780814  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890831  Modified array declarations.  (WRB)
C   890831  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DPOFA
      INTEGER LDA,N,INFO
      DOUBLE PRECISION A(LDA,*)
C
      DOUBLE PRECISION DDOT,T
      DOUBLE PRECISION S
      INTEGER J,JM1,K
C***FIRST EXECUTABLE STATEMENT  DPOFA
         DO 30 J = 1, N
            INFO = J
            S = 0.0D0
            JM1 = J - 1
            IF (JM1 .LT. 1) GO TO 20
            DO 10 K = 1, JM1
               T = A(K,J) - DDOT(K-1,A(1,K),1,A(1,J),1)
               T = T/A(K,K)
               A(K,J) = T
               S = S + T*T
   10       CONTINUE
   20       CONTINUE
            S = A(J,J) - S
            IF (S .LE. 0.0D0) GO TO 40
            A(J,J) = SQRT(S)
   30    CONTINUE
         INFO = 0
   40 CONTINUE
      RETURN
      END


text:
      subroutine parse_blocks(wlist_p,wlist,luerr)
      ! sort wlist into levels according to blocks indicated
      ! by the separators ( and )

      implicit none

      include 'def_word_list.h'

      integer, parameter ::
     &      ntest = 00

      character, parameter ::
     &      char_block_open  = '(',
     &      char_block_close = ')',
     &      char_comnt_1 = '#',
     &      char_comnt_2 = '!'

      type(word_list), intent(inout) ::
     &      wlist
      type(word_list), intent(inout) ::
     &      wlist_p
      integer, intent(in) ::
     &      luerr

      type(word_list_entry), pointer ::
     &      wl_pnt, wlp_pnt

      integer ::
     &      level
      logical ::
     &      branch
      character ::
     &      sep_rem

      ! set to head
      wlist%current   => wlist%head
      wlist_p%current => wlist_p%head

      branch = .false.
      level = 0
      main_lp: do

        ! check whether this is a comment
        if (wlist%current%word(1:1).eq.char_comnt_1.or.
     &      wlist%current%word(1:1).eq.char_comnt_2) then

          ! skip until newline appears (marked by 'E')
          comment_lp: do
            sep_rem = wlist%current%sep
            if (.not.associated(wlist%current%next)) exit main_lp
            wlist%current => wlist%current%next
            if (sep_rem.eq.'E') cycle main_lp
          end do comment_lp

        end if

        ! remove empty lines and empty entries that indicate
        ! ending lines after other seprators
        if (.not.(len_trim(wlist%current%word).eq.0.and.
     &                     wlist%current%sep.eq.'E')) then

          call new_word_list_entry(wlist_p,branch)
          ! copy entry
          wlist_p%current%word = wlist%current%word
          wlist_p%current%sep  = wlist%current%sep
          wlist_p%current%line = wlist%current%line
          wlist_p%current%col  = wlist%current%col

          ! open new block?
          if (wlist%current%sep.eq.char_block_open) then
            branch = .true.
            level = level+1
          else if (wlist%current%sep.eq.char_block_close) then
            branch = .false.
            level = level-1
            ! move parsed list pointer one level up
            wlist_p%current => wlist_p%current%up
          else
            branch = .false.
          end if

          if (level.lt.0) exit main_lp ! error!

        end if

        ! advance original linear list
        if (.not.associated(wlist%current%next)) exit main_lp
        wlist%current => wlist%current%next

      end do main_lp

      if (level.ne.0) then
        write(luerr,*) 'Error while parsing blocks:'
        if (level.lt.0)
     &    write(luerr,*) 'A ")" occurred without matching "("'
        if (level.gt.0)
     &    write(luerr,*) 'It seems that not all "(" are matched by ")"'
        write(luerr,'(x,a,i4,a,i3)') 'I noted this in line',
     &         wlist%current%line,' near column ',wlist%current%col
        if (associated(wlist%head).and.associated(wlist%head%next))
     &      write(luerr,*) 'Current context starts with "',
     &         trim(wlist%head%word),' ',
     &         trim(wlist%head%next%word),' ... "'
        if (ntest.ge.100) call print_word_list(luerr,wlist)
        call quit(0,'parse_blocks','Error while parsing blocks')
      end if

      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      program hpftest37
      real tc(250,250), ts(250,250), north(250)
CHPF$ TEMPLATE t(250,250)
CHPF$ ALIGN tc(i,j), ts(i,j) with t(i,j)
CHPF$ ALIGN north(i) with t(1,i)
CHPF$ PROCESSORS p(4)
CHPF$ DISTRIBUTE t(block, *) ONTO p
      print *, 'hpftest37 running'
      print *, 'thermo'
c
c initialization
      print *, 'initializing'
chpf$ independent(i)
      do i=1,250
         north(i) = 250.0
      enddo
chpf$ independent(i)
      do i=1,250
         tc(1,i) = north(i)
         ts(1,i) = north(i)
      enddo
chpf$ independent(i,j)
      do i=2,250
         do j=1,250
            tc(i,j) = 10.0
            ts(i,j) = 10.0
         enddo
      enddo
c
c iterations... should be a test, but reductions are not accepted
      print *, 'running'
      do k=1,500
c
c computation and copy back (dataparallel semantic)
chpf$ independent(i,j)
         do i=2,249
            do j=2,249
               ts(i,j) = 0.25*
     $              (tc(i-1,j) + tc(i+1,j) + tc(i,j-1) + tc(i,j+1))
            enddo
         enddo
chpf$ independent(i,j)
         do i=2,249
            do j=2,249
               tc(i,j) = ts(i,j)
            enddo
         enddo
      enddo
c
c print results
      print *, 'results:'
 10   format(F8.2, F8.2, F8.2, F8.2, F8.2, F8.2, F8.2, F8.2)
      do i=25,200,25
         write (6,10) tc(i,50), tc(i,75), tc(i,100), tc(i,125),
     $        tc(i,150), tc(i,175), tc(i,200), tc(i,225)
      enddo
      print *, 'hpftest37 ended'
      end


text:
      subroutine pvswitch4(q1,F0,F1,F2,F3,F4,F4a)
C----from F0,F1,F2,F3,F4 calculate F4a which has all denominators
C----shifted by q1
      implicit none
      include 'TRydef.f'
      integer n1,n2,n3,n4,ep
      double complex F0(-2:0),F1(y1max,-2:0),F2(y2max,-2:0),
     . F3(y3max,-2:0),F4(y4max,-2:0),F4a(y4max,-2:0)
      double precision q1(4)
      do ep=-2,0
      do n1=1,4
      do n2=n1,4
      do n3=n2,4
      do n4=n3,4
      F4a(y4(n1,n2,n3,n4),ep)=F4(y4(n1,n2,n3,n4),ep)
     . -q1(n1)*F3(y3(n2,n3,n4),ep)
     . -q1(n2)*F3(y3(n3,n4,n1),ep)
     . -q1(n3)*F3(y3(n4,n1,n2),ep)
     . -q1(n4)*F3(y3(n1,n2,n3),ep)
     . +q1(n1)*q1(n2)*F2(y2(n3,n4),ep)
     . +q1(n1)*q1(n3)*F2(y2(n2,n4),ep)
     . +q1(n1)*q1(n4)*F2(y2(n2,n3),ep)
     . +q1(n2)*q1(n3)*F2(y2(n1,n4),ep)
     . +q1(n2)*q1(n4)*F2(y2(n1,n3),ep)
     . +q1(n3)*q1(n4)*F2(y2(n1,n2),ep)

     . -q1(n1)*q1(n2)*q1(n3)*F1(n4,ep)
     . -q1(n2)*q1(n3)*q1(n4)*F1(n1,ep)
     . -q1(n3)*q1(n4)*q1(n1)*F1(n2,ep)
     . -q1(n4)*q1(n1)*q1(n2)*F1(n3,ep)

     . +q1(n1)*q1(n2)*q1(n3)*q1(n4)*F0(ep)

      enddo
      enddo
      enddo
      enddo
      enddo
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      module ifbc_def_m
c
        use global_const_m
        use conpar_m
        use if_velocity_m
c
        implicit none
c
        integer, parameter :: nifbc = 1
        integer, parameter :: ivapor_frac = 1
        integer, parameter :: iwt = 2
c
        real*8, pointer :: if_normal(:,:)
        real*8, pointer :: ifbc(:,:)
c
      end module ifbc_def_m
c
      module ifbc_m
c
        use ifbc_def_m
        use workfc_m
        use number_def_m
        use pointer_data
c
        implicit none
c
        contains
c
        subroutine ifbc_malloc
          allocate(if_normal(nshg,nsd))
          allocate(ifbc(nshg,nifbc+1))
          ifbc(:,ivapor_frac) = zero
          ifbc(:,iwt) = one
        end subroutine ifbc_malloc
c
        subroutine ifbc_mfree
          deallocate (if_normal)
          deallocate (ifbc)
        end subroutine ifbc_mfree
c
        subroutine ifbc_set
     &  (
     &    y
     &,   BC
     &,   iBC
     &,   ilwork
     &,   nlwork
     &  )
c
          use elmpar_m
          use blkdat_m
          use shpdat_m
          use genpar_m
          use matdat_def_m
c
          include "mpif.h"
c
          real*8, intent(inout) :: y(nshg,ndof)
          real*8,  intent(inout) ::  BC(nshg,ndofbc)
          integer, intent(inout) :: iBC(nshg)
          integer, intent(in)    :: ilwork(nlwork)
          integer, intent(in)    :: nlwork
c
          integer :: i0,i1,n,iel,iblk,materif0,materif1, npro, ierr, isclr
          integer, pointer :: ienif0(:,:),ienif1(:,:)
          logical :: is_ideal_gas,is_ideal_gas_mixture
c
          if (numpe > 1) then
            call commu (ifbc(:,ivapor_frac), ilwork, 1, 'in ')
            call commu (ifbc(:,iwt), ilwork, 1, 'in ')
            call MPI_BARRIER (MPI_COMM_WORLD,ierr)
            call commu (ifbc(:,ivapor_frac), ilwork, 1, 'out')
            call commu (ifbc(:,iwt), ilwork, 1, 'out ')
            call MPI_BARRIER (MPI_COMM_WORLD,ierr)
          endif
c
          where(ifbc(:,ivapor_frac) .ne. zero)
            BC(:,6+isclr) = ifbc(:,ivapor_frac)/ifbc(:,iwt)
          endwhere
c
          return
c
          do iblk = 1,nelblif
c
            npro  = lcblkif(1,iblk+1) - lcblkif(1,iblk)
c
            ienif0 => mienif0(iblk)%p
            ienif1 => mienif1(iblk)%p
c
            do iel = 1,npro
c              do n = 1,nshl0
              do n = 1,3  ! only triangles on the surface
                i0 = ienif0(iel,n)
                i1 = ienif1(iel,n)
c
c... set vapor fraction
c
      isclr = 1
                BC(i0,6+isclr) = ifbc(i0,ivapor_frac)/ifbc(i0,iwt)
C                BC(i0,6+isclr) = zero
                BC(i1,6+isclr) = one
                y(i0,5+isclr) = BC(i0,6+isclr)
                y(i1,5+isclr) = BC(i1,6+isclr)
c      write(*,*) 'i0, i1: ',i0,i1,iBC(i0),iBC(i1),btest(ibc(i0),6),btest(ibc(i1),6)
C      write(*,*) '[',myrank,'] i0: ',i0,ifbc(i0,ivapor_frac),ifbc(i0,iwt),BC(i0,6+isclr),btest(ibc(i0),6)
c      write(*,*) 'i1: ',i1,ifbc(i1,ivapor_frac),ifbc(i1,iwt),BC(i1,6+isclr),btest(ibc(i1),6)
c
              enddo
            enddo
          enddo
c
        end subroutine ifbc_set
c
      end module ifbc_m

max_stars_repo_name: SCOREC/phasta

max_issues_repo_name: SCOREC/phasta

max_forks_repo_name: SCOREC/phasta


text:
program potencial
    use iso_fortran_env

    implicit none

    integer, parameter   :: n = 100
    real*8, dimension(:), allocatable   :: x, y, z, m, ep
    real*8, parameter      :: G = 4.299e-6, eps=0.0
    real*8                 :: dist
    integer                :: i, j

    allocate(x(n), y(n), z(n), m(n), ep(n))

    !-------------leo las posiciones de las particulas---------------------
    open(20, file='../tests/test_data/mock_particles.dat', status='old')
    do i=1, n
            read(20,*) x(i), y(i), z(i), m(i)
    end do
     close(20)

    !======================================================================
    !Calcula la energia potencial específica de cada una de las particulas

    do i = 1, n
            ep(i) = 0.
            do j = 1, n
                    if (i /= j) then

                            dist = sqrt((x(i)-x(j))**2 + (y(i)-y(j))**2 + (z(i)-z(j))**2 + eps**2)
                            ep(i) = ep(i) + m(j)/dist

                    end if
            end do
            ep(i) = ep(i)*G
    end do
    !$OMP END DO
    !$OMP END PARALLEL

    !======================================================================

    !-------------escribo los potenciales de las particulas---------------------
    open(30, file='../tests/test_data/fpotential_test.dat', status='unknown')
    do i = 1,n
            write(30,*) ep(i)
    end do
     close(30)

    deallocate(x, y, z, m, ep)

    end program

max_stars_repo_name: vcristiani/galaxy-chop

max_issues_repo_name: vcristiani/galaxy-chop

max_forks_repo_name: vcristiani/galaxy-chop


text:
C$Procedure      SPKW08 ( Write SPK segment, type 8 )

      SUBROUTINE SPKW08 (  HANDLE,  BODY,    CENTER,  FRAME,
     .                     FIRST,   LAST,    SEGID,   DEGREE,
     .                     N,       STATES,  EPOCH1,  STEP     )

C$ Abstract
C
C     Write a type 8 segment to an SPK file.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     NAIF_IDS
C     SPC
C     SPK
C     TIME
C
C$ Keywords
C
C     EPHEMERIS
C     FILES
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE 'spk08.inc'

      INTEGER               HANDLE
      INTEGER               BODY
      INTEGER               CENTER
      CHARACTER*(*)         FRAME
      DOUBLE PRECISION      FIRST
      DOUBLE PRECISION      LAST
      CHARACTER*(*)         SEGID
      INTEGER               DEGREE
      INTEGER               N
      DOUBLE PRECISION      STATES ( 6, * )
      DOUBLE PRECISION      EPOCH1
      DOUBLE PRECISION      STEP

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     MAXDEG     P   Maximum degree of interpolating polynomials.
C     TOLSCL     P   Scale factor used to compute time bound tolerance.
C     HANDLE     I   Handle of an SPK file open for writing.
C     BODY       I   NAIF code for an ephemeris object.
C     CENTER     I   NAIF code for center of motion of BODY.
C     FRAME      I   Reference frame name.
C     FIRST      I   Start time of interval covered by segment.
C     LAST       I   End time of interval covered by segment.
C     SEGID      I   Segment identifier.
C     DEGREE     I   Degree of interpolating polynomials.
C     N          I   Number of states.
C     STATES     I   Array of states.
C     EPOCH1     I   Epoch of first state in STATES array.
C     STEP       I   Time step separating epochs of states.
C     MAXDEG     P   Maximum allowed degree of interpolating polynomial.
C
C$ Detailed_Input
C
C     HANDLE         is the file handle of an SPK file that has been
C                    opened for writing.
C
C     BODY           is the NAIF integer code for an ephemeris object
C                    whose state relative to another body is described
C                    by the segment to be created.
C
C     CENTER         is the NAIF integer code for the center of motion
C                    of the object identified by BODY.
C
C     FRAME          is the NAIF name for a reference frame
C                    relative to which the state information for BODY
C                    is specified.
C
C     FIRST,
C     LAST           are, respectively, the start and stop times of
C                    the time interval over which the segment defines
C                    the state of BODY.
C
C     SEGID          is the segment identifier.  An SPK segment
C                    identifier may contain up to 40 characters.
C
C     DEGREE         is the degree of the Lagrange polynomials used to
C                    interpolate the states.  All components of the
C                    state vectors are interpolated by polynomials of
C                    fixed degree.
C
C     N              is the number of states in the input state vector
C                    array.
C
C     STATES         contains a time-ordered array of geometric states
C                    ( x, y, z, dx/dt, dy/dt, dz/dt, in kilometers and
C                    kilometers per second ) of BODY relative to CENTER,
C                    specified relative to FRAME.
C
C     EPOCH1         is the epoch corresponding to the first state in
C                    the state array.  Because extra states are needed
C                    at the beginning and end of the segment in order
C                    for the interpolation method to work, EPOCH1 will
C                    normally precede FIRST.
C
C     STEP           is the time step separating the epochs of adjacent
C                    states in the input state array.  STEP is specified
C                    in seconds.
C
C$ Detailed_Output
C
C     None.  See $Particulars for a description of the effect of this
C     routine.
C
C$ Parameters
C
C     See the include file spk08.inc for declarations of the
C     parameters described below.
C
C     MAXDEG         is the maximum degree of Lagrange polynomials that
C                    can be used to interpolate states from the segment
C                    written by this routine.
C
C     TOLSCL         is a tolerance scale factor (also called a
C                    "relative tolerance") used for time coverage
C                    bound checking. TOLSCL is unitless. TOLSCL
C                    produces a tolerance value via the formula
C
C                       TOL = TOLSCL * MAX( ABS(FIRST), ABS(LAST) )
C
C                    where FIRST and LAST are the coverage time bounds
C                    of a type 8 segment, expressed as seconds past
C                    J2000 TDB.
C
C                    The resulting parameter TOL is used as a tolerance
C                    for comparing the input segment descriptor time
C                    bounds to the first and last epoch covered by the
C                    sequence of time intervals defined by the inputs
C                    to SPKW08:
C
C                       EPOCH1
C                       STEP
C                       N
C
C$ Exceptions
C
C     If any of the following exceptions occur, this routine will return
C     without creating a new segment.
C
C     1) If FRAME is not a recognized name, the error
C        SPICE(INVALIDREFFRAME) is signaled.
C
C     2) If the last non-blank character of SEGID occurs past index 40,
C        the error SPICE(SEGIDTOOLONG) is signaled.
C
C     3) If SEGID contains any nonprintable characters, the error
C        SPICE(NONPRINTABLECHARS) is signaled.
C
C     4) If DEGREE is not at least 1 or is greater than MAXDEG, the
C        error SPICE(INVALIDDEGREE) is signaled.
C
C     5) If the number of states N is not at least DEGREE+1, the error
C        SPICE(TOOFEWSTATES) will be signaled.
C
C     6) If FIRST is greater than LAST then the error
C        SPICE(BADDESCRTIMES) will be signaled.
C
C     7) If STEP is non-positive, the error SPICE(INVALIDSTEPSIZE) will
C        be signaled.
C
C     8) If the start time of the first record exceeds the descriptor
C        begin time by more than a computed tolerance, or if the end
C        time of the last record precedes the descriptor end time by
C        more than a computed tolerance, the error SPICE(COVERAGEGAP)
C        is signaled. See the Parameters section above for a
C        description of the tolerance.
C
C$ Files
C
C     A new type 8 SPK segment is written to the SPK file attached
C     to HANDLE.
C
C$ Particulars
C
C     This routine writes an SPK type 08 data segment to the open SPK
C     file according to the format described in the type 08 section of
C     the SPK Required Reading. The SPK file must have been opened with
C     write access.
C
C$ Examples
C
C     Suppose that you have states and are prepared to produce
C     a segment of type 08 in an SPK file.
C
C     The following code fragment could be used to add the new segment
C     to a previously opened SPK file attached to HANDLE. The file must
C     have been opened with write access.
C
C            C
C            C     Create a segment identifier.
C            C
C                  SEGID = 'MY_SAMPLE_SPK_TYPE_8_SEGMENT'
C
C            C
C            C     Write the segment.
C            C
C                  CALL SPKW08 (  HANDLE,  BODY,    CENTER,  FRAME,
C                .                FIRST,   LAST,    SEGID,   DEGREE,
C                .                N,       STATES,  EPOCH1,  STEP     )
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     K.R. Gehringer (JPL)
C     N.J. Bachman   (JPL)
C     J.M. Lynch     (JPL)
C     W.L. Taber     (JPL)
C
C$ Version
C
C-    SPICELIB Version 3.0.0, 11-JAN-2013 (NJB)
C
C        Relaxed test on relationship between the time bounds of the
C        input record set (determined by EPOCH1, STEP, and N) and the
C        descriptor bounds FIRST and LAST. Now the descriptor bounds
C        may extend beyond the time bounds of the record set by a ratio
C        computed using the parameter TOLSCL (see Parameters above for
C        details). MAXDEG was increased to 27.
C
C-    SPICELIB Version 2.0.0, 19-SEP-1995 (WLT)
C
C        The routine was upgraded to support non-inertial reference
C        frames.
C
C-    SPICELIB Version 1.0.1, 05-OCT-1993 (KRG)
C
C        Removed all references to a specific method of opening the SPK
C        file in the $ Brief_I/O, $ Detailed_Input, $ Particulars and
C        $ Examples sections of the header. It is assumed that a person
C        using this routine has some knowledge of the DAF system and the
C        methods for obtaining file handles.
C
C-    SPICELIB Version 1.0.0, 08-AUG-1993 (NJB) (JML) (WLT)
C
C-&

C$ Index_Entries
C
C     write spk type_8 ephemeris data segment
C
C-&

C$ Revisions
C
C-    SPICELIB Version 2.0.0, 19-SEP-1995 (WLT)
C
C        The routine was upgraded to support non-inertial reference
C        frames.
C
C-    SPICELIB Version 1.0.1, 05-OCT-1993 (KRG)
C
C        Removed all references to a specific method of opening the SPK
C        file in the $ Brief_I/O, $ Detailed_Input, $ Particulars and
C        $ Examples sections of the header. It is assumed that a person
C        using this routine has some knowledge of the DAF system and the
C        methods for obtaining file handles.
C
C-    SPICELIB Version 1.0.0, 08-AUG-1993 (NJB) (JML) (WLT)
C
C-&

C
C     SPICELIB functions
C
      INTEGER               LASTNB

      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C

C
C     SIDLEN is the maximum number of characters allowed in an
C     SPK segment identifier.
C
C     NS is the size of a packed SPK segment descriptor.
C
C     ND is the number of double precision components in an SPK
C     segment descriptor.
C
C     NI is the number of integer components in an SPK segment
C     descriptor.
C
C     DTYPE is the data type.
C
C     FPRINT is the integer value of the first printable ASCII
C     character.
C
C     LPRINT is the integer value of the last printable ASCII character.
C
C
      INTEGER               SIDLEN
      PARAMETER           ( SIDLEN =  40 )

      INTEGER               NS
      PARAMETER           ( NS     =   5 )

      INTEGER               ND
      PARAMETER           ( ND     =   2 )

      INTEGER               NI
      PARAMETER           ( NI     =   6 )

      INTEGER               DTYPE
      PARAMETER           ( DTYPE  =   8 )

      INTEGER               FPRINT
      PARAMETER           ( FPRINT =  32 )

      INTEGER               LPRINT
      PARAMETER           ( LPRINT = 126 )

      INTEGER               STATSZ
      PARAMETER           ( STATSZ =   6 )

      INTEGER               TIMLEN
      PARAMETER           ( TIMLEN =  40 )

C
C     Local variables
C
      CHARACTER*(TIMLEN)    ETSTR

      DOUBLE PRECISION      DCD    ( ND )
      DOUBLE PRECISION      DESCR  ( NS )
      DOUBLE PRECISION      LTIME
      DOUBLE PRECISION      TOL

      INTEGER               CHRCOD
      INTEGER               I
      INTEGER               ICD    ( NI )
      INTEGER               REFCOD

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'SPKW08' )

C
C     Get the NAIF integer code for the reference frame.
C
      CALL NAMFRM ( FRAME, REFCOD )

      IF ( REFCOD .EQ. 0 ) THEN

         CALL SETMSG ( 'The reference frame # is not supported.'   )
         CALL ERRCH  ( '#', FRAME                                  )
         CALL SIGERR ( 'SPICE(INVALIDREFFRAME)'                    )
         CALL CHKOUT ( 'SPKW08'                                    )
         RETURN

      END IF

C
C     Check to see if the segment identifier is too long.
C
      IF ( LASTNB(SEGID) .GT. SIDLEN ) THEN

         CALL SETMSG ( 'Segment identifier contains more than ' //
     .                 '40 characters.'                          )
         CALL SIGERR ( 'SPICE(SEGIDTOOLONG)'                     )
         CALL CHKOUT ( 'SPKW08'                                  )
         RETURN

      END IF

C
C     Now check that all the characters in the segment identifier
C     can be printed.
C
      DO I = 1, LASTNB(SEGID)

         CHRCOD = ICHAR( SEGID(I:I) )

         IF ( ( CHRCOD .LT. FPRINT ) .OR. ( CHRCOD .GT. LPRINT ) ) THEN

            CALL SETMSG ( 'The segment identifier contains '
     .               //   'nonprintable characters: ICHAR(SEGID(#:#)) '
     .               //   ' = #'                                  )
            CALL ERRINT ( '#',   I                                )
            CALL ERRINT ( '#',   I                                )
            CALL ERRINT ( '#',   CHRCOD                           )
            CALL SIGERR ( 'SPICE(NONPRINTABLECHARS)'              )
            CALL CHKOUT ( 'SPKW08'                                )
            RETURN

         END IF

      END DO

C
C     Make sure that the degree of the interpolating polynomials is
C     in range.
C
      IF (  ( DEGREE .LT. 1 ) .OR. ( DEGREE .GT. MAXDEG )  )  THEN

         CALL SETMSG ( 'The interpolating polynomials have degree #; '//
     .                 'the valid degree range is [1, #].'             )
         CALL ERRINT ( '#', DEGREE                                     )
         CALL ERRINT ( '#', MAXDEG                                     )
         CALL SIGERR ( 'SPICE(INVALIDDEGREE)'                          )
         CALL CHKOUT ( 'SPKW08'                                        )
         RETURN

      END IF

C
C     Make sure that the number of states is sufficient to define a
C     polynomial whose degree is DEGREE.
C
      IF ( N .LE. DEGREE ) THEN

         CALL SETMSG ( 'At least # states are required to define a '  //
     .                 'polynomial of degree #.  Number of states '   //
     .                 'supplied:  #.'                                 )
         CALL ERRINT ( '#', DEGREE + 1                                 )
         CALL ERRINT ( '#', DEGREE                                     )
         CALL ERRINT ( '#', N                                          )
         CALL SIGERR ( 'SPICE(TOOFEWSTATES)'                           )
         CALL CHKOUT ( 'SPKW08'                                        )
         RETURN

      END IF

C
C     The segment stop time should be greater than the begin time.
C
      IF ( FIRST .GE. LAST ) THEN

         CALL SETMSG ( 'The segment start time: # is greater than ' //
     .                 'or equal to the segment end time: #'       )
         CALL ERRDP  ( '#', FIRST                                  )
         CALL ERRDP  ( '#', LAST                                   )
         CALL SIGERR ( 'SPICE(BADDESCRTIMES)'                      )
         CALL CHKOUT ( 'SPKW08'                                    )
         RETURN

      END IF

C
C     The step size must be positive.
C
      IF ( STEP .LE. 0.D0 ) THEN

         CALL SETMSG ( 'The step size must be > 0 but was #. ' )
         CALL ERRDP  ( '#', STEP                               )
         CALL SIGERR ( 'SPICE(INVALIDSTEPSIZE)'                )
         CALL CHKOUT ( 'SPKW08'                                )
         RETURN

      END IF

C
C     Compute the tolerance to use for descriptor time bound checks.
C
      TOL = TOLSCL * MAX( ABS(FIRST), ABS(LAST) )

      IF ( FIRST .LT. EPOCH1 - TOL ) THEN

         CALL SETMSG ( 'The segment descriptor start time # is too '
     .   //            'much less than the beginning time of the  '
     .   //            'segment data # (in seconds past J2000: #). '
     .   //            'The difference is # seconds; the  '
     .   //            'tolerance is # seconds.'                     )
         CALL ETCAL  ( FIRST,  ETSTR                                 )
         CALL ERRCH  ( '#',    ETSTR                                 )
         CALL ETCAL  ( EPOCH1, ETSTR                                 )
         CALL ERRCH  ( '#',    ETSTR                                 )
         CALL ERRDP  ( '#',    FIRST                                 )
         CALL ERRDP  ( '#',    EPOCH1-FIRST                          )
         CALL ERRDP  ( '#',    TOL                                   )
         CALL SIGERR ( 'SPICE(COVERAGEGAP)'                          )
         CALL CHKOUT ( 'SPKW08'                                      )
         RETURN

      END IF

C
C     The end time of the final record must be greater than or
C     equal to the end time of the segment.
C
      LTIME = EPOCH1  +  ( (N-1) * STEP )

      IF ( LAST .GT. LTIME + TOL ) THEN

         CALL SETMSG ( 'The segment descriptor end time # is too '
     .   //            'much greater than the end time of the segment '
     .   //            'data # (in seconds past J2000: #). The '
     .   //            'difference is # seconds; the tolerance is # '
     .   //            'seconds.'                                    )
         CALL ETCAL  ( LAST,  ETSTR                                  )
         CALL ERRCH  ( '#',   ETSTR                                  )
         CALL ETCAL  ( LTIME, ETSTR                                  )
         CALL ERRCH  ( '#',   ETSTR                                  )
         CALL ERRDP  ( '#',   LAST                                   )
         CALL ERRDP  ( '#',   LAST-LTIME                             )
         CALL ERRDP  ( '#',   TOL                                    )
         CALL SIGERR ( 'SPICE(COVERAGEGAP)'                          )
         CALL CHKOUT ( 'SPKW08'                                      )
         RETURN

      END IF

C
C     If we made it this far, we're ready to start writing the segment.
C
C     Store the start and end times to be associated
C     with this segment.
C
      DCD(1) = FIRST
      DCD(2) = LAST

C
C     Create the integer portion of the descriptor.
C
      ICD(1) = BODY
      ICD(2) = CENTER
      ICD(3) = REFCOD
      ICD(4) = DTYPE

C
C     Pack the segment descriptor.
C
      CALL DAFPS ( ND, NI, DCD, ICD, DESCR )

C
C     Begin a new segment.
C
      CALL DAFBNA ( HANDLE, DESCR, SEGID )

      IF ( FAILED() ) THEN
         CALL CHKOUT ( 'SPKW08' )
         RETURN
      END IF

C
C     The type 8 segment structure is eloquently described by this
C     diagram from the SPK Required Reading:
C
C        +-----------------------+
C        | State 1               |
C        +-----------------------+
C        | State 2               |
C        +-----------------------+
C                    .
C                    .
C                    .
C        +-----------------------+
C        | State N               |
C        +-----------------------+
C        | Epoch of state 1 (ET) |
C        +-----------------------+
C        | Step size             |
C        +-----------------------+
C        | Polynomial degree     |
C        +-----------------------+
C        | Number of states      |
C        +-----------------------+
C
C
      CALL DAFADA ( STATES,          N*STATSZ )
      CALL DAFADA ( EPOCH1,          1        )
      CALL DAFADA ( STEP,            1        )
      CALL DAFADA ( DBLE(DEGREE),    1        )
      CALL DAFADA ( DBLE(N),         1        )

C
C     As long as nothing went wrong, end the segment.
C
      IF ( .NOT. FAILED() ) THEN
         CALL DAFENA
      END IF

      CALL CHKOUT ( 'SPKW08' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      PROGRAM formel05

      integer un
      parameter (un=1)

      j = un
      I = C(un)
c      I = C(1)
      call foo(C)

      end


text:
      subroutine jww(p7,p3,p4,p1,za,zb,zab,jw)
      implicit none
      include 'constants.f'
      include 'cmplxmass.f'
      include 'masses.f'
      include 'zprods_decl.f'
      include 'sprods_com.f'
      include 'ewcharge.f'
      include 'zcouple.f'
      integer jdu1,p1,p2,p3,p4,p7,ro
      double complex zab(mxpart,4,mxpart),zab2,zba2,jw(4,2,2),
     & rxw,propw34,propz17,gamv(2,2),gamvl3(2,2),gamvl4(2,2)
      double precision qn,p17(4),p34(4),t3,s134,s347,
     & s34,s137,s147,s17,q3,q4,l3,l4
C     amp(jdu1,jdu2,h17,h28)
      parameter(qn=0d0)
C-----Begin statement functions
      zab2(p1,p2,p3,p4)=za(p1,p2)*zb(p2,p4)+za(p1,p3)*zb(p3,p4)
      zba2(p1,p2,p3,p4)=zb(p1,p2)*za(p2,p4)+zb(p1,p3)*za(p3,p4)
      t3(p1,p2,p3)=s(p1,p2)+s(p2,p3)+s(p3,p1)
C-----end statement functions

      rxw=sqrt((cone-cxw)/cxw)
      s17=s(p1,p7)
      s34=s(p3,p4)
      s134=t3(p1,p3,p4)
      s137=t3(p1,p3,p7)
      s147=t3(p1,p4,p7)
      s347=t3(p3,p4,p7)
      p17(:)=0.5d0*dble(zab(p1,:,p1)+zab(p7,:,p7))
      p34(:)=0.5d0*dble(zab(p3,:,p3)+zab(p4,:,p4))

      propz17=s17-dcmplx(zmass**2,-zmass*zwidth)
      propw34=s34-dcmplx(wmass**2,-wmass*wwidth)

      q3=qn
      l3=ln
      q4=qe
      l4=le
      do jdu1=1,2
      gamV(jdu1,1)=Q(jdu1)/s17+L(jdu1)*rxw/propz17
      gamV(jdu1,2)=Q(jdu1)/s17+R(jdu1)*rxw/propz17
      gamvl3(jdu1,1)=Q(jdu1)*q3/s17+L(jdu1)*l3/propz17
      gamvl3(jdu1,2)=Q(jdu1)*q3/s17+R(jdu1)*l3/propz17
      gamvl4(jdu1,1)=Q(jdu1)*q4/s17+L(jdu1)*l4/propz17
      gamvl4(jdu1,2)=Q(jdu1)*q4/s17+R(jdu1)*l4/propz17
      enddo

      do ro=1,4
      jw(ro,1,1)= + propw34**(-1) * ( za(p7,p3)*zb(p7,p4)*zab(p7,ro,p1)
     &    *cxw**(-1)*s347**(-1) + za(p7,p3)*zb(p3,p4)*zab(p3,ro,p1)*
     &    cxw**(-1)*s347**(-1) )
      jw(ro,1,1) = jw(ro,1,1) + gamv(1,1)*propw34**(-1) * ( 2.D0*za(p7,
     &    p3)*zb(p1,p4)*p17(ro) - 2.D0*za(p7,p3)*zb(p1,p4)*p34(ro) - 2.D
     &    0*zab2(p7,p3,p4,p1)*zab(p3,ro,p4) + 2.D0*zab2(p3,p1,p7,p4)*
     &    zab(p7,ro,p1) )
      jw(ro,1,1) = jw(ro,1,1) - 2.D0*gamvl3(1,1)*za(p7,p3)*zb(p7,p1)*
     & zab(p7,ro,p4)*s137**(-1) + 2.D0*gamvl3(1,1)*za(p7,p3)*zb(p1,p3)*
     &    zab(p3,ro,p4)*s137**(-1) - 2.D0*gamvl4(1,1)*za(p7,p1)*zb(p1,
     &    p4)*zab(p3,ro,p1)*s147**(-1) - 2.D0*gamvl4(1,1)*za(p7,p4)*zb(
     &    p1,p4)*zab(p3,ro,p4)*s147**(-1)
      jw(ro,1,2)= + gamv(1,2)*propw34**(-1) * ( 2.D0*za(p1,p3)*zb(p7,p4
     &    )*p17(ro) - 2.D0*za(p1,p3)*zb(p7,p4)*p34(ro) + 2.D0*zab2(p3,
     &    p1,p7,p4)*zab(p1,ro,p7) - 2.D0*zba2(p7,p3,p4,p1)*zab(p3,ro,p4
     &    ) )
      jw(ro,1,2) = jw(ro,1,2) + 2.D0*gamvl3(1,2)*za(p1,p3)*zb(p7,p1)*
     & zab(p1,ro,p4)*s137**(-1) + 2.D0*gamvl3(1,2)*za(p1,p3)*zb(p7,p3)*
     &    zab(p3,ro,p4)*s137**(-1) + 2.D0*gamvl4(1,2)*za(p7,p1)*zb(p7,
     &    p4)*zab(p3,ro,p7)*s147**(-1) - 2.D0*gamvl4(1,2)*za(p1,p4)*zb(
     &    p7,p4)*zab(p3,ro,p4)*s147**(-1)
      jw(ro,2,1)= + propw34**(-1) * (  - za(p1,p3)*zb(p1,p4)*zab(p7,ro,
     &    p1)*cxw**(-1)*s134**(-1) + za(p3,p4)*zb(p1,p4)*zab(p7,ro,p4)*
     &    cxw**(-1)*s134**(-1) )
      jw(ro,2,1) = jw(ro,2,1) + gamv(2,1)*propw34**(-1) * ( 2.D0*za(p7,
     &    p3)*zb(p1,p4)*p17(ro) - 2.D0*za(p7,p3)*zb(p1,p4)*p34(ro) - 2.D
     &    0*zab2(p7,p3,p4,p1)*zab(p3,ro,p4) + 2.D0*zab2(p3,p1,p7,p4)*
     &    zab(p7,ro,p1) )
      jw(ro,2,1) = jw(ro,2,1) - 2.D0*gamvl3(2,1)*za(p7,p3)*zb(p7,p1)*
     & zab(p7,ro,p4)*s137**(-1) + 2.D0*gamvl3(2,1)*za(p7,p3)*zb(p1,p3)*
     &    zab(p3,ro,p4)*s137**(-1) - 2.D0*gamvl4(2,1)*za(p7,p1)*zb(p1,
     &    p4)*zab(p3,ro,p1)*s147**(-1) - 2.D0*gamvl4(2,1)*za(p7,p4)*zb(
     &    p1,p4)*zab(p3,ro,p4)*s147**(-1)
      jw(ro,2,2)= + gamv(2,2)*propw34**(-1) * ( 2.D0*za(p1,p3)*zb(p7,p4
     &    )*p17(ro) - 2.D0*za(p1,p3)*zb(p7,p4)*p34(ro) + 2.D0*zab2(p3,
     &    p1,p7,p4)*zab(p1,ro,p7) - 2.D0*zba2(p7,p3,p4,p1)*zab(p3,ro,p4
     &    ) )
      jw(ro,2,2) = jw(ro,2,2) + 2.D0*gamvl3(2,2)*za(p1,p3)*zb(p7,p1)*
     & zab(p1,ro,p4)*s137**(-1) + 2.D0*gamvl3(2,2)*za(p1,p3)*zb(p7,p3)*
     &    zab(p3,ro,p4)*s137**(-1) + 2.D0*gamvl4(2,2)*za(p7,p1)*zb(p7,
     &    p4)*zab(p3,ro,p7)*s147**(-1) - 2.D0*gamvl4(2,2)*za(p1,p4)*zb(
     &    p7,p4)*zab(p3,ro,p4)*s147**(-1)
      enddo
      jw(:,:,:)=jw(:,:,:)/sqrt(2d0*cxw)
      return
      end


! RUN: %S/test_errors.sh %s %t %flang -fopenmp
! REQUIRES: shell
! OpenMP Version 4.5
! Various checks with the ordered construct

SUBROUTINE WORK(I)
  INTEGER I
END SUBROUTINE WORK

SUBROUTINE ORDERED_GOOD(N)
  INTEGER N, I, A(10), B(10), C(10)
  !$OMP SIMD
  DO I = 1,N
    IF (I <= 10) THEN
      !$OMP ORDERED SIMD
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  END DO
  !$OMP END SIMD
END SUBROUTINE ORDERED_GOOD

SUBROUTINE ORDERED_BAD(N)
  INTEGER N, I, A(10), B(10), C(10)

  !$OMP DO SIMD
  DO I = 1,N
    IF (I <= 10) THEN
      !ERROR: The only OpenMP constructs that can be encountered during execution of a 'SIMD' region are the `ATOMIC` construct, the `LOOP` construct, the `SIMD` construct and the `ORDERED` construct with the `SIMD` clause.
      !ERROR: The ORDERED clause must be present on the loop construct if any ORDERED region ever binds to a loop region arising from the loop construct.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  END DO
  !$OMP END DO SIMD

  !$OMP PARALLEL DO
  DO I = 1,N
    IF (I <= 10) THEN
      !ERROR: The ORDERED clause must be present on the loop construct if any ORDERED region ever binds to a loop region arising from the loop construct.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  END DO
  !$OMP END PARALLEL DO

  !$OMP CRITICAL
  DO I = 1,N
    IF (I <= 10) THEN
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  END DO
  !$OMP END CRITICAL

  !$OMP CRITICAL
    WRITE(*,*) I
    !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
    !$OMP ORDERED
    CALL WORK(I)
    !$OMP END ORDERED
  !$OMP END CRITICAL

  !$OMP ORDERED
    WRITE(*,*) I
    IF (I <= 10) THEN
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  !$OMP END ORDERED

  !$OMP TASK
    C =  C - A * B
    !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
    !$OMP ORDERED
    CALL WORK(I)
    !$OMP END ORDERED
  !$OMP END TASK

  !$OMP TASKLOOP
  DO I = 1,N
    IF (I <= 10) THEN
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    ENDIF
  END DO
  !$OMP END TASKLOOP

  !$OMP CRITICAL
    C =  C - A * B
    !$OMP MASTER
    DO I = 1,N
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    END DO
    !$OMP END MASTER
  !$OMP END CRITICAL

  !$OMP ORDERED
    C =  C - A * B
    !$OMP MASTER
    DO I = 1,N
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    END DO
    !$OMP END MASTER
  !$OMP END ORDERED

  !$OMP TASK
    C =  C - A * B
    !ERROR: `MASTER` region may not be closely nested inside of `WORKSHARING`, `LOOP`, `TASK`, `TASKLOOP`, or `ATOMIC` region.
    !$OMP MASTER
    DO I = 1,N
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    END DO
    !$OMP END MASTER
  !$OMP END TASK

  !$OMP TASKLOOP
  DO J= 1,N
    C =  C - A * B
    !ERROR: `MASTER` region may not be closely nested inside of `WORKSHARING`, `LOOP`, `TASK`, `TASKLOOP`, or `ATOMIC` region.
    !$OMP MASTER
    DO I = 1,N
      !ERROR: `ORDERED` region may not be closely nested inside of `CRITICAL`, `ORDERED`, explicit `TASK` or `TASKLOOP` region.
      !$OMP ORDERED
      CALL WORK(I)
      !$OMP END ORDERED
    END DO
    !$OMP END MASTER
  END DO
  !$OMP END TASKLOOP

END SUBROUTINE ORDERED_BAD


text:
module module_out

use module_types
use module_status
use module_utils
use module_datetime
use module_options
use netcdf
use module_float
use module_forcing, only: NLAYER,LAYER

implicit none

character(len=maxlen)                    :: output_mode = ""
character(len=maxlen)                    :: output_name = ""
integer                                  :: output_id   = 34
integer                                  :: output_record = 0
integer                                  :: output_nid
!integer                                  :: output_strid
integer                                  :: output_idl
integer                                  :: output_timeid
!integer                                  :: output_dateid
integer                                  :: output_lonid
integer                                  :: output_latid
integer                                  :: output_zid
integer                                  :: output_uid
integer                                  :: output_vid
integer                                  :: output_rid
integer                                  :: output_nfloats
integer                                  :: output_xoid
integer                                  :: output_yoid
integer                                  :: output_zoid
integer                                  :: output_toid
integer                                  :: output_status
integer                                  :: output_kid
real(dp)                                 :: output_missing

integer                                  :: input_id
integer                                  :: input_timeid

character(180)                           :: out_units
character(180)                           :: model_units
character(20)                            :: model_calendar

contains
! ...
! ====================================================================
! ...
subroutine trajectory_open(filename,Nfloats,missing)

character(len=*), intent(in)               :: filename
integer, intent(in)                        :: Nfloats
real(dp), intent(in)                       :: missing

! ... Local variables
! ...
integer i,natts
character(len=maxlen) ext
character(len=1000) lcom

STATUS_ERROR = 0; STATUS_TEXT = ""

output_nfloats = Nfloats
output_missing = missing
output_record = 0

! ... Check for extension
! ...
do i=len(filename),1,-1
  if (filename(i:i).eq.'.') then
    ext = filename(i+1:)
    exit
  endif
enddo

if (lowercase(ext).eq.'nc') then
  output_mode = 'nc'
  output_name = trim(filename)

  STATUS_ERROR = NF90_CREATE(filename,NF90_CLOBBER,output_id); call cdf_error()

  STATUS_ERROR = NF90_DEF_DIM(output_id,'floats',Nfloats,output_nid);      call cdf_error()
!  STATUS_ERROR = NF90_DEF_DIM(output_id,'string19',19,output_strid);       call cdf_error()
  STATUS_ERROR = NF90_DEF_DIM(output_id,'time',NF90_UNLIMITED,output_idl); call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'release_longitude',NF90_DOUBLE,(/output_nid/),output_xoid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_xoid,'units','degrees_east')
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'release_latitude',NF90_DOUBLE,(/output_nid/),output_yoid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_yoid,'units','degrees_north')
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'release_depth',NF90_DOUBLE,(/output_nid/),output_zoid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_zoid,'units','meters')
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'model_level',NF90_INT,(/output_nid/),output_kid)
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'release_time',NF90_DOUBLE,(/output_nid/),output_toid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_toid,'units','seconds since initial time')
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'time',NF90_DOUBLE,(/output_idl/),output_timeid)
  call cdf_error()
  call cdf_copyatts(.false.,input_id,input_timeid,output_id,output_timeid,natts)
  !STATUS_ERROR = NF90_PUT_ATT(output_id,output_timeid,'units','Julian days')
  !call cdf_error()

!  STATUS_ERROR = NF90_DEF_VAR(output_id,'date',NF90_CHAR,(/output_strid,output_idl/),output_dateid)
!  call cdf_error()
!  STATUS_ERROR = NF90_PUT_ATT(output_id,output_dateid,'units','ISO 8601 date and time')
!  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'lon',NF90_REAL,(/output_nid,output_idl/),output_lonid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_lonid,'units','degrees_east')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_lonid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'lat',NF90_REAL,(/output_nid,output_idl/),output_latid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_latid,'units','degrees_north')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_latid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'depth',NF90_REAL,(/output_nid,output_idl/),output_zid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_zid,'units','meters')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_zid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'u',NF90_REAL,(/output_nid,output_idl/),output_uid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_uid,'units','m/s')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_uid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'v',NF90_REAL,(/output_nid,output_idl/),output_vid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_vid,'units','m/s')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_vid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'distance',NF90_REAL,(/output_nid,output_idl/),output_rid)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_rid,'units','m')
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_rid,'_FillValue',sngl(missing))
  call cdf_error()

  STATUS_ERROR = NF90_DEF_VAR(output_id,'exitcode',NF90_INT,(/output_nid/),output_status)
  call cdf_error()
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_status,'long_name','Float status at the end of the simulation')
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_status,"is-1",'float has not been released')
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_status,"is0",'float was moving')
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_status,"is1",'float left the model area')
  STATUS_ERROR = NF90_PUT_ATT(output_id,output_status,"is2",'float was stranded')
  call cdf_error()

  call get_commandline(lcom)
  STATUS_ERROR = NF90_PUT_ATT(output_id,0,'CommandLine',TRIM(lcom))
  call cdf_error()

  STATUS_ERROR = NF90_ENDDEF(output_id)
  call cdf_error()

endif

end subroutine trajectory_open
! ...
! ====================================================================
! ...
subroutine param_write(p1,p2,p3,p4,p5,p6,p7,p8)

real(dp), intent(in)         :: p1,p2,p3,p4,p5,p6,p7,p8

STATUS_ERROR = NF90_REDEF(output_id); call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'alpha',p1);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'A11',p2);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'A12',p3);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'A21',p4);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'A22',p5);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'mu',p6);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'K1',p7);   call cdf_error()
STATUS_ERROR = NF90_PUT_ATT(output_id,0,'K0',p8);   call cdf_error()
STATUS_ERROR = NF90_ENDDEF(output_id); call cdf_error()

end subroutine param_write
! ...
! ====================================================================
! ...
subroutine release_write()

STATUS_ERROR = NF90_PUT_VAR(output_id,output_xoid,FLT%release_lon);   call cdf_error()
STATUS_ERROR = NF90_PUT_VAR(output_id,output_yoid,FLT%release_lat);   call cdf_error()
STATUS_ERROR = NF90_PUT_VAR(output_id,output_zoid,FLT%release_depth); call cdf_error()
STATUS_ERROR = NF90_PUT_VAR(output_id,output_toid,FLT%release_time);  call cdf_error()
STATUS_ERROR = NF90_PUT_VAR(output_id,output_kid,LAYER(FLT%k));  call cdf_error()

end subroutine release_write
! ...
! ====================================================================
! ...
subroutine exitcode_write()

! ... Local variables:
! ...
integer i, code(FLT%Nfloats)

do i=1,FLT%Nfloats
  if (.not.FLT%released(i)) then
    code(i) = -1
  else if (FLT%outside(i)) then
    code(i) = 1
  else if (FLT%stranded(i)) then
    code(i) = 2
  else
    code(i) = 0
  endif
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_status,code);   call cdf_error()

end subroutine exitcode_write
! ...
! ====================================================================
! ...
subroutine trajectory_write(time)

real(dp), intent(in)                                :: time

! ... Local variables
! ...
integer i
real(dp) time_out(1), xx(output_nfloats)

STATUS_ERROR = 0; STATUS_TEXT = ""

output_record = output_record + 1

time_out = time_transform([time],model_units,model_calendar,out_units,model_calendar)
STATUS_ERROR = NF90_PUT_VAR(output_id,output_timeid,time_out,[output_record],[1])

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%lon(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_lonid,xx,[1,output_record],[output_nfloats,1])

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%lat(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_latid,xx,[1,output_record],[output_nfloats,1])
call cdf_error()

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%z(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_zid,xx,[1,output_record],[output_nfloats,1])
call cdf_error()

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%dist(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_rid,xx,[1,output_record],[output_nfloats,1])
call cdf_error()

end subroutine trajectory_write
! ....
! ====================================================================
! ...
subroutine velocity_write()

! ... Local variables
! ...
integer i
real(dp) xx(output_nfloats)

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%u(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_uid,xx,[1,output_record],[output_nfloats,1])
call cdf_error()

xx(:) = output_missing
do i=1,FLT%Nfloats
  if (FLT%released(i)) xx(i) = FLT%v(i)
enddo
STATUS_ERROR = NF90_PUT_VAR(output_id,output_vid,xx,[1,output_record],[output_nfloats,1])
call cdf_error()

end subroutine velocity_write
! ....
! ====================================================================
! ...
subroutine trajectory_close()

! ... Local variables
! ...

STATUS_ERROR = 0; STATUS_TEXT = ""

if (trim(output_mode).eq.'text') then
  close(output_id)
else if (trim(output_mode).eq.'json') then
  close(output_id)
else if (trim(output_mode).eq.'nc') then
  STATUS_ERROR = NF90_CLOSE(output_id)
  if (STATUS_ERROR.ne.NF90_NOERR) then
    STATUS_TEXT = trim(NF90_STRERROR(STATUS_ERROR))
    return
  endif
endif

return

end subroutine trajectory_close
! ....
! ====================================================================
! ...
!subroutine cdf_error()
!
!if (STATUS_ERROR.ne.NF90_NOERR) then
!  STATUS_TEXT = trim(NF90_STRERROR(STATUS_ERROR))
!  return
!endif
!
!end subroutine cdf_error
! ....
! ====================================================================
! ...

end module module_out

max_stars_repo_name: quimbp/cosmo

max_issues_repo_name: quimbp/cosmo

max_forks_repo_name: quimbp/cosmo


text:
!! R440 proc-component-def-stmt
!    is PROCEDURE ( [ proc-interface ] ) ,
!            proc-component-attr-spec-list :: proc-decl-list
!
module whut

  type point
     integer :: x,y
  end type point

  type, extends(point) :: feeble
     procedure (), pointer, nopass :: foible
     procedure (bob), pass(a), private, pointer :: foo => bar
     procedure (bob), pass, public, pointer :: foo => bar
end type

interface
   function bob(a)
     import feeble
     type(feeble) :: a
   end function bob
end interface

contains

  function boo(a)
    type(point) :: a
  end function boo

end module whut

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf