text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130
Subroutine init_Ac
  use Global_Variables
  use salmon_parallel, only: nproc_group_global, nproc_id_global
  use salmon_communication, only: comm_bcast, comm_is_root
  use Ac_alocal_laser
  implicit none
  integer :: iter, npower
  real(8) :: tt

  javt = 0d0
  Ac_ext = 0d0
  Ac_ind = 0d0

  if(rlaser_int_wcm2_1 < 0d0)then
    f0_1 = amplitude1
  else
    f0_1=5.338d-9*sqrt(rlaser_int_wcm2_1)      ! electric field in a.u.
  end if
  if(rlaser_int_wcm2_2 < 0d0)then
    f0_2 = amplitude2
  else
    f0_2=5.338d-9*sqrt(rlaser_int_wcm2_2)      ! electric field in a.u.
  end if

!  f0_1=5.338d-9*sqrt(IWcm2_1)      ! electric field in a.u.
!  f0_2=5.338d-9*sqrt(IWcm2_2)      ! electric field in a.u.
!  omega_1=omegaev_1/(2d0*Ry)  ! frequency in a.u.
!  tpulse_1=tpulsefs_1/0.02418d0 ! pulse duration in a.u.

!  omega_2=omegaev_2/(2d0*Ry)  ! frequency in a.u.
!  tpulse_2=tpulsefs_2/0.02418d0 ! pulse duration in a.u.
!  T1_T2=T1_T2fs/0.02418d0 ! pulse duration in a.u.
  javt=0.d0
  Ac_ext=0.d0

  select case(AE_shape1)
  case('impulse')
    Ac_ext(:,1)=epdir_re1(1)*e_impulse
    Ac_ext(:,2)=epdir_re1(2)*e_impulse
    Ac_ext(:,3)=epdir_re1(3)*e_impulse
  case('Acos2','Acos3','Acos4','Acos6','Acos8')
    select case(ae_shape1)
    case('Acos2'); npower = 2
    case('Acos3'); npower = 3
    case('Acos4'); npower = 4
    case('Acos6'); npower = 6
    case('Acos8'); npower = 8
    case default
      stop 'Error in init_Ac.f90'
    end select

   !do iter=0,Nt+1
     !tt=iter*dt - 0.5d0*pulse_tw1
    do iter=-1,Nt+1
      if(iter==-1 .and. t1_delay.ge.0d0) cycle !only for restart of field using rt_wfn_k
      tt=iter*dt - 0.5d0*pulse_tw1 - t1_delay
      if (abs(tt)<0.5d0*pulse_tw1) then
        Ac_ext(iter,:)=-f0_1/omega1*(cos(pi*tt/pulse_tw1))**npower &
          *aimag( (epdir_re1(:) + zI*epdir_im1(:)) &
          *exp(zI*(omega1*tt+phi_CEP1*2d0*pi))  &
          )
      end if
    enddo
    T1_T2 = T1_T2 + t1_delay

  case('Ecos2')
    if(phi_CEP1 /= 0.75d0)then
      call Err_finalize("Error: phi_cep1 should be 0.75 when ae_shape1 is 'Ecos2'.")
    end if
    do iter=0,Nt+1
      tt=iter*dt - 0.5d0*pulse_tw1
      if (abs(tt)<0.5d0*pulse_tw1) then
        Ac_ext(iter,:)=-epdir_re1(:)*f0_1/(8d0*pi**2*omega1 - 2d0*pulse_tw1**2*omega1**3) &
          *( &
          (-4d0*pi**2+pulse_tw1**2*omega1**2 + pulse_tw1**2*omega1**2*cos(2d0*pi*tt/pulse_tw1))*cos(omega1*tt) &
          +2d0*pi*(2d0*pi*cos(pulse_tw1*omega1/2d0) &
          +pulse_tw1*omega1*sin(2d0*pi*tt/pulse_tw1)*sin(omega1*tt)))
      end if
    enddo

  case('Esin2sin')
    do iter=0,Nt+1
      tt=iter*dt
      if (tt<pulse_tw1) then
        Ac_ext(iter,:)=epdir_re1(:)*f0_1*(&
          &-(cos(omega1*tt+phi_CEP1*2d0*pi)-cos(phi_CEP1*2d0*pi))/(2*omega1)&
          &+(cos((omega1+2*Pi/pulse_tw1)*tt+phi_CEP1*2d0*pi)-cos(phi_CEP1*2d0*pi))/(4*(omega1+2*Pi/pulse_tw1))&
          &+(cos((omega1-2*Pi/pulse_tw1)*tt+phi_CEP1*2d0*pi)-cos(phi_CEP1*2d0*pi))/(4*(omega1-2*Pi/pulse_tw1))&
          &-0.5*(cos(omega1*pulse_tw1+phi_CEP1*2d0*pi)-cos(phi_CEP1*2d0*pi))&
          &/(omega1*((omega1*pulse_tw1/(2*pi))**2-1.d0))&
          &*(3.d0*(tt/pulse_tw1)**2-2.d0*(tt/pulse_tw1)**3)&
          &)
      else
        Ac_ext(iter,:)=Ac_ext(iter-1,:)
      endif
    enddo

  case('Asin2cos')
    ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! pump laser
    do iter=0,Nt+1
      tt=iter*dt
      if (tt<pulse_tw1) then
        Ac_ext(iter,:)=-epdir_re1(:)*f0_1/omega1*(sin(pi*tt/pulse_tw1))**2*cos(omega1*tt+phi_CEP1*2d0*pi)
      end if
    enddo

  case('input')
    Ac_ext=0d0
    if(comm_is_root(nproc_id_global))then
      open(899,file='input_Ac.dat')
      do iter=0,Nt+1
        read(899,*)Ac_ext(iter,1),Ac_ext(iter,2),Ac_ext(iter,3)
      end do
      close(899)
    end if
    call comm_bcast(Ac_ext,nproc_group_global)

  case('Asin2_cw')
    ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! pump laser
    do iter=0,Nt+1
      tt=iter*dt
      if (tt<pulse_tw1*0.5d0) then
        Ac_ext(iter,:)=-Epdir_re1(:)*f0_1/omega1*(sin(pi*tt/pulse_tw1))**2*cos(omega1*tt+phi_CEP1*2d0*pi)
      else
        Ac_ext(iter,:)=-Epdir_re1(:)*f0_1/omega1*cos(omega1*tt+phi_CEP1*2d0*pi)
      end if
    enddo
  case('none')
    !there is no pump
  case default
    call Err_finalize("Invalid pulse_shape_1 parameter!")
  end select

! Probe
  select case(ae_shape2)
  case('impulse')
    do iter=0,Nt+1
      tt=iter*dt
      if(tt > T1_T2)then
        Ac_ext(iter,1)=Ac_ext(iter,1) + epdir_re2(1)*e_impulse
        Ac_ext(iter,2)=Ac_ext(iter,2) + epdir_re2(2)*e_impulse
        Ac_ext(iter,3)=Ac_ext(iter,3) + epdir_re2(3)*e_impulse
      end if
    end do
  case('Acos2','Acos3','Acos4','Acos6','Acos8')
    select case(ae_shape2)
    case('Acos2'); npower = 2
    case('Acos3'); npower = 3
    case('Acos4'); npower = 4
    case('Acos6'); npower = 6
    case('Acos8'); npower = 8
    case default
      stop 'Error in init_Ac.f90'
    end select

    do iter=0,Nt+1
      tt=iter*dt - 0.5d0*pulse_tw1 - T1_T2
      if (abs(tt)<0.5d0*pulse_tw2) then
        Ac_ext(iter,:)=Ac_ext(iter,:) &
          -f0_2/omega2*(cos(pi*tt/pulse_tw2))**npower &
          *aimag( (epdir_re2(:) + zI*epdir_im2(:)) &
          *exp(zI*(omega2*tt+phi_CEP2*2d0*pi))  &
          )
      end if
    end do

  case('Ecos2')
    if(phi_CEP2 /= 0.75d0)then
      call Err_finalize("Error: phi_cep2 should be 0.75 when ae_shape2 is 'Ecos2'.")
    end if
    do iter=0,Nt+1
      tt=iter*dt - 0.5d0*pulse_tw1 - T1_T2
      if (abs(tt)<0.5d0*pulse_tw2) then
        Ac_ext(iter,:)=Ac_ext(iter,:) &
          -epdir_re2(:)*f0_2/(8d0*pi**2*omega2 - 2d0*pulse_tw2**2*omega2**3) &
          *( &
          (-4d0*pi**2+pulse_tw2**2*omega2**2 + pulse_tw2**2*omega2**2*cos(2d0*pi*tt/pulse_tw2))*cos(omega2*tt) &
          +2d0*pi*(2d0*pi*cos(pulse_tw2*omega2/2d0) &
          +pulse_tw2*omega2*sin(2d0*pi*tt/pulse_tw2)*sin(omega2*tt)))
      end if
    enddo

  case('Esin2sin')
      ! probe laser
    do iter=0,Nt+1
      tt=iter*dt
      if(tt-T1_T2 <= 0d0)then

      else if ( (tt-T1_T2>0d0) .and. (tt-T1_T2<pulse_tw2) ) then
        Ac_ext(iter,:)=Ac_ext(iter,:)+Epdir_re2(:)*f0_2*(&
          &-(cos(omega2*(tt-T1_T2)+phi_CEP2*2d0*pi)-cos(phi_CEP2*2d0*pi))/(2*omega2)&
          &+(cos((omega2+2*Pi/pulse_tw2)*(tt-T1_T2)+phi_CEP2*2d0*pi)-cos(phi_CEP2*2d0*pi))&
          &/(4*(omega2+2*Pi/pulse_tw2))&
          &+(cos((omega2-2*Pi/pulse_tw2)*(tt-T1_T2)+phi_CEP2*2d0*pi)-cos(phi_CEP2*2d0*pi))&
          /(4*(omega2-2*Pi/pulse_tw2))&
          &-0.5*(cos(omega2*pulse_tw2+phi_CEP2*2d0*pi)-cos(phi_CEP2*2d0*pi))/(omega2*((omega2*pulse_tw2/(2*pi))**2-1.d0))&
          &*(3.d0*((tt-T1_T2)/pulse_tw2)**2-2.d0*((tt-T1_T2)/pulse_tw2)**3)&
          &)
      else
        Ac_ext(iter,:)=Ac_ext(iter-1,:)
      endif
    enddo

  case('Asin2cos')
      ! pulse shape : A(t)=f0/omega*sin(Pi t/T)**2 *cos (omega t+phi_CEP*2d0*pi)
    ! probe laser
    do iter=0,Nt+1
      tt=iter*dt
      if ( (tt-T1_T2>0d0) .and. (tt-T1_T2<pulse_tw2) ) then
        Ac_ext(iter,:)=Ac_ext(iter,:) &
          &-Epdir_re2(:)*f0_2/omega2*(sin(pi*(tt-T1_T2)/pulse_tw2))**2*cos(omega2*(tt-T1_T2)+phi_CEP2*2d0*pi)
      endif
    enddo
  case('input')
    !There is no probe
  case('Asin2_cw')
    !There is no probe
  case('none')
  case default
    call Err_finalize("Invalid pulse_shape_2 parameter!")
  end select
  Ac_ind=0.d0

  Ac_tot=Ac_ind+Ac_ext

  !hidden option (AY)
  if(alocal_laser=='y') call init_Ac_alocal

  return
End Subroutine init_Ac
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
program ch0902
    ! Integer Overflow and the i Edit Descriptor
    implicit none

    integer :: big = 10
    integer :: i

    do i = 1, 40
        print 100, i, big
        big = big*10
    end do

100 format (' ', i3, ' ', i12)
end program

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran


text:
      SUBROUTINE PLT_LWF_LABEL
     &          (xlabel,ylabel,
     &           file_id,prgm_id,
     &           case_id,prfl_id,
     &           xmtr_id,freq,tlat,tlon,bearng,
     &           power,incl,headng,talt,ralt,nravg,
     &           plt_label,mxlbl,nrlbl)

c***********************************************************************
c  Subroutine       plt_lwf_label
c***********************************************************************
c
c  Program Source:
c     Naval Ocean Systems Center - Code 542
c
c  Date:
c     26 Jan 1996
c
c  Function:
c     Places labels for LWF plots on a graph and returns them to the
c     calling routine.
c
c  Parameters passed:
c     xlabel        [r] position of the first and final line of the label
c     ylabel        [r]
c
c     file_id       [s] file identification
c     prgm_id       [s] program which generated the data
c     case_id       [s] user identification of the data
c     prfl_id       [s] ionospheric profile identification
c
c     xmtr_id       [s] transmitter identification
c     freq          [r] frequency
c     tlat          [r] transmitter latitude
c     tlon          [r] transmitter longitude
c     bearng        [r] bearing of the path
c
c     power         [r] transmitter power
c     incl          [r] antenna inclination
c     headng        [r] antenna heading
c     talt          [r] transmitter altitude
c     ralt          [r] receiver    altitude
c
c     nravg         [i] number of points for running average
c
c  Parameters returned:
c     plt_label     [s] array of data labels
c     mxlbl         [i] dimension of PLT_LABEL
c     nrlbl         [i] number    of PLT_LABELs generated
c
c  Common blocks referenced:
c     None
c
c  Functions and subroutines referenced:
c     None
c
c  References:
c     None
c
c  Change History:
c
c     12 Mar 1997   Modified to use new string routines in GRF.
c
c*******************!***************************************************

      character*(*) plt_label(mxlbl)
      character*  8 prgm_id
      character* 20 xmtr_id
      character* 40 prfl_id
      character* 60 frmt1,frmt2
      character* 80 case_id
      character*120 file_id(3)
      integer       str_length
      real          freq,tlat,tlon,bearng,
     &              power,incl,headng,talt,ralt,
     &              xlabel,ylabel

      numx=STR_LENGTH(xmtr_id)
      nums=numx-7
      if (numx .eq. 0) then
         write(frmt1,'(''(a,2x,a,a)'')')
         write(frmt2,'(''(7x,2f6.1,f7.1,f6.1,i5,i3,f6.1,2f5.1)'')')
      else
     &if (nums .eq. 0) then
         write(frmt1,'(''(a,2x,a,a)'')')
         write(frmt2,'(''(a,2f6.1,f7.1,f6.1,i5,i3,f6.1,2f5.1)'')')
      else
     &if (nums .lt. 0) then
         write(frmt1,'(''(a,2x,a,a)'')')
         write(frmt2,'(''(a,'',i2.2,''x,2f6.1,'',
     &                 ''f7.1,f6.1,i5,i3,f6.1,2f5.1)'')') -nums
      else
         write(frmt1,'(''(a,'',i2.2,''x,2x,a,a)'')') nums
         write(frmt2,'(''(a,2f6.1,f7.1,f6.1,i5,i3,f6.1,2f5.1)'')')
      end if

      nrlbl=1
      plt_label(nrlbl)=case_id
      nrlbl=nrlbl+1
      write(plt_label(nrlbl),frmt1)
     &     'xmtr_id',
     &     'freq  tlat   tlon  brng  pwr ',
     &     'in   hdg talt ralt'
      nrlbl=nrlbl+1
      if (numx .eq. 0) then
         write(plt_label(nrlbl),frmt2)
     &         freq,tlat,tlon,bearng,int(power),
     &         int(incl),headng,talt,ralt
      else
         write(plt_label(nrlbl),frmt2)
     &         xmtr_id(:numx),
     &         freq,tlat,tlon,bearng,int(power),
     &         int(incl),headng,talt,ralt
      end if
      nrlbl=nrlbl+1
      write(plt_label(nrlbl),'(''Prfl:  '',a)') prfl_id
      nrlbl=nrlbl+1
      write(plt_label(nrlbl),'(''File:  '',a)') file_id(1)
      nrlbl=nrlbl+1
      write(plt_label(nrlbl),'(''File:  '',a)') file_id(2)
      nrlbl=nrlbl+1
      write(plt_label(nrlbl),'(''Prgm:  '',a)') prgm_id

c     Add labels for running average if necessary
      if (nravg .gt. 1) then
         nrlbl=nrlbl+1
         if (nravg .lt. 10) then
            write(plt_label(nrlbl),
     &          '(''nAvg:  '',i1)') nravg
         else
            write(plt_label(nrlbl),
     &          '(''nAvg:  '',i2)') nravg
         end if
      end if

      xp=xlabel
      yp=ylabel

      do n=1,nrlbl
         call GRF_STRING (xp,yp,.1,plt_label(n),0.,'LB',retX)
         yp=yp-.15
      end do

c     Return the final label position
      xlabel=xp
      ylabel=yp

      RETURN
      END      ! PLT_LWF_LABEL

ext: for

max_stars_repo_name: spinorkit/LWPC

max_issues_repo_name: spinorkit/LWPC

max_forks_repo_name: spinorkit/LWPC


text:
.ND
.ce 6
ADDENDUM NUMBER TWO
.sp 1
to
.sp 1
FOURTH BERKELEY SOFTWARE DISTRIBUTION
.sp 2
LICENSE (6/83) AGREEMENT BETWEEN
.sp 2
THE REGENTS OF THE UNIVERSITY OF CALIFORNIA AND FOREIGN LICENSEES
.sp 2
.IP 1.
The above referenced License Agreement and Addendum Number
One ("Agreement") entered into between
The Regents of the University of California (``Regents'') and the below-named
Licensee controls Licensee's use of release 4.3 of Regents' Fourth Berkeley
Software Distribution (``4BSD'') which was distributed upon execution of
the Agreement.  The Tahoe version of 4.3BSD is hereby made available to
Licensee under the same terms and conditions set forth in the Agreement.
.sp 2
.PP
IN WITNESS WHEREOF, Regents and Licensee have executed this
Addendum Number Two on the dates specified below.  The effective date of this
amendment shall be the date of execution by Regents.
.sp 2
.TS
expand tab(~);
cw(2.9375i) cw(0.75i) cw(2.9375i).
THE REGENTS OF THE UNIVERSITY~~
\v'-10p'\l'2.9375i'~~\l'2.9375i'
\v'-7p'OF CALIFORNIA~~
\v'-17p'\l'2.9375i'~~\l'2.9375i'
\v'-15p'(LICENSOR)~~(LICENSEE)
.T&
lw(2.9375i) cw(0.75i) lw(2.9375i).
~
~
~
By\l'2.75i'~~By\l'2.75i'
~
Title\l'2.625i'~~Title\l'2.625i'
~
Date\l'2.625i'~~Date\l'2.625i'
.TE

max_stars_repo_name: weiss/original-bsd

max_issues_repo_name: JamesLinus/original-bsd

max_forks_repo_name: JamesLinus/original-bsd


text:
      subroutine cmo_create_cmo(imsgin,xmsgin,cmsgin,msgtype,nwds,
     *                          ierror_return)
C
C
C#######################################################################
C
C     PURPOSE -
C
C        This routine creates a Mesh Object.
C
C     INPUT ARGUMENTS -
C
C        imsgin()  - Integer array of command input tokens
C        xmsgin()  - Real array of command input tokens
C        cmsgin()  - Character array of command input tokens
C        msgtype() - Integer array of command input token types
C        nwds      - Number of command input tokens
C
C     OUTPUT ARGUMENTS -
C
C        ierror_return - Error Return Code (==0 ==> OK, <>0 ==> Error)
C
C     CHANGE HISTORY -
C
C        $Log: cmo_create_cmo.f,v $
C        Revision 2.00  2007/11/05 19:45:48  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.1   27 Jul 2001 11:27:40   tam
CPVCS    changed call to cmo_mesh_type to cmo_set_mesh_type
CPVCS
CPVCS       Rev 1.0   24 Jan 2000 13:24:38   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.6   Mon Apr 14 16:40:10 1997   pvcs
CPVCS    No change.
CPVCS
CPVCS       Rev 1.5   03/15/95 15:22:36   ejl
CPVCS    Finished installing the defaults.
CPVCS
CPVCS       Rev 1.4   02/23/95 14:52:14   ejl
CPVCS    Fixed so defaults will work:
CPVCS    COMMAND////input//.
CPVCS
CPVCS       Rev 1.3   02/22/95 09:31:00   ejl
CPVCS
CPVCS
CPVCS       Rev 1.2   02/16/95 16:11:32   ejl
CPVCS
CPVCS
CPVCS       Rev 1.1   02/16/95 14:40:16   ejl
CPVCS    Fixed error in if statement in defaults.
CPVCS
CPVCS
CPVCS       Rev 1.0   02/16/95 13:54:44   dcg
CPVCS    Original Version
C
C#######################################################################
C
      implicit none
C
C#######################################################################
C
      integer nwds, imsgin(nwds), msgtype(nwds)
      integer ierror_return
      REAL*8 xmsgin(nwds)
      character*32 cmsgin(nwds)
C
C
C#######################################################################
C
C     LOCAL VARIABLE DEFINITION
C
C
      integer len, ierror, number_nodes, number_elements

      integer icharlnf

      character*32 cmo_name
      character*32 isubname
C
C
C#######################################################################
C BEGIN begin
C
C
C.... Mesh Object Name.
C
      isubname='cmo_create'
      cmo_name=cmsgin(3)
C
      len=icharlnf(cmo_name)
C
      if((cmo_name(1:len).eq.'-cmo-') .or.
     *   (cmo_name(1:len).eq.'-def-')) then
C
C....    Use the Current Mesh Object.
C
         call cmo_get_name(cmo_name,ierror_return)
C
      endif
C
C.... Create the Mesh Object from the Default Mesh Object.
C
      number_elements=0
      number_nodes=0
      call cmo_create(cmo_name,ierror_return)
C
      if((ierror_return.eq.0) .and. (nwds.gt.3)) then
C
         if((nwds.ge.4) .and. (cmsgin(4)(1:5).ne.'-def-') .and.
     *                        (cmsgin(4)(1:9).ne.'-default-')) then
            number_nodes=imsgin(4)
         else
            number_nodes=0
         endif
C
         if((nwds.ge.5) .and. (cmsgin(5)(1:5).ne.'-def-') .and.
     *                        (cmsgin(5)(1:9).ne.'-default-')) then
            number_elements=imsgin(5)
         else
            number_elements=0
         endif
C
         if(msgtype(6).eq.3) then
C
C....       Mesh Type Specified.
C
            call cmo_set_mesh_type(cmo_name, cmsgin(6), ierror_return)
C
         else
C
C....       Store data in the Mesh Object.
C
            if((nwds.ge.6) .and. (cmsgin(6)(1:5).ne.'-def-') .and.
     *                           (cmsgin(6)(1:9).ne.'-default-')) then
               call cmo_set_info('ndimensions_geom',cmo_name,imsgin(6),
     *                           1,1,ierror)
            endif
C
            if((nwds.ge.7) .and. (cmsgin(7)(1:5).ne.'-def-') .and.
     *                           (cmsgin(7)(1:9).ne.'-default-')) then
               call cmo_set_info('ndimensions_topo',cmo_name,imsgin(7),
     *                           1,1,ierror)
            endif
C
            if((nwds.ge.8) .and. (cmsgin(8)(1:5).ne.'-def-') .and.
     *                           (cmsgin(8)(1:9).ne.'-default-')) then
               call cmo_set_info('nodes_per_element',cmo_name,imsgin(8),
     *                           1,1,ierror)
            endif
C
            if((nwds.ge.9) .and. (cmsgin(9)(1:5).ne.'-def-') .and.
     *                           (cmsgin(9)(1:9).ne.'-default-')) then
               call cmo_set_info('faces_per_element',cmo_name,imsgin(9),
     *                           1,1,ierror)
            endif
C
            if((nwds.ge.10) .and. (cmsgin(10)(1:5).ne.'-def-') .and.
     *                            (cmsgin(10)(1:9).ne.'-default-')) then
               call cmo_set_info('edges_per_element',cmo_name,
     *                           imsgin(10),1,1,ierror)
            endif
C
         endif
C
C....    Set all Memory Managed array lengths.
C

         call cmo_memory(cmo_name,number_nodes,number_elements,
     *                   ierror_return)
C
      endif

      if (ierror_return.eq.0) then
      call cmo_set_info('nelements',cmo_name,number_elements,1,1,ierror)
      if(ierror.ne.0) call x3d_error(isubname,'cmo_set_info nelements')

      call cmo_set_info('nnodes',cmo_name,number_nodes,1,1,ierror)
      if(ierror.ne.0) call x3d_error(isubname,'cmo_set_info nodes')

      endif
C
      return
      end

ext: f

max_forks_repo_name: daniellivingston/LaGriT


text:
PROGRAM Numeric_Inquiry

  ! ------------
  ! Module usage
  ! ------------

  USE Type_Kinds     , ONLY: Single, Double
  USE Message_Handler, ONLY: Program_Message

  ! -----------------------
  ! Disable implicit typing
  ! -----------------------

  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------

  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'Numeric_Inquiry'
  CHARACTER(*), PARAMETER :: PROGRAM_VERSION_ID = &

  REAL(Single), PARAMETER :: SP_ONE = 1.0_Single
  REAL(Double), PARAMETER :: DP_ONE = 1.0_Double

  INTEGER, PARAMETER :: N_TEST_NUMBERS = 5
  REAL(Single), PARAMETER :: SP_TEST_NUMBER(N_TEST_NUMBERS) = &
    (/ 1.234567890123456e-16_Single, &
       1.234567890123456e-01_Single, &
       1.234567890123456e+01_Single, &
       1.234567890123456e+16_Single, &
       0.01_Single /)
  REAL(Double), PARAMETER :: DP_TEST_NUMBER(N_TEST_NUMBERS) = &
    (/ 1.234567890123456e-16_Double, &
       1.234567890123456e-01_Double, &
       1.234567890123456e+01_Double, &
       1.234567890123456e+16_Double, &
       0.01_Double /)

  ! ------------
  ! Derived type
  ! ------------

  TYPE :: SP_NI_type
    REAL(Single) :: Number
    INTEGER      :: Digits
    REAL(Single) :: Epsilon
    REAL(Single) :: Huge
    INTEGER      :: MaxExponent
    INTEGER      :: MinExponent
    INTEGER      :: Precision
    INTEGER      :: Radix
    INTEGER      :: Range
    REAL(Single) :: Tiny
    INTEGER      :: Exponent
    REAL(Single) :: Fraction
    REAL(Single) :: mNearest
    REAL(Single) :: pNearest
    REAL(Single) :: RRSpacing
    REAL(Single) :: Scale_x_0
    REAL(Single) :: Set_Exponent_x_0
    REAL(Single) :: Spacing
    REAL(Single) :: Epsilon_on_Spacing
    REAL(Single) :: Epsilon_Delta
  END TYPE SP_NI_type

  TYPE :: DP_NI_type
    REAL(Double) :: Number
    INTEGER      :: Digits
    REAL(Double) :: Epsilon
    REAL(Double) :: Huge
    INTEGER      :: MaxExponent
    INTEGER      :: MinExponent
    INTEGER      :: Precision
    INTEGER      :: Radix
    INTEGER      :: Range
    REAL(Double) :: Tiny
    INTEGER      :: Exponent
    REAL(Double) :: Fraction
    REAL(Double) :: mNearest
    REAL(Double) :: pNearest
    REAL(Double) :: RRSpacing
    REAL(Double) :: Scale_x_0
    REAL(Double) :: Set_Exponent_x_0
    REAL(Double) :: Spacing
    REAL(Double) :: Epsilon_on_Spacing
    REAL(Double) :: Epsilon_Delta
  END TYPE DP_NI_type

  ! ---------
  ! Variables
  ! ---------
  INTEGER :: i
  TYPE(SP_NI_type) :: SP_NI(N_TEST_NUMBERS)
  TYPE(DP_NI_type) :: DP_NI(N_TEST_NUMBERS)

  ! Program header
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to display numeric inquiry intrinsic results.', &
                        '$Revision$' )

  ! Output data for single type
  WRITE( *, '( /5x, "Single kind type : ", i5 )' ) Single
  DO i = 1, N_TEST_NUMBERS
    SP_NI(i)%Number = SP_TEST_NUMBER(i)
    SP_NI(i)%Digits             = DIGITS(       SP_NI(i)%Number )
    SP_NI(i)%Epsilon            = EPSILON(      SP_NI(i)%Number )
    SP_NI(i)%Huge               = HUGE(         SP_NI(i)%Number )
    SP_NI(i)%MaxExponent        = MAXEXPONENT(  SP_NI(i)%Number )
    SP_NI(i)%MinExponent        = MINEXPONENT(  SP_NI(i)%Number )
    SP_NI(i)%Precision          = PRECISION(    SP_NI(i)%Number )
    SP_NI(i)%Radix              = RADIX(        SP_NI(i)%Number )
    SP_NI(i)%Range              = RANGE(        SP_NI(i)%Number )
    SP_NI(i)%Tiny               = TINY(         SP_NI(i)%Number )
    SP_NI(i)%Exponent           = EXPONENT(     SP_NI(i)%Number )
    SP_NI(i)%Fraction           = FRACTION(     SP_NI(i)%Number )
    SP_NI(i)%mNearest           = NEAREST(      SP_NI(i)%Number, -SP_NI(i)%Number )
    SP_NI(i)%pNearest           = NEAREST(      SP_NI(i)%Number,  SP_NI(i)%Number )
    SP_NI(i)%RRSpacing          = RRSPACING(    SP_NI(i)%Number )
    SP_NI(i)%Scale_x_0          = SCALE(        SP_NI(i)%Number, 0 )
    SP_NI(i)%Set_Exponent_x_0   = SET_EXPONENT( SP_NI(i)%Number, 0 )
    SP_NI(i)%Spacing            = SPACING(      SP_NI(i)%Number )
    SP_NI(i)%Epsilon_on_Spacing = SP_NI(i)%Epsilon / SP_NI(i)%Spacing
    SP_NI(i)%Epsilon_Delta      = SP_NI(i)%Epsilon * REAL(SP_NI(i)%Radix,Single)**(SP_NI(i)%Exponent-1)
  END DO

  WRITE( *, '( /2x, "Number            = ", 10(1x,es20.13) )' ) SP_NI%Number
  WRITE( *, '(  2x, "--------------------", 10(a)          )' ) RESHAPE( (/'---------------------'/), &
                                                                         (/ N_TEST_NUMBERS /), &
                                                                         (/'---------------------'/))
  WRITE( *, '(  2x, "DIGITS()          = ", 10(16x,i5)     )' ) SP_NI%Digits
  WRITE( *, '(  2x, "EPSILON()         = ", 10(1x,es20.13) )' ) SP_NI%Epsilon
  WRITE( *, '(  2x, "HUGE()            = ", 10(1x,es20.13) )' ) SP_NI%Huge
  WRITE( *, '(  2x, "MAXEXPONENT()     = ", 10(16x,i5)     )' ) SP_NI%MaxExponent
  WRITE( *, '(  2x, "MINEXPONENT()     = ", 10(16x,i5)     )' ) SP_NI%MinExponent
  WRITE( *, '(  2x, "PRECISION()       = ", 10(16x,i5)     )' ) SP_NI%Precision
  WRITE( *, '(  2x, "RADIX()           = ", 10(16x,i5)     )' ) SP_NI%Radix
  WRITE( *, '(  2x, "RANGE()           = ", 10(16x,i5)     )' ) SP_NI%Range
  WRITE( *, '(  2x, "TINY()            = ", 10(1x,es20.13) )' ) SP_NI%Tiny
  WRITE( *, '(  2x, "EXPONENT()        = ", 10(16x,i5)     )' ) SP_NI%Exponent
  WRITE( *, '(  2x, "FRACTION()        = ", 10(1x,es20.13) )' ) SP_NI%Fraction
  WRITE( *, '(  2x, "NEAREST(-)        = ", 10(1x,es20.13) )' ) SP_NI%mNearest
  WRITE( *, '(  2x, "NEAREST(+)        = ", 10(1x,es20.13) )' ) SP_NI%pNearest
  WRITE( *, '(  2x, "RRSPACING()       = ", 10(1x,es20.13) )' ) SP_NI%RRSpacing
  WRITE( *, '(  2x, "SCALE(x,0)        = ", 10(1x,es20.13) )' ) SP_NI%Scale_x_0
  WRITE( *, '(  2x, "SET_EXPONENT(x,0) = ", 10(1x,es20.13) )' ) SP_NI%Set_Exponent_x_0
  WRITE( *, '(  2x, "SPACING()         = ", 10(1x,es20.13) )' ) SP_NI%Spacing
  WRITE( *, '(  2x, "EPSILON/SPACING   = ", 10(1x,es20.13) )' ) SP_NI%Epsilon_on_Spacing
  WRITE( *, '(  2x, "EPSILON_DELTA     = ", 10(1x,es20.13) )' ) SP_NI%Epsilon_Delta

  ! Output data for double type
  WRITE( *, '(//5x, "Double kind type : ", i5 )' ) Double
  DO i = 1, N_TEST_NUMBERS
    DP_NI(i)%Number = DP_TEST_NUMBER(i)
    DP_NI(i)%Digits             = DIGITS(       DP_NI(i)%Number )
    DP_NI(i)%Epsilon            = EPSILON(      DP_NI(i)%Number )
    DP_NI(i)%Huge               = HUGE(         DP_NI(i)%Number )
    DP_NI(i)%MaxExponent        = MAXEXPONENT(  DP_NI(i)%Number )
    DP_NI(i)%MinExponent        = MINEXPONENT(  DP_NI(i)%Number )
    DP_NI(i)%Precision          = PRECISION(    DP_NI(i)%Number )
    DP_NI(i)%Radix              = RADIX(        DP_NI(i)%Number )
    DP_NI(i)%Range              = RANGE(        DP_NI(i)%Number )
    DP_NI(i)%Tiny               = TINY(         DP_NI(i)%Number )
    DP_NI(i)%Exponent           = EXPONENT(     DP_NI(i)%Number )
    DP_NI(i)%Fraction           = FRACTION(     DP_NI(i)%Number )
    DP_NI(i)%mNearest           = NEAREST(      DP_NI(i)%Number, -DP_NI(i)%Number )
    DP_NI(i)%pNearest           = NEAREST(      DP_NI(i)%Number,  DP_NI(i)%Number )
    DP_NI(i)%RRSpacing          = RRSPACING(    DP_NI(i)%Number )
    DP_NI(i)%Scale_x_0          = SCALE(        DP_NI(i)%Number, 0 )
    DP_NI(i)%Set_Exponent_x_0   = SET_EXPONENT( DP_NI(i)%Number, 0 )
    DP_NI(i)%Spacing            = SPACING(      DP_NI(i)%Number )
    DP_NI(i)%Epsilon_on_Spacing = DP_NI(i)%Epsilon / DP_NI(i)%Spacing
    DP_NI(i)%Epsilon_Delta      = DP_NI(i)%Epsilon * REAL(DP_NI(i)%Radix,Double)**(DP_NI(i)%Exponent-1)
  END DO
  WRITE( *, '( /2x, "Number            = ", 10(1x,es27.20) )' ) DP_NI%Number
  WRITE( *, '(  2x, "--------------------", 10(a)         )' ) RESHAPE( (/'----------------------------'/), &
                                                                        (/ N_TEST_NUMBERS /), &
                                                                        (/'----------------------------'/))
  WRITE( *, '(  2x, "DIGITS()          = ", 10(23x,i5)     )' ) DP_NI%Digits
  WRITE( *, '(  2x, "EPSILON()         = ", 10(1x,es27.20) )' ) DP_NI%Epsilon
  WRITE( *, '(  2x, "HUGE()            = ", 10(1x,es27.20) )' ) DP_NI%Huge
  WRITE( *, '(  2x, "MAXEXPONENT()     = ", 10(23x,i5)     )' ) DP_NI%MaxExponent
  WRITE( *, '(  2x, "MINEXPONENT()     = ", 10(23x,i5)     )' ) DP_NI%MinExponent
  WRITE( *, '(  2x, "PRECISION()       = ", 10(23x,i5)     )' ) DP_NI%Precision
  WRITE( *, '(  2x, "RADIX()           = ", 10(23x,i5)     )' ) DP_NI%Radix
  WRITE( *, '(  2x, "RANGE()           = ", 10(23x,i5)     )' ) DP_NI%Range
  WRITE( *, '(  2x, "TINY()            = ", 10(1x,es27.20) )' ) DP_NI%Tiny
  WRITE( *, '(  2x, "EXPONENT()        = ", 10(23x,i5)     )' ) DP_NI%Exponent
  WRITE( *, '(  2x, "FRACTION()        = ", 10(1x,es27.20) )' ) DP_NI%Fraction
  WRITE( *, '(  2x, "NEAREST(-)        = ", 10(1x,es27.20) )' ) DP_NI%mNearest
  WRITE( *, '(  2x, "NEAREST(+)        = ", 10(1x,es27.20) )' ) DP_NI%pNearest
  WRITE( *, '(  2x, "RRSPACING()       = ", 10(1x,es27.20) )' ) DP_NI%RRSpacing
  WRITE( *, '(  2x, "SCALE(x,0)        = ", 10(1x,es27.20) )' ) DP_NI%Scale_x_0
  WRITE( *, '(  2x, "SET_EXPONENT(x,0) = ", 10(1x,es27.20) )' ) DP_NI%Set_Exponent_x_0
  WRITE( *, '(  2x, "SPACING()         = ", 10(1x,es27.20) )' ) DP_NI%Spacing
  WRITE( *, '(  2x, "EPSILON/SPACING   = ", 10(1x,es27.20) )' ) DP_NI%Epsilon_on_Spacing
  WRITE( *, '(  2x, "EPSILON_DELTA     = ", 10(1x,es27.20) )' ) DP_NI%Epsilon_Delta

END PROGRAM Numeric_Inquiry

!-------------------------------------------------------------------------------
!                          -- MODIFICATION HISTORY --
!-------------------------------------------------------------------------------
!
!
! $Date: 2004/08/31 19:25:54 $
!
! $Revision$
!
! $Name:  $
!
! $State: Exp $
!
! $Log: Numeric_Inquiry.f90,v $
! Revision 1.3  2004/08/31 19:25:54  paulv
! - Altered code to report data separately for Single and Double kind types.
!
! Revision 1.2  2004/08/13 20:26:28  paulv
! - Fixed a bug using the NEAREST intrinsic.
!
! Revision 1.1  2003/07/03 21:40:13  paulv
! Initial checkin.
!
!

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


! Copyright (c) 2018-2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

! Check modfile generation for external interface
module m
  interface
    integer function f(x)
    end function
    subroutine s(y, z)
      logical y
      complex z
    end subroutine
  end interface
end

!Expect: m.mod
!module m
! interface
!  function f(x)
!   real(4)::x
!   integer(4)::f
!  end
! end interface
! interface
!  subroutine s(y,z)
!   logical(4)::y
!   complex(4)::z
!  end
! end interface
!end


text:
C    @(#)mod_cbus.f	20.3 2/13/96
        subroutine mod_cbus (ic, icb)
        integer ic, icb
c
c       Perform continuation bus modification changes from chgcrd(ic)
c       onto bctbl(*,icb).
c
        include 'ipfinc/parametr.inc'

        include 'ipfinc/cbus.inc'
        include 'ipfinc/changr.inc'
        include 'ipfinc/prt.inc'
c
        character fieldfmt(7)*10, month*1
        integer field(4,7)

        data field /
     &      1, 21, 25,  2, 1, 26, 30,  3, 1, 31, 34,  4, 1, 35, 38, 5,
     &      1, 43, 47,  6, 1, 48, 52, 11, 1, 53, 57, 12 /

        data fieldfmt /
     &      '(bz, f5.0)', '(bz, f5.0)', '(bz, f4.0)', '(bz, f4.0)',
     &      '(bz, f5.0)', '(bz, f5.0)', '(bz, f5.0)' /

        read (chgcrd(ic)(122:125), '(bz, i4)') nc
c
c       Perform table driven change modifications
c
        do i = 1, 7

           i1 = field(1,i)
           i2 = field(2,i)
           i3 = field(3,i)
           i4 = field(4,i)

           if (chgcrd(ic)(i2:i3) .ne. ' ') then
              if (index (fieldfmt(i), 'f') .ne. 0) then
c
c                Floating point decoding
c
                 read (chgcrd(ic)(i2:i3), fieldfmt(i), err=900)
     &              bctbl(i4,icb)
              else
c
c                Integer or character decoding
c
                 read (chgcrd(ic)(i2:i3), fieldfmt(i), err=900)
     &              kbctbl(i4,icb)
              endif
           endif
        enddo
c
c       Implement any changes in datein
c
        if (chgcrd(ic)(75:77) .ne. ' ') then
           read (chgcrd(ic)(75:77), 140, err=900) month, kdin
  140      format (a1, i2)

           kdin = intdte(month, kdin)
           kbctbl(7, icb) = kdin
        endif

        go to 920

  900   write (errbuf(1), 910) chgcrd(ic)(1:80)
  910   format (' Illegal data in field : (',a80,')')
        call prterx ('W',1)
        chgcrd(ic)(126:126) = 'E'

  920   continue

        return
        end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
            SUBROUTINE  NAMOPN ( snbufr, snoutf, sfoutf,
     +				 snprmf, sfprmf, timstn,
     +                    	 lunin, lunsn, lunsf,
     +			         jenlun, jenstn, jensid,
     +				 nprmsn, nprmsf, sncfac, sfcfac,
     +				 parmsn, parmsf, snbprm, npbfr,
     +				 nprms2, prmsf2, sfcfc2,
     +				 tblde, nde, degree, bfropn, iret )
C************************************************************************
C* NAMOPN								*
C*									*
C* This subroutine opens the files needed by NAMSND.  The output files	*
C* are opened or created and the number of parameters is determined.	*
C* Conversion factors for the parameters are read from the parameter	*
C* packing file.  Then, the input file is opened.			*
C*									*
C* The name of the input BUFR file is given in snbufr.                  *
C*									*
C* If the surface file name is followed by +, then a second surface	*
C* file is opened.  Its name is that of the first with the suffix	*
C* _aux.  The packing table for this auxiliary file is expected to 	*
C* have the name of that specified in SFPRMF with the suffix _aux.	*
C* The auxiliary parameters are diagnosed parameters only.		*
C*									*
C* If the parameter file names (SNPRMF and SFPRMF) are blank, then	*
C* a dump of the parameter names from the BUFR table stored in the	*
C* first BUFR message is written to the file bufr_table.dump.		*
C*									*
C* NOTE ABOUT THE PARAMETER TABLE ENTRIES:				*
C*									*
C* Diagnosed or extra computed parameters are added at the end of the	*
C* parameter list for surface data.  They may be added anywhere for	*
C* profile data, but the last two columns must contain the actual	*
C* BUFR parameters in correct order.  The penultimate column is the	*
C* scaling factor, and the last column is the corresponding parameter	*
C* name expected by the conversion codes.  Note that the names in the	*
C* last column need not agree with the names in the first column for	*
C* sounding parameter files.  Surface parameter files do not have the	*
C* extra last column of names.						*
C*									*
C* An ASCII file for a particular station is requested by placing	*
C* the string |SSS=##### following the specification for SNBUFR.	*
C* SSS is the three character station ID and ##### is the station	*
C* number.  The ASCII file will be named prof.SSS.			*
C*									*
C* NAMOPN  ( SNBUFR, SNOUTF, SFOUTF, SNPRMF, SFPRMF, TIMSTN,		*
C*	     LUNIN, LUNSN, LUNSF, JENLUN, JENSTN, JENSID,		*
C*	     NPRMSN, NPRMSF, SNCFAC, SFCFAC,				*
C*	     PARMSN, PARMSF, SNBPRM, NPBFR, NPRMS2, PRMSF2, SFCFC2,	*
C*	     TBLDE, NDE, DEGREE, BFROPN, IRET )				*
C*									*
C* Input parameters:							*
C*	SNBUFR		CHAR*		Model sounding file name	*
C*	SNOUTF		CHAR*		Output sounding file name	*
C*	SFOUTF		CHAR*		Output surface file name	*
C*	SNPRMF		CHAR*		Sounding parm packing file name *
C*	SFPRMF		CHAR*		Surface parm packing file name  *
C*	TIMSTN		CHAR*		Maximum times / stations	*
C*									*
C* Output parameters:							*
C*	LUNIN		INTEGER		Model snd file unit number	*
C*	LUNSN		INTEGER		Output snd data file number	*
C*	LUNSF (2)	INTEGER		Output sfc data file numbers	*
C*	JENLUN		INTEGER		Unit number of ASCII file	*
C*	JENSTN		INTEGER		Station number for ASCII file	*
C*	JENSID		CHAR*3		Station ID for ASCII file	*
C*	NPRMSN		INTEGER		Number of sounding parms	*
C*	NPRMSF		INTEGER		Number of surface parms		*
C*	SNCFAC (NPRMSN) REAL		Conversion factors for snd data	*
C*	SFCFAC (NPRMSF) REAL		Conversion factors for sfc data *
C*	PARMSN (NPRMSN) CHAR*4		Sounding parm list		*
C*	PARMSF (NPRMSF) CHAR*4		Surface parm list		*
C*	SNBPRM (NPBFR)  CHAR*4		BUFR profile parm list		*
C*	NPBFR		INTEGER		Number in profile parm list	*
C*	NPRMS2		INTEGER		Number of auxiliary sfc parms	*
C*	PRMSF2 (NPRMS2) CHAR*4		Auxiliary surface parm list	*
C*	SFCFC2 (NPRMS2) REAL		Conversion factors for sfc data *
C*	TBLDE(*)	CHAR*		Table D entry list		*
C*	NDE		INTEGER		Number of Table D entries	*
C*	DEGREE		LOGICAL		Flag for stn lat/lon in degrees *
C*	BFROPN		LOGICAL		Flag that BUFR file is open	*
C*	IRET		INTEGER		Return code			*
C*					  3 = wrote bufr_table.dump	*
C*					  0 = normal			*
C*					 -1 = file not open		*
C*					 -3 = error opening SNBUFR      *
C*					 -4 = error opening SNPRMF	*
C*					 -5 = wrong # of snd parms	*
C*					 -6 = wrong sounding parms	*
C*					 -7 = cannot create snd file	*
C*					 -8 = error opening SFPRMF	*
C*					 -9 = wrong # of sfc parms	*
C*					-10 = wrong surface parms	*
C*					-11 = cannot create sfc file	*
C*					-13 = NAMCFC failed		*
C*					-30 = error opening aux file	*
C*					-31 = wrong # of aux parms	*
C*					-32 = wrong auxiliary parms	*
C*					-33 = cannot create aux file	*
C*					-55 = error opening SNOUTF      *
C*					-56 = error opening SFOUTF      *
C*					-57 = error opening aux SFOUTF  *
C**									*
C* Log:									*
C* K. Brill/NMC		12/93						*
C* K. Brill/NMC		05/94	Open BUFR file with DA_OPEN		*
C* K. Brill/NMC		 7/94	Added SNBPRM				*
C* K. Brill/NMC		10/94	Changes for auxiliary surface file	*
C* K. Brill/EMC		 4/96	Add SFCFC2				*
C* K. Brill/EMC		 7/98	Changes for using J. Woollen's decoder	*
C* K. Brill/EMC		 8/98	Add NPBFR, Cleanup/document error codes	*
C* K. Brill/EMC		 9/98	Added BFROPN in call			*
C* D. Kidwell/NCEP	11/98	SNEFIL -> SNBUFR; corrected prologue    *
C* K. Brill/EMC		11/98	Pass logical flag to NAMCFC		*
C* D. Kidwell/NCEP	12/98	SNM... -> NAM..., SNMODL -> NAMSND      *
C* D. Kidwell/NCEP	12/98	Check for file existence, allow ' +'    *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	snbufr, snoutf, sfoutf, snprmf, sfprmf, timstn
	CHARACTER*(*)   parmsn (*), parmsf (*), prmsf2 (*), snbprm (*)
	CHARACTER*(*)	jensid, tblde (*)
	REAL		sncfac (*), sfcfac (*), sfcfc2 (*)
	INTEGER		lunsf (2)
	LOGICAL		degree, bfropn
C*
        INTEGER         iscl (MMPARM), ioff (MMPARM), ibit (MMPARM)
	INTEGER		ilist (2)
        CHARACTER       reqprm (MMPARM)*4, filprm (MMPARM)*4
	CHARACTER	sfout*256, sfout2*256, sfprm2*256, newfil*256
	CHARACTER	snef(2)*132, bfrfl*132, csn (2)*16
	LOGICAL		pkflg, mrgdat, exist
C----------------------------------------------------------------------
	iret  = 0
	bfropn = .false.
C*
	nde = 0
	lunsn = 0
	lunsf (1) = 0
	lunsf (2) = 0
	nprmsn = 0
	nprmsf = 0
	nprms2 = 0
	jenlun = 0
	jenstn = 0
	jensid = ' '
	CALL ST_LSTR ( snprmf, lnsn, ier )
	CALL ST_LSTR ( sfprmf, lnsf, ier )
C
C*	Find out if an ASCII output file is needed.
C
	iq = INDEX ( snbufr, '|' )
	IF ( iq .ne. 0 ) THEN
	    CALL ST_CLST ( snbufr, '|', ' ', 2, snef, num, ier )
	    IF ( lnsf .eq. 0 .and. lnsn .eq. 0 ) THEN
		CALL NAMINV ( snef (1), iret )
		IF ( iret .ne. 0 ) THEN
		    CALL ER_WMSG ( 'NAMSND', iret, ' ', ier )
		END IF
		iret = 3
		RETURN
	    END IF
	    CALL ST_CLST ( snef (2), '=', ' ', 2, csn, num, ier )
	    jensid = csn (1) (1:3)
	    CALL ST_NUMB ( csn (2), jenstn, ier )
	    snef (2) = 'prof.' // jensid
	    CALL ST_LCUC ( jensid, jensid, ier )
	    CALL FL_SUNK ( snef (2), jenlun, ier )
	    IF ( snoutf .eq. ' ' ) THEN
C
C*      	Get parameter information
C
        	CALL DP_FILE ( snprmf, npsn, filprm, iscl, ioff,
     +			       ibit, pkflg, iret )
	    	nprmsn = npsn
	    	DO ip = 1, nprmsn
	            parmsn (ip) = filprm (ip)
	    	END DO
	        CALL NAMCFC ( snprmf, .true.,
     +			      npbfr, sncfac, snbprm, iret )
		IF ( iret .ne. 0 ) RETURN
	    END IF
	    IF ( sfoutf .eq. ' ' ) THEN
	    	CALL DP_FILE ( sfprmf, npsf, filprm, iscl,
     +			       ioff, ibit, pkflg, ier )
	    	nprmsf = npsf
	    	DO ip = 1, nprmsf
	            parmsf (ip) = filprm (ip)
	    	END DO
	    	CALL NAMCFC ( sfprmf, .false.,
     +			      n, sfcfac, reqprm, iret )
		IF ( iret .ne. 0 ) RETURN
	    END IF
	ELSE
	    snef (1) = snbufr
	    IF ( lnsf .eq. 0 .and. lnsn .eq. 0 ) THEN
		CALL NAMINV ( snef (1), iret )
		IF ( iret .ne. 0 ) THEN
		    CALL ER_WMSG ( 'NAMSND', iret, ' ', ier )
		END IF
		iret = 3
		RETURN
	    END IF
	END IF
C
C*	Open the sounding output file.
C
	IF ( snoutf .ne. ' ' ) THEN
	    CALL FL_INQR ( snoutf, exist, newfil, ier )
	    IF ( exist ) THEN
	        CALL SN_OPNF ( snoutf, .true., lunsn, iflsrc, npsn,
     +			       filprm, ivert, mrgdat, ier )
	        IF ( ier .eq. 0 ) THEN
	    	    CALL DP_FILE ( snprmf, npreq, reqprm, iscl,
     +			           ioff, ibit, pkflg, ier )
	    	    IF ( ier .ne. 0 ) THEN
		        iret = -4
		        RETURN
	    	    END IF
	    	    IF ( npsn .ne. npreq ) THEN
		        iret = -5
		        RETURN
	    	    END IF
	    	    DO i = 1, npsn
		        IF ( filprm (i) .ne. reqprm (i) ) THEN
		    	    iret = -6
		    	    RETURN
		        END IF
	    	    END DO
		  ELSE
		    iret = -55
		    RETURN
		END IF
	    ELSE
C
C*	        Create a new sounding data file.
C
C*          	Get the number of times and the number of stations to
C*          	add.
C
            	CALL ST_ILST  ( timstn, '/', IMISSD, 2, ilist, n,
     +				ier )
            	IF  ( ilist (1) .eq. IMISSD )  THEN
                    maxtim = 49
            	ELSE
                    maxtim = ilist (1)
            	END IF
            	IF  ( ilist (2) .eq. IMISSD )  THEN
                    maxstn = 700
            	ELSE
                    maxstn = ilist (2)
            	END IF
	    	CALL SN_CRFP ( snoutf, snprmf, 0, maxstn, maxtim,
     +			       .true., lunsn, npsn, filprm, pkflg,
     +			       ier )
	    	IF ( ier .ne. 0 ) THEN
		    iret = -7
		    RETURN
	    	END IF
	    	CALL ER_WMSG ( 'NAMSND', 2, snoutf, ier )
	    END IF
	    nprmsn = npsn
	    DO ip = 1, nprmsn
	        parmsn (ip) = filprm (ip)
	    END DO
	    CALL NAMCFC ( snprmf, .true.,
     +			  npbfr, sncfac, snbprm, iret )
	    IF ( iret .ne. 0 ) RETURN
	END IF
C
C*	Open the surface output file(s).
C
	iplus = INDEX ( sfoutf, '+' )
	IF ( iplus .gt. 1 ) THEN
	    im = iplus - 1
	    sfout = sfoutf ( 1:im )
	    CALL ST_LSTR ( sfout ( 1:im ), lens, ier )
	    sfout2 = sfout ( 1:lens ) // '_aux'
	    CALL ST_LSTR ( sfprmf, lens, ier )
	    sfprm2 = sfprmf ( 1:lens ) // '_aux'
	ELSE
	    sfout = sfoutf
	    sfout2 = ' '
	END IF
	IF ( sfout .ne. ' ' ) THEN
	    CALL FL_INQR ( sfout, exist, newfil, ier )
	    IF ( exist ) THEN
	        CALL SF_OPNF ( sfout, .true., lunsf (1), iflsrc,
     +			       npsf, filprm, ier )
	        IF ( ier .eq. 0 ) THEN
	    	    CALL DP_FILE ( sfprmf, npreq, reqprm, iscl,
     +			           ioff, ibit, pkflg, ier )
	    	    IF ( ier .ne. 0 ) THEN
		        iret = -8
		        RETURN
	    	    END IF
	    	    IF ( npsf .ne. npreq ) THEN
		        iret = -9
		        RETURN
	    	    END IF
	    	    DO i = 1, npsf
		        IF ( filprm (i) .ne. reqprm (i) ) THEN
		    	    iret = -10
		    	    RETURN
		        END IF
	    	    END DO
		  ELSE
		    iret = -56
		    RETURN
		END IF
	    ELSE
C
C*	    	Create a new surface data file.
C
C*          	Get the number of times and the number of stations to
C*          	add.
C
            	CALL ST_ILST  ( timstn, '/', IMISSD, 2, ilist,
     +				n, ier )
            	IF  ( ilist (1) .eq. IMISSD )  THEN
                    maxtim = 49
            	ELSE
                    maxtim = ilist (1)
            	END IF
            	IF  ( ilist (2) .eq. IMISSD )  THEN
                    maxstn = 700
            	ELSE
                    maxstn = ilist (2)
            	END IF
	    	CALL SF_CRFP ( sfout, sfprmf, 0, maxstn, maxtim,
     +			       .true., lunsf (1), npsf, filprm, pkflg,
     +			       ier )
	    	IF ( ier .ne. 0 ) THEN
		    iret = -11
		    RETURN
	    	END IF
	    	CALL ER_WMSG ( 'NAMSND', 2, sfout, ier )
	    END IF
	    nprmsf = npsf
	    DO ip = 1, nprmsf
	        parmsf (ip) = filprm (ip)
	    END DO
	    CALL NAMCFC ( sfprmf, .false.,
     +			  n, sfcfac, reqprm, iret )
	    IF ( iret .ne. 0 ) RETURN
	END IF
C
C*	Open a second file if needed.
C
	IF ( sfout2 .ne. ' ' ) THEN
	    CALL FL_INQR ( sfout2, exist, newfil, ier )
	    IF ( exist ) THEN
	        CALL SF_OPNF ( sfout2, .true., lunsf (2), iflsrc,
     +			       npsf, filprm, ier )
	        IF ( ier .eq. 0 ) THEN
	    	    CALL DP_FILE ( sfprm2, npreq, reqprm, iscl,
     +			           ioff, ibit, pkflg, ier )
	    	    IF ( ier .ne. 0 ) THEN
		        iret = -30
		        RETURN
	    	    END IF
	    	    IF ( npsf .ne. npreq ) THEN
		        iret = -31
		        RETURN
	    	    END IF
	    	    DO i = 1, npsf
		        IF ( filprm (i) .ne. reqprm (i) ) THEN
		    	    iret = -32
		    	    RETURN
		        END IF
	    	    END DO
		  ELSE
		    iret = -57
		    RETURN
		END IF
	    ELSE
C
C*	    	Create a new surface data file.
C
C*          	Get the number of times and the number of stations to
C*          	add.
C
            	CALL ST_ILST  ( timstn, '/', IMISSD, 2, ilist,
     +				n, ier )
            	IF  ( ilist (1) .eq. IMISSD )  THEN
                    maxtim = 49
            	ELSE
                    maxtim = ilist (1)
            	END IF
            	IF  ( ilist (2) .eq. IMISSD )  THEN
                    maxstn = 700
            	ELSE
                    maxstn = ilist (2)
            	END IF
	    	CALL SF_CRFP ( sfout2, sfprm2, 0, maxstn, maxtim,
     +			       .true., lunsf (2), npsf, filprm, pkflg,
     +			       ier )
	    	IF ( ier .ne. 0 ) THEN
		    iret = -33
		    RETURN
	    	END IF
	    	CALL ER_WMSG ( 'NAMSND', 2, sfout2, ier )
	    END IF
	    nprms2 = npsf
	    DO ip = 1, nprms2
	        prmsf2 (ip) = filprm (ip)
	    END DO
	    CALL NAMCFC ( sfprm2, .false.,
     +			  n, sfcfc2, reqprm, iret )
	    IF ( iret .ne. 0 ) RETURN
	END IF
C
C*	Open the input file.
C
	bfrfl = snef (1)
	CALL SS_ENVR ( bfrfl, bfrfl, ier )
        CALL NAMPRP ( bfrfl, lunin, tblde, nde, degree, ier )
	IF ( ier .ne. 0 ) THEN
	    CALL ER_WMSG ( 'NAMSND', ier, ' ', ier2 )
	    iret = -3
	    RETURN
	END IF
	bfropn = .true.
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
*
*     ------------------------------------------------------------------
*	C F P F
*     ------------------------------------------------------------------
*
      SUBROUTINE CFPF(N,IVI,ILI,ISI,IVJ,ILJ,ISJ,COEFP)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
*
*     THIS IS A DUMMY SUBROUTINE TO CALCULATE CFP OF F-ELECTRONS.  IT IS
*     VALID ONLY FOR ONE OR TWO ELECTRONS IN THE F-SHELL UNDER
*     CONSIDERATION.
*
      COEFP=1.D0
      RETURN
      END


text:
* =foo0.f in Burley's g77 test suite.
      subroutine sub(a)
      common /info/ iarray(1000)
      equivalence (m,iarray(100)), (n,iarray(200))
      real a(m,n)
      a(1,1) = a(2,2)
      end


text:
   ! This is a quick test of hyperBOB using the Eggholder function
   ! see https://en.wikipedia.org/wiki/Test_functions_for_optimization
   ! This uses Ver 1.3 of hyperBOB which allows xMat_hyperBOB to be bigger than needed.

   program EggholderTest
      use hyperBOB_module
      use mpi
      implicit none
      integer mpierr, myRank, nunit, nBOBs, i
      integer nx, npt, iprint, maxFun, seed, err
      real(wp) rhoBMin, rhoEMax, fMin, maxTime
      real(wp), allocatable:: x(:), xL(:), xU(:)
      real(wp) comp_start, comp_end
      character(len=40) fmtString
      external EggholderFcn

      ! Initialize MPI
      call MPI_INIT( mpierr )
      call MPI_COMM_RANK(MPI_COMM_WORLD, myRank, mpiErr) ! get rank of this process in world
      if (myRank == 0) then
         write(*,'(/,a)') 'Starting Eggholder Test'
         write(*,'(a,$)') 'Enter maxFun evaluations: '
         read(*,*) maxFun
      end if
      ! boradcast maxFun to all processes
      call MPI_BCAST(maxFun, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, mpierr)
      nx = 2
      ! allocate space
      allocate( x(nx), xL(nx), xU(nx) )
      npt = 2*nx + 1
      xL = -512.; xU = 512.
      rhoBMin = 512.
      rhoEMax = 0.001
      iprint = 0
      !maxFun = 2000
      maxTime = 0 ! this is not implemented yet
      seed = 5 ! random number start
      comp_start = MPI_Wtime()
      ! intialize hyperBOB
      call hyperBOB_initialize(nx+5)   ! add 5 extra rows
      call hyperBOB(nx, npt, x, xL, xU, rhoBMin, rhoEMax, iprint, maxFun, maxTime, EggholderFcn, seed, fMin, err)
      ! Save xMat
      if (myRank == 0) then
         open(newunit=nunit, file='allSolns.dat',status='unknown')
         fmtString = '(i4, #######(1x,g13.5))'
         nBOBs = size(xMat_hyperBOB, dim=2)
         write(fmtString(6:12),'(i7)') nBobs
         write(nunit, fmtString) 0, fVal_hyperBOB(1:nBobs)
         do i=1,nx
            write(nunit, fmtString) i, xMat_hyperBOB(i,1:nBobs)
         end do
         close(unit=nunit)
      end if
      ! deallocate space in hyperBOB
      call hyperBOB_cleanUp()

      comp_end = MPI_Wtime()

      if(myRank == 0) then
         write(*,'(a,g13.5)') 'fMin = ', fMin
         write(*,'(a,2(1x,g13.5))') 'Solution: ', x(1), x(2)
         write(*,'(a,g13.5)') 'Total time (min): ', (comp_end-comp_start)/60.
      end if

      call MPI_FINALIZE(mpierr)

      deallocate( x, xL, xU )
      stop

   end program EggholderTest

   subroutine EggholderFcn(n, xy, f)
      use hyperBOB_module
      ! This function's solution: f(512, 404.2319) = -959.6407
      integer, intent(in)              :: n  ! deminsion of problem
      real(wp),dimension(:),intent(in) :: xy ! Value of x were f is to be evaluated at
      real(wp), intent(out)            :: f  ! function value
      ! local declarations
      integer i
      real(wp) x, y
      x = xy(1)
      y = xy(2)

      f = -(y + 47.)*sin(sqrt(abs(x/2. + (y + 47.)))) - x*sin(sqrt(abs(x - (y + 47.))))

      return
   end subroutine EggholderFcn

max_stars_repo_name: maxEntropyProd/hyperBOB

max_issues_repo_name: maxEntropyProd/hyperBOB

max_forks_repo_name: maxEntropyProd/hyperBOB


text:
PROGRAM Decoherence

  USE KindModule, ONLY: &
    DP, Zero, One, Two, Pi, TwoPi
  USE UnitsModule, ONLY: &
    Kilometer, MeV
  USE ProgramHeaderModule, ONLY: &
    iX_B0, iX_E0, iX_B1, iX_E1, &
    iE_B0, iE_E0, iE_B1, iE_E1, &
    iZ_B0, iZ_E0, iZ_B1, iZ_E1
  USE ProgramInitializationModule, ONLY: &
    InitializeProgram, &
    FinalizeProgram
  USE ReferenceElementModuleX, ONLY: &
    InitializeReferenceElementX, &
    FinalizeReferenceElementX
  USE ReferenceElementModuleX_Lagrange, ONLY: &
    InitializeReferenceElementX_Lagrange, &
    FinalizeReferenceElementX_Lagrange
  USE ReferenceElementModuleE, ONLY: &
    InitializeReferenceElementE, &
    FinalizeReferenceElementE
  USE ReferenceElementModuleE_Lagrange, ONLY: &
    InitializeReferenceElementE_Lagrange, &
    FinalizeReferenceElementE_Lagrange
  USE ReferenceElementModule, ONLY: &
    InitializeReferenceElement, &
    FinalizeReferenceElement
  USE ReferenceElementModule_Lagrange, ONLY: &
    InitializeReferenceElement_Lagrange, &
    FinalizeReferenceElement_Lagrange
  USE GeometryComputationModule, ONLY: &
    ComputeGeometryX
  USE GeometryComputationModuleE, ONLY: &
    ComputeGeometryE
  USE GeometryFieldsModule, ONLY: &
    uGF
  USE GeometryFieldsModuleE, ONLY: &
    uGE
  USE InputOutputModuleHDF, ONLY: &
    WriteFieldsHDF
  USE EquationOfStateModule, ONLY: &
    InitializeEquationOfState, &
    FinalizeEquationOfState
  USE TwoMoment_ClosureModule, ONLY: &
    InitializeClosure_TwoMoment
  USE InitializationModule, ONLY: &
    InitializeFields_Decoherence

  IMPLICIT NONE

  INTEGER       :: nNodes
  INTEGER       :: nE, bcE, nX(3), bcX(3)
  REAL(DP)      :: eL, eR, xL(3), xR(3)

  nX  = [ 128, 1, 1 ]
  xL  = [ 1.0d0 * Kilometer, 0.0_DP, 0.0_DP ]
  xR  = [ 1.0d2 * Kilometer, Pi,     TwoPi  ]
  bcX = [ 1, 0, 0 ]

  nE  = 40
  eL  = 0.0d0 * MeV
  eR  = 2.5d2 * MeV
  bcE = 0

  nNodes = 1

  CALL InitializeProgram &
         ( ProgramName_Option &
             = 'Decoherence', &
           nX_Option &
             = nX, &
           swX_Option &
             = [ 1, 1, 1 ], &
           bcX_Option &
             = bcX, &
           xL_Option &
             = xL, &
           xR_Option &
             = xR, &
           nE_Option &
             = nE, &
           swE_Option &
             = 0, &
           bcE_Option &
             = bcE, &
           eL_Option &
             = eL, &
           eR_Option &
             = eR, &
           nNodes_Option &
             = nNodes, &
           CoordinateSystem_Option &
             = 'SPHERICAL', &
           ActivateUnits_Option &
             = .TRUE., &
           BasicInitialization_Option &
             = .TRUE. )

  ! --- Position Space Reference Element and Geometry ---

  CALL InitializeReferenceElementX

  CALL InitializeReferenceElementX_Lagrange

  CALL ComputeGeometryX &
         ( iX_B0, iX_E0, iX_B1, iX_E1, uGF )

  ! --- Energy Space Reference Element and Geometry ---

  CALL InitializeReferenceElementE

  CALL InitializeReferenceElementE_Lagrange

  CALL ComputeGeometryE &
         ( iE_B0, iE_E0, iE_B1, iE_E1, uGE )

  ! --- Phase Space Reference Element ---

  CALL InitializeReferenceElement

  CALL InitializeReferenceElement_Lagrange

  ! --- Initialize Equation of State ---

  CALL InitializeEquationOfState &
         ( EquationOfState_Option &
             = 'TABLE', &
           EquationOfStateTableName_Option &
             = 'EquationOfStateTable.h5' )

  ! --- Initialize Moment Closure ---

  CALL InitializeClosure_TwoMoment

  CALL InitializeFields_Decoherence

  CALL WriteFieldsHDF &
         ( Time = 0.0_DP, &
           WriteGF_Option = .TRUE., &
           WriteFF_Option = .TRUE., &
           WriteRF_Option = .TRUE. )

  !!! Stuff will go here

  CALL FinalizeEquationOfState

  CALL FinalizeReferenceElementX

  CALL FinalizeReferenceElementX_Lagrange

  CALL FinalizeReferenceElementE

  CALL FinalizeReferenceElementE_Lagrange

  CALL FinalizeReferenceElement

  CALL FinalizeReferenceElement_Lagrange

  CALL FinalizeProgram

END PROGRAM Decoherence

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


text:
      SUBROUTINE MF_FINLLG
     &          (xintg,nderiv)

      implicit complex (a-h,o-z)

      real           omega,wavenr,
     &               delta(4)

      common/mf_flag/iexact,iovflo,kexact,lowg,nodivd,noevct,
     &               nofinl,nointg,nomesh,notlin
     &      /mf_lgfl/lgflag(4)
     &      /mf_mode/theta,c,s,csq,ssq,omega,wavenr,ideriv
     &      /mf_rmtx/x(4),dxdc(4),dxdh(4),dhdxdc(4)

      dimension      xintg(8),r(4),drdc(4)

      data           delta/1.0,0.0,0.0,1.0/

      do i=1,4
         if (lgflag(i) .eq. 1) then
            r(i)=EXP(xintg(i))
            x(i)=(r(i)+delta(i))/c
            xintg(i)=x(i)
            if (nderiv .eq. 1) then
               drdc(i)=r(i)*xintg(i+4)
               dxdc(i)=(drdc(i)-x(i))/c
               xintg(i+4)=dxdc(i)
            end if
         end if
      end do
      RETURN
      END      ! MF_FINLLG


text:
! Calculate the exponential integral using the algorithm of
! [1] Tseng, Lee, Journal of Hydrology, 205 (1998) 38-51
module m_exponential_integral

   implicit none

   real, parameter :: series_laguerre = 4.0

contains

   ! Calculate the exponential integral E_1(x):
   !
   !        inf
   !          /     -t
   !         |     e
   ! E (x) = | dt -----
   !  1      |      t
   !        /
   !         x
   !
   ! Input:  arg - position at which exponential integral is evaluated (arg > 0)
   ! Output: res - E_1(arg)
   pure subroutine calculateExponentialIntegral(arg, res)

      implicit none

      real, intent(in)  :: arg
      real, intent(inout) :: res

      ! For arguments smaller than 4 the series expansion is used
      if (arg < series_laguerre) then
         res = seriesExpansion(arg)

         ! otherwise a Gauss-Laguerre expansion is better
      else
         res = exp(-arg)*gauss_laguerre(arg)
      endif

   end subroutine calculateExponentialIntegral

   ! Series expansion of the exponential integral
   !
   !                          n_cut
   !                          -----     n  n
   !                           \    (-1)  x
   ! E (x) = -gamma - ln(x) -   )   --------
   !  1                        /     n * n!
   !                          -----
   !                          n = 1
   !
   ! where gamma is the Euler constant.
   ! n_cut is set to 25
   ! Input: arg - argument for which the exponential integral is approximated
   ! Return: approximation by series expansion for E_1(arg)
   pure real function seriesExpansion(arg)

      implicit none

      real, intent(in) :: arg

      real    :: res, fact  ! result of the summation, 1 / n
      integer :: i          ! counter variable

      real, parameter :: EULER_GAMMA = 0.57721566490153286060651209008241 ! Euler constant
      integer, parameter :: ITERATION = 25 ! Cutoff for series expansion

      ! initialize summation result
      res = 0.0

      ! perform the summation
      do i = ITERATION, 2, -1
         ! calculate 1/n
         fact = 1.0/i
         ! add next term of summation
         res = arg*fact*(fact - res)
      end do

      ! calculate the final result
      seriesExpansion = -EULER_GAMMA - log(arg) + arg*(1.0 - res)

   end function seriesExpansion

   ! The Gauss Laguerre expansion of the exponential integral can be written as
   !
   !              N
   ! E (arg)    -----    a
   !  1          \        n
   ! ------- =    )   --------
   !   -arg      /    x  + arg
   !  e         -----  n
   !             n=1
   !
   ! where the a_n and x_n are determined by least quadrature and are given in [1]
   ! Input: arg - point at which Gaussian Laguerre quadrature is calculated
   ! Return: E_1(arg) in this approximation
   pure real function gauss_laguerre(arg)

      implicit none

      real, intent(in) :: arg

      ! the quadrature constants a_n and x_n from [1]
      real, parameter :: a(1:15) = (/ &
                         0.2182348859400869e+00, 0.3422101779228833e+00, 0.2630275779416801e+00, &
                         0.1264258181059305e+00, 0.4020686492100091e-01, 0.8563877803611838e-02, &
                         0.1212436147214252e-02, 0.1116743923442519e-03, 0.6459926762022901e-05, &
                         0.2226316907096273e-06, 0.4227430384979365e-08, 0.3921897267041089e-10, &
                         0.1456515264073126e-12, 0.1483027051113301e-15, 0.1600594906211133e-19/)
      real, parameter :: x(1:15) = (/ &
                         0.9330781201728180e-01, 0.4926917403018839e+00, 0.1215595412070949e+01, &
                         0.2269949526203743e+01, 0.3667622721751437e+01, 0.5425336627413553e+01, &
                         0.7565916226613068e+01, 0.1012022856801911e+02, 0.1313028248217572e+02, &
                         0.1665440770832996e+02, 0.2077647889944877e+02, 0.2562389422672878e+02, &
                         0.3140751916975394e+02, 0.3853068330648601e+02, 0.4802608557268579e+02/)

      ! Calculate the summation
      gauss_laguerre = sum(a/(x + arg))

   end function gauss_laguerre

end module m_exponential_integral

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
C    Copyright(C) 2014-2017 National Technology & Engineering Solutions of
C    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    Redistribution and use in source and binary forms, with or without
C    modification, are permitted provided that the following conditions are
C    met:
C
C    * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C    * Redistributions in binary form must reproduce the above
C      copyright notice, this list of conditions and the following
C      disclaimer in the documentation and/or other materials provided
C      with the distribution.
C
C    * Neither the name of NTESS nor the names of its
C      contributors may be used to endorse or promote products derived
C      from this software without specific prior written permission.
C
C    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C    A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C    OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C $Id: undelm.f,v 1.2 1998/07/14 18:20:13 gdsjaar Exp $
C $Log: undelm.f,v $
C Revision 1.2  1998/07/14 18:20:13  gdsjaar
C Removed unused variables, cleaned up a little.
C
C Changed BLUE labels to GREEN to help visibility on black background
C (indirectly requested by a couple users)
C
C Revision 1.1.1.1  1990/11/30 11:17:29  gdsjaar
C FASTQ Version 2.0X
C
c Revision 1.1  90/11/30  11:17:27  gdsjaar
c Initial revision
c
C
CC* FILE: [.PAVING]UNDELM.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE UNDELM (MXND, MLN, LNODES, XN, YN, NUID, LXK, KXL, NXL,
     &   LXN, NNN, LLL, KKK, NAVAIL, IAVAIL, N0, N1, N2, N3, L1, L2, L3,
     &   K1, K2, NOROOM, ERR, GRAPH, VIDEO)
C***********************************************************************
C
C  SUBROUTINE UNDELM = UNDELETES AN ELEMENT BY EXPANDING N1 INTO A
C                      NEW ELEMENT
C
C***********************************************************************
C
      DIMENSION LXK(4, MXND), NXL(2, 3*MXND), KXL(2, 3*MXND)
      DIMENSION LXN(4, MXND), XN(MXND), YN(MXND), NUID(MXND)
      DIMENSION LNODES (MLN, MXND)
C
      LOGICAL ERR, NOROOM, GRAPH, VIDEO
C
      ERR = .FALSE.
C
C  MAKE SURE THAT N2 HAS AT LEAST FOUR LINES ATTACHED TO IT
C
      IF (LXN (4, N2) .EQ. 0) THEN
         ERR = .TRUE.
         CALL MESAGE ('** N2 IN UNDELM CANNOT BE USED'//
     &      ' TO EXPAND AN ELEMNT **')
         GOTO 140
      ENDIF
C
C  ERASE THE LINE L3
C
      IF ((GRAPH) .OR. (VIDEO)) THEN
         CALL LCOLOR ('BLACK')
         CALL D2NODE (MXND, XN, YN, N0, N2)
         IF (GRAPH) THEN
            CALL LCOLOR ('WHITE')
         ELSE
            CALL LCOLOR ('YELOW')
         ENDIF
         CALL SFLUSH
         IF (VIDEO) CALL SNAPIT (3)
      ENDIF
C
C  DEFINE THE NEW NODE AND THE TWO NEW LINES
C
      NNN = NNN + 1
      IF (NNN .GT. MXND) THEN
         NOROOM = .TRUE.
         GOTO 140
      ENDIF
      NNEW = NNN
      XN (NNEW) = (XN (N0) + XN (N2)) * .5
      YN (NNEW) = (YN (N0) + YN (N2)) * .5
C
      LLL = LLL + 2
      L4 = LLL -1
      L5 = LLL
      NXL (1, L4) = N1
      NXL (2, L4) = NNEW
      NXL (1, L5) = N3
      NXL (2, L5) = NNEW
C
C  NOW CHANGE LINE L3'S END POINT FROM N2 TO NNEW
C
      IF (NXL (1, L3) .EQ. N2) THEN
         NXL (1, L3) = NNEW
      ELSEIF (NXL (2, L3) .EQ. N2) THEN
         NXL (2, L3) = NNEW
      ELSE
         CALL MESAGE ('** PROBLEMS IN UNDLEM WITH L3''S END POINT **')
         ERR = .TRUE.
         GOTO 140
      ENDIF
C
C  NOW UPDATE THE LXN ARRAYS
C
      LXN (1, NNEW) = L3
      LXN (2, NNEW) = L4
      LXN (3, NNEW) = L5
      LXN (4, NNEW) = 0
      CALL FIXLXN (MXND, LXN, NXL, NUID, NAVAIL, IAVAIL, NNN, LLL,
     &   NNN, LLL, ERR, NOROOM)
      IF ((NOROOM) .OR. (ERR)) GOTO 140
C
C  REMOVE L3 FROM THE LIST OF LINES FOR N2
C
      CALL DELLXN (MXND, LXN, NUID, NAVAIL, IAVAIL, N2,
     &   L3, NNN, ERR, NOROOM)
      IF ((NOROOM) .OR. (ERR)) THEN
         CALL MESAGE ('** PROBLEMS IN UNDELM UNHOOKING L3 FROM N2 **')
         GOTO 140
      ENDIF
C
C  ADD LINE L4 TO N1
C
      CALL ADDLXN (MXND, LXN, NUID, NAVAIL, IAVAIL,
     &   N1, L4, NNN, ERR, NOROOM)
      IF ((NOROOM) .OR. (ERR)) THEN
         CALL MESAGE ('** PROBLEMS IN UNDELM HOOKING L4 TO N1 **')
         GOTO 140
      ENDIF
C
C  ADD LINE L5 TO N3
C
      CALL ADDLXN (MXND, LXN, NUID, NAVAIL, IAVAIL,
     &   N3, L5, NNN, ERR, NOROOM)
      IF ((NOROOM) .OR. (ERR)) THEN
         CALL MESAGE ('** PROBLEMS IN UNDELM HOOKING L5 TO N3 **')
         GOTO 140
      ENDIF
C
C  NOW ADD THE NEW ELEMENT
C
      KKK = KKK + 1
      LXK (1, KKK) = L1
      LXK (2, KKK) = L2
      LXK (3, KKK) = L5
      LXK (4, KKK) = L4
C
C  NOW FIX THE KXL ARRAY FOR LINE L1
C
      IF (KXL (1, L1) .EQ. K2) THEN
         KXL (1, L1) = KKK
      ELSEIF (KXL (2, L1) .EQ. K2) THEN
         KXL (2, L1) = KKK
      ELSE
         CALL MESAGE ('** PROBLEMS IN UNDELM REPLACING K2 FOR L1 **')
         ERR = .TRUE.
         GOTO 140
      ENDIF
C
C  NOW FIX THE KXL ARRAY FOR LINE L2
C
      IF (KXL (1, L2) .EQ. K1) THEN
         KXL (1, L2) = KKK
      ELSEIF (KXL (2, L2) .EQ. K1) THEN
         KXL (2, L2) = KKK
      ELSE
         CALL MESAGE ('** PROBLEMS IN UNDELM REPLACING K1 FOR L2 **')
         ERR = .TRUE.
         GOTO 140
      ENDIF
C
C  ADD THE KXL ENTRIES FOR THE NEW LINES
C
      KXL (1, L4) = K2
      KXL (2, L4) = KKK
      KXL (1, L5) = K1
      KXL (2, L5) = KKK
C
C  NOW FIX THE LXK ARRAY FOR THE ELEMENT K1
C
      DO 100 I = 1, 4
         IF (LXK (I, K1) .EQ. L2) THEN
            LXK (I, K1) = L5
            GOTO 110
         ENDIF
  100 CONTINUE
      CALL MESAGE ('** PROBLEMS IN UNDELM REPLACING L2 WITH L5 IN '//
     &   'K1 **')
      ERR = .TRUE.
      GOTO 140
  110 CONTINUE
C
C  NOW FIX THE LXK ARRAY FOR THE ELEMENT K2
C
      DO 120 I = 1, 4
         IF (LXK (I, K2) .EQ. L1) THEN
            LXK (I, K2) = L4
            GOTO 130
         ENDIF
  120 CONTINUE
      CALL MESAGE ('** PROBLEMS IN UNDELM REPLACING L1 WITH L4 IN '//
     &   'K2 **')
      ERR = .TRUE.
      GOTO 140
  130 CONTINUE
C
C  NOW REDRAW THE LINES
C
      IF ((GRAPH) .OR. (VIDEO)) THEN
         CALL D2NODE (MXND, XN, YN, N0, NNEW)
         CALL D2NODE (MXND, XN, YN, N1, NNEW)
         CALL D2NODE (MXND, XN, YN, N3, NNEW)
         CALL SFLUSH
         IF (VIDEO) CALL SNAPIT (3)
      ENDIF
C
      LNODES (4, NNEW) = 2
      CALL MARKSM (MXND, MLN, LXK, KXL, NXL, LXN, LNODES,
     &   NNEW, ERR)
      IF (ERR) GOTO 140
      CALL MARKSM (MXND, MLN, LXK, KXL, NXL, LXN, LNODES,
     &   N0, ERR)
      IF (ERR) GOTO 140
      CALL MARKSM (MXND, MLN, LXK, KXL, NXL, LXN, LNODES,
     &   N1, ERR)
      IF (ERR) GOTO 140
      CALL MARKSM (MXND, MLN, LXK, KXL, NXL, LXN, LNODES,
     &   N2, ERR)
      IF (ERR) GOTO 140
      CALL MARKSM (MXND, MLN, LXK, KXL, NXL, LXN, LNODES,
     &   N3, ERR)
      IF (ERR) GOTO 140
  140 CONTINUE
      RETURN
C
      END

max_stars_repo_name: jbakosi/seacas

max_issues_repo_name: jbakosi/seacas

max_forks_repo_name: jbakosi/seacas


text:
!==============================================================================!
  subroutine Balance_Mass(flow)
!------------------------------------------------------------------------------!
!   Modifies the fluxes at outflow boundaries to conserve the mass.            !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Const_Mod
  use Comm_Mod
  use Grid_Mod,  only: Grid_Type, Grid_Mod_Bnd_Cond_Type,  &
                       INFLOW, OUTFLOW, CONVECT, PRESSURE
  use Field_Mod, only: Field_Type, Field_Mod_Alias_Momentum
  use Var_Mod,   only: Var_Type
  use Face_Mod,  only: Face_Type
  use Bulk_Mod,  only: Bulk_Type
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Field_Type), target :: flow
!-----------------------------------[Locals]-----------------------------------!
  type(Grid_Type), pointer :: grid
  type(Bulk_Type), pointer :: bulk
  type(Var_Type),  pointer :: u, v, w
  type(Face_Type), pointer :: m_flux
  integer                  :: s, c1, c2
  real                     :: fac
!==============================================================================!

  ! Take aliases
  grid   => flow % pnt_grid
  bulk   => flow % bulk
  m_flux => flow % m_flux
  call Field_Mod_Alias_Momentum(flow, u, v, w)

  !--------------------------------------!
  !   Calculate the inflow mass fluxes   !
  !--------------------------------------!
  bulk % mass_in = 0.0
  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2 < 0) then
      m_flux % n(s) = flow % density_f(s) * ( u % n(c2)*grid % sx(s)    &
                                            + v % n(c2)*grid % sy(s)    &
                                            + w % n(c2)*grid % sz(s) )

      if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. INFLOW) then
        bulk % mass_in = bulk % mass_in - m_flux % n(s)
      end if

      if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. PRESSURE  &
         .and. m_flux % n(s) < 0.0) then
        bulk % mass_in = bulk % mass_in - m_flux % n(s)
      end if

      if(Grid_Mod_Bnd_Cond_Type(grid,c2) .eq. CONVECT  &
         .and. m_flux % n(s) < 0.0) then
        bulk % mass_in = bulk % mass_in - m_flux % n(s)
      end if
    end if
  end do
  call Comm_Mod_Global_Sum_Real(bulk % mass_in)

  !---------------------------------------!
  !   Calculate the outflow mass fluxes   !
  !     then correct it to satisfy the    !
  !          overall mass balance         !
  !---------------------------------------!
  bulk % mass_out = 0.0
  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2  < 0) then
      if(Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. OUTFLOW) then
        u % n(c2) = u % n(c1)
        v % n(c2) = v % n(c1)
        w % n(c2) = w % n(c1)
        m_flux % n(s) = flow % density_f(s) * ( u % n(c2)*grid % sx(s)    &
                                              + v % n(c2)*grid % sy(s)    &
                                              + w % n(c2)*grid % sz(s) )
        bulk % mass_out = bulk % mass_out + m_flux % n(s)
      end if

      if(Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. CONVECT  &
         .and. m_flux % n(s) > 0.0) then
        m_flux % n(s) = flow % density_f(s) * ( u % n(c2)*grid % sx(s)    &
                                              + v % n(c2)*grid % sy(s)    &
                                              + w % n(c2)*grid % sz(s) )
        bulk % mass_out = bulk % mass_out + m_flux % n(s)
      end if

      m_flux % n(s) = flow % density_f(s) * ( u % n(c2)*grid % sx(s)    &
                                            + v % n(c2)*grid % sy(s)    &
                                            + w % n(c2)*grid % sz(s) )

      if(Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. PRESSURE  &
         .and. m_flux % n(s) >0.0) then
        bulk % mass_out = bulk % mass_out + m_flux % n(s)
      end if

    end if
  end do
  call Comm_Mod_Global_Sum_Real(bulk % mass_out)  ! not checked

  fac = bulk % mass_in / (bulk % mass_out + TINY)

  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)
    if(c2  < 0) then
      if(Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. OUTFLOW  .or.  &
         Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. CONVECT  .or.  &
         Grid_Mod_Bnd_Cond_Type(grid, c2) .eq. PRESSURE) then
        u % n(c2) = u % n(c2) * fac
        v % n(c2) = v % n(c2) * fac
        w % n(c2) = w % n(c2) * fac
        m_flux % n(s) = m_flux % n(s) * fac
        bulk % mass_out = bulk % mass_out + m_flux % n(s)
      end if
    end if
  end do

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
c Copyright 2007-2015 Free Software Foundation, Inc.

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

      real*8 a,b
      complex*16 c

      a = 1000
      b = -50
      c = cmplx(a,b)
      write(*,*) s
      stop
      end

max_stars_repo_name: armoredsoftware/protocol

max_issues_repo_name: armoredsoftware/protocol

max_forks_repo_name: armoredsoftware/protocol


text:
module M {

  active component Commands {

    async input port p: P

    command reg port cmdRegOut
    command recv port cmdIn
    command resp port cmdResponseOut

    @ Command C1
    async command C1
    async command C2 priority 10
    async command C3 drop
    async command C4 priority 10 drop

    @ Command C5
    sync command C5(a: U32)

    @ Command C6
    guarded command C6(a: F32)

  }

}

ext: fpp

max_stars_repo_name: kevin-f-ortega/fpp

max_issues_repo_name: kevin-f-ortega/fpp

max_forks_repo_name: kevin-f-ortega/fpp


text: This Memorials memorial Address(Benches bench, 38.566, 121.749812) is located on the Covell Greenbelt near Covell Park and the intersection of Cabrillo and Del Oro Avenues.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      FUNCTION K91R(J,M,K,L,R, JA)
C BY JONATHAN PINKLEY
C MODIFIED FROM K15C BY JM 3/27/79
      DIMENSION IPOL(11,4), QC(4), QN(4), E(11)
      k91r=ja    ! Added 7/27/93 to report own old value
      IF (M .NE. 1) GO TO 30
C INITIAL BELIEFS
      X = .999
      PX = .001
      Y = .001
      PY = .999
      Z = .999
      PZ = .001
      W = .001
      PW = .999
      QC(1) = 1.999
      QC(2) =1.999
      QC(3) = 0.001
      QC(4) = 0.001
      DO 10 N = 1, 4
   10 QN(N) = 2
C DEFINE POLICIES(FIRST,WHAT IF OUTCOME=1)
      DATA IPOL /4*0, 7*1, 0, 3*1, 3*0, 4*1, 3*0, 1, 2*0, 1, 0, 0, 1, 1,
     1 2*0, 1, 0, 0, 1, 0, 0, 1, 0, 1/
       IOLD=0
      K91R = 0
      N = 0
      GO TO 100
C UPDATE STATS OF HIS CONTINGENCIES
C N IS OUTCOME OF M-2
   30 IF (M .LE. 2) GO TO 100
      IF (J .EQ. 0) QC(N) = QC(N) + 1
      QN(N) = QN(N) + 1
C REVERSE  Y AND Z
      GO TO (40, 60, 50, 70), N
   40 X = QC(1) / QN(1)
      PX = 1 - X
      GO TO 100
   50 Y = QC(3) / QN(3)
      PY = 1 - Y
      GO TO 100
   60 Z = QC(2) / QN(2)
      PZ = 1 - Z
      GO TO 100
   70 W = QC(4) / QN(4)
      GO TO 100
C CALC EXPECTATIONS OF POLICIES
  100 E(1) = (3*Z) / (Z + PX)
      E(2) = (3*(Y*Z + W*PZ) + 5*Z*PX + PX*PZ) / (Y*Z + W*PZ + PX + Z*
     1 PX + PX*PZ)
      E(3) = (3*W*Y + 5*W*PX + PX*PZ) / (W*Y + 2*W*PX + PX*PZ)
      E(4) = (3*W*PY + 5*Z*PX + PX*PY) / (W*PY + PX*PY + Z*PX + PX*PY)
      E(5) = (3*Z + 5*X*Z + Z*PX) / (1 - X*Y - W*PX + 2*Z)
      E(6) = (8*W*Z + Z*PX) / (2*W*Z + W*PY + Z*PX)
      E(7) = (3*Z*PY + 5*X*Z + Z*PY) / (2*Z*PY + PW*PY + X*Z)
      E(8) = (3*(Y*Z + W*PZ) + 5*(Z*PW + W*X) + 1 - X*Y - Z*PY) /(Y *Z +
     1 W*PZ + 2 - 2*X*Y - W*PX + Z*PW + W*X - Z*PY)
      E(9) = (3*W*Y + 5*W + 1 - X*Y - Z*PY) / (2*W + 1 - X*Y - Z*PY)
      E(10) = (3*W*PY + 5*(Z*PW + W*X) + PY) / (PY + Z*PW + W*X + PY)
      E(11) = (5*W + PY) / (W + PY)
C FIND POL WITH MAX E
      IBEST = 1
      BESTE = E(1)
      DO 80 I = 2, 11
        IF (E(I) .LE. BESTE) GO TO 80
        IBEST = I
        BESTE = E(I)
   80 CONTINUE
C CALC OUTCOME FOR USE IN CHOICE AND NEXT MV STATS
   90 N = 2 *  IOLD + J + 1
C CHOICE(CHOSEN POLICY,PREV OUTCOME)
      K91R = IPOL(IBEST,N)
       IOLD=K91R
      RETURN
      END FUNCTION

max_stars_repo_name: Axelrod-Python/TourExec

max_issues_repo_name: Axelrod-Python/TourExec

max_forks_repo_name: Axelrod-Python/TourExec


text:
      PROGRAM LOOPWHILE
        INTEGER I

C       FORTRAN 77 does not have a while loop, so we use GOTO statements
C       with conditions instead. This is one of two easy ways to do it.
        I = 1024
   10   CONTINUE
C       Check condition.
        IF (I .GT. 0) THEN
C         Handle I.
          WRITE (*,*) I
          I = I / 2
C         Jump back to before the IF block.
          GOTO 10
        ENDIF
        STOP
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 09.01.09.01/05                        *
C  *    TEST TITLE : Error 407 handled by user             *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      CALL INITGL ('09.01.09.01/05')
      CALL E09407 (1)
      CALL ENDIT
      END


text:
      SUBROUTINE CHEKSTAB(LUN)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    CHEKSTAB
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1994-01-06
C
C ABSTRACT: THIS SUBROUTINE CHECKS THAT AN INTERNAL BUFR TABLE
C   REPRESENTATION IS SELF-CONSISTENT AND FULLY DEFINED.  IF ANY ERRORS
C   ARE FOUND, THEN AN APPROPRIATE CALL IS MADE TO BUFR ARCHIVE LIBRARY
C   SUBROUTINE BORT.
C
C PROGRAM HISTORY LOG:
C 1994-01-06  J. WOOLLEN -- ORIGINAL AUTHOR
C 1995-06-28  J. WOOLLEN -- INCREASED THE SIZE OF INTERNAL BUFR TABLE
C                           ARRAYS IN ORDER TO HANDLE BIGGER FILES
C 1998-07-08  J. WOOLLEN -- REPLACED CALL TO CRAY LIBRARY ROUTINE
C                           "ABORT" WITH CALL TO NEW INTERNAL BUFRLIB
C                           ROUTINE "BORT"
C 1999-11-18  J. WOOLLEN -- THE NUMBER OF BUFR FILES WHICH CAN BE
C                           OPENED AT ONE TIME INCREASED FROM 10 TO 32
C                           (NECESSARY IN ORDER TO PROCESS MULTIPLE
C                           BUFR FILES UNDER THE MPI)
C 2003-11-04  J. ATOR    -- ADDED DOCUMENTATION
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED HISTORY
C                           DOCUMENTATION; OUTPUTS MORE COMPLETE
C                           DIAGNOSTIC INFO WHEN ROUTINE TERMINATES
C                           ABNORMALLY
C
C USAGE:    CALL CHEKSTAB (LUN)
C   INPUT ARGUMENT LIST:
C     LUN      - INTEGER: I/O STREAM INDEX INTO INTERNAL MEMORY ARRAYS
C
C REMARKS:
C    THIS ROUTINE CALLS:        BORT     NEMTAB   NEMTBB   NEMTBD
C    THIS ROUTINE IS CALLED BY: MAKESTAB
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /TABABD/ NTBA(0:NFILES),NTBB(0:NFILES),NTBD(0:NFILES),
     .                MTAB(MAXTBA,NFILES),IDNA(MAXTBA,NFILES,2),
     .                IDNB(MAXTBB,NFILES),IDND(MAXTBD,NFILES),
     .                TABA(MAXTBA,NFILES),TABB(MAXTBB,NFILES),
     .                TABD(MAXTBD,NFILES)

      CHARACTER*600 TABD
      CHARACTER*128 TABB
      CHARACTER*128 TABA
      CHARACTER*128 BORT_STR
      CHARACTER*24  UNIT
      CHARACTER*8   NEMO,NEMS(MAXCD)
      CHARACTER*1   TAB
      DIMENSION     IRPS(MAXCD),KNTS(MAXCD)

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

C  THERE MUST BE ENTRIES IN TABLES A, B, AND D
C  -------------------------------------------

      IF(NTBA(LUN).EQ.0) GOTO 900
      IF(NTBB(LUN).EQ.0) GOTO 901
      IF(NTBD(LUN).EQ.0) GOTO 902

C  MAKE SURE EACH TABLE A ENTRY DEFINED AS A SEQUENCE
C  --------------------------------------------------

      DO I=1,NTBA(LUN)
      NEMO = TABA(I,LUN)(4:11)
      CALL NEMTAB(LUN,NEMO,IDN,TAB,IRET)
      IF(TAB.NE.'D') GOTO 903
      ENDDO

C  CHECK TABLE B CONTENTS
C  ----------------------

      DO ITAB=1,NTBB(LUN)
      CALL NEMTBB(LUN,ITAB,UNIT,ISCL,IREF,IBIT)
      ENDDO

C  CHECK TABLE D CONTNETS
C  ----------------------

      DO ITAB=1,NTBD(LUN)
      CALL NEMTBD(LUN,ITAB,NSEQ,NEMS,IRPS,KNTS)
      ENDDO

C  EXITS
C  -----

      RETURN
900   CALL BORT
     . ('BUFRLIB: CHEKSTAB - EMPTY TABLE A IN INTERNAL BUFR TABLES')
901   CALL BORT
     . ('BUFRLIB: CHEKSTAB - EMPTY TABLE B IN INTERNAL BUFR TABLES')
902   CALL BORT
     . ('BUFRLIB: CHEKSTAB - EMPTY TABLE D IN INTERNAL BUFR TABLES')
903   WRITE(BORT_STR,'("BUFRLIB: CHEKSTAB - TABLE A ENTRY: ",A," NOT '//
     . 'DEFINED AS A SEQUENCE")') NEMO
      CALL BORT(BORT_STR)
      END


text:
!-----------------------------------------------------------------------
!
! Header with common subroutine
!
      include 't320-basis-f.h'
!-----------------------------------------------------------------------
      subroutine feval(x1,x2,val)
      real*8 x1,x2,val

      val=x1*x1+x2*x2+x1*x2+1

      end
!-----------------------------------------------------------------------
      subroutine dfeval(x1,x2,val)
      real*8 x1,x2,val

      val=2*x1+x2

      end
!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceed/fortran.h'

      integer ceed,err
      integer input,output
      integer p,q,d
      parameter(p=6)
      parameter(q=4)
      parameter(d=2)

      real*8 qref(d*q)
      real*8 qweight(q)
      real*8 interp(p*q)
      real*8 grad(d*p*q)
      real*8 xq(d*q)
      real*8 xr(d*p)
      real*8 iinput(p)
      real*8 ooutput(d*q)
      real*8 val,diff
      real*8 x1,x2
      integer*8 ioffset,ooffset

      integer b,i

      character arg*32

      xq=(/2.d-1,6.d-1,1.d0/3.d0,2.d-1,2.d-1,2.d-1,   1.d0/3.d0,6.d-1/)
      xr=(/0.d0,5.d-1,1.d0,0.d0,5.d-1,0.d0,0.d0,0.d0,   0.d0,5.d-1,5.d-1,1.d0/)

      call getarg(1,arg)

      call buildmats(qref,qweight,interp,grad)

      call ceedinit(trim(arg)//char(0),ceed,err)

      call ceedbasiscreateh1(ceed,ceed_triangle,1,p,q,interp,grad,qref,qweight,&
     & b,err)

      do i=1,p
        x1=xr(0*p+i)
        x2=xr(1*p+i)
        call feval(x1,x2,val)
        iinput(i)=val
      enddo

      call ceedvectorcreate(ceed,p,input,err)
      ioffset=0
      call ceedvectorsetarray(input,ceed_mem_host,ceed_use_pointer,iinput,&
     & ioffset,err)
      call ceedvectorcreate(ceed,q*d,output,err)
      call ceedvectorsetvalue(output,0.d0,err)

      call ceedbasisapply(b,1,ceed_notranspose,ceed_eval_grad,input,output,err)

      call ceedvectorgetarrayread(output,ceed_mem_host,ooutput,ooffset,err)
      do i=1,q
        x1=xq(0*q+i)
        x2=xq(1*q+i)
        call dfeval(x1,x2,val)
        diff=val-ooutput(0*q+i+ooffset)
        if (abs(diff)>1.0d-10) then
! LCOV_EXCL_START
          write(*,'(A,I1,A,F12.8,A,F12.8)')  '[',i,'] ',ooutput(i+ooffset),&
     &     ' != ',val
! LCOV_EXCL_STOP
        endif
        call dfeval(x2,x1,val)
        diff=val-ooutput(1*q+i+ooffset)
        if (abs(diff)>1.0d-10) then
! LCOV_EXCL_START
          write(*,'(A,I1,A,F12.8,A,F12.8)')  '[',i,'] ',ooutput(i+ooffset),&
     &     ' != ',val
! LCOV_EXCL_STOP
        endif
      enddo
      call ceedvectorrestorearrayread(output,ooutput,ooffset,err)

      call ceedvectordestroy(input,err)
      call ceedvectordestroy(output,err)
      call ceedbasisdestroy(b,err)
      call ceeddestroy(ceed,err)

      end
!-----------------------------------------------------------------------

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
! PWR008: Declare the intent for each procedure parameter
! https://www.appentra.com/knowledge/checks/pwr008

SUBROUTINE example(a, b)
    IMPLICIT NONE
    INTEGER :: a
    INTEGER :: b

    a = 5
    b = a * 2 ! This is allowed because no intent has been declared
END SUBROUTINE example

max_stars_repo_name: teamappentra/check-examples

max_issues_repo_name: teamappentra/check-examples

max_forks_repo_name: teamappentra/check-examples


text:
C Copyright (c) 2007 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
C retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

      SUBROUTINE STRAN(ISRCHR,NISR,SOLEA,SOLEB,
     &                 IDBLKA,IDBLKB,
     &                 ITT,iblk,TIMES,CENTER,
     &                 INSUB,ICOMPL,
     &                 XB,YB,ZB,ICONB,DUME)
C
      include 'aexds1.blk'
      include 'aexds2.blk'
      include 'bmesh.blk'
      include 'ebbyeb.blk'
      include 'ex2tp.blk'
      include 'steps.blk'
      include 'tapes.blk'
C
      DIMENSION SOLEB(NUMEBB,*),SOLEA(NUMEBA,*),ISRCHR(NISR,*),
     &          ITT(NVAREL,*),TIMES(*),CENTER(NUMEBB,*)
      DIMENSION XB(*),YB(*),ZB(*),XX(27),YY(27),ZZ(27),ICONB(NELNDB,*)
      DIMENSION DUME(*)
C
C
      IF (ISTEP .EQ. -1)THEN
        NTM = NTIMES
      ELSE
        NTM = 1
      END IF
C
      DO 5 IST = 1, NTM
        IF (ISTEP .EQ. -1)THEN
          ISTP = IST
        ELSE
          ISTP = ISTEP
        END IF
C
        DO 10 IVAR = 1, NVAREL
          IF (ITT(IVAR,iblk) .EQ. 0)GO TO 10
C
C If first time into subroutine for this element block,
C initialize the SOLEB array
C If not first time into subroutine for this element block
C retrieve SOLEB from storage in EXODUS
C
          IF (INSUB .EQ. 1)THEN
            CALL INIELT(SOLEB,IVAR,TIMES,ISTP,IDBLK,CENTER,DUME)
          ELSE
            CALL EXGEV(NTP4EX,IST,IVAR,IDBLKB,NUMEBB,SOLEB(1,IVAR),
     &                 IERR)
          END IF
C
          CALL EXGEV(NTP2EX,ISTP,IVAR,IDBLKA,NUMEBA,SOLEA(1,IVAR),IERR)
          DO 20 IELB = 1, NUMEBB
            IELA = ISRCHR(1,IELB)
            IF (IELA .NE. 0)THEN
              SOLEB(IELB,IVAR) = SOLEA(IELA,IVAR)
            END IF
   20     CONTINUE
C
C If there is more searching to do (i.e. many blocks to one)
C use EXODUS as temporary storage
C don't bother to perform needed adjustments yet
C
          IF (ICOMPL .NE. 1)THEN
            CALL EXPEV(NTP4EX,IST,IVAR,IDBLKB,NUMEBB,SOLEB(1,IVAR),
     &                 IERR)
          ELSE
C
C Make needed adjustments to element varaible data and
C write element vars out to EXODUS data base
C
C ELMASS is special
C
            IF (NAMVAR(NVARGP+IVAR)(1:6) .EQ. 'ELMASS')THEN
C
C ELMASS was changed to nodal density prior to processing.
C need to go back from density to element mass now
C
              NNODES = 4
              DO 100 IEL = 1, NUMEBB
                DO 105 I = 1, NNODES
                  XX(I) = XB(ICONB(I,IEL))
                  YY(I) = YB(ICONB(I,IEL))
                  ZZ(I) = ZB(ICONB(I,IEL))
 105            CONTINUE
                CALL VOL(ITYPE,XX,YY,ZZ,VOLUME)
                SOLEB(IEL,IVAR) = SOLEB(IEL,IVAR) * VOLUME
 100          CONTINUE
            END IF
            CALL EXPEV(NTP4EX,IST,IVAR,IDBLKB,NUMEBB,SOLEB(1,IVAR),
     &                 IERR)
          END IF
C
   10   CONTINUE
    5 CONTINUE
      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
      SUBROUTINE vegas(region,ndim,fxn,init,ncall,itmx,nprn,tgral,sd,
     *chi2a)
      INTEGER init,itmx,ncall,ndim,nprn,NDMX,MXDIM
      REAL tgral,chi2a,sd,region(2*ndim),fxn,ALPH,TINY
      PARAMETER (ALPH=1.5,NDMX=50,MXDIM=10,TINY=1.e-30)
      EXTERNAL fxn
CU    USES fxn,ran2,rebin
      INTEGER i,idum,it,j,k,mds,nd,ndo,ng,npg,ia(MXDIM),kg(MXDIM)
      REAL calls,dv2g,dxg,f,f2,f2b,fb,rc,ti,tsi,wgt,xjac,xn,xnd,xo,
     *d(NDMX,MXDIM),di(NDMX,MXDIM),dt(MXDIM),dx(MXDIM),r(NDMX),x(MXDIM),
     *xi(NDMX,MXDIM),xin(NDMX),ran2
      DOUBLE PRECISION schi,si,swgt
      COMMON /ranno/ idum
      SAVE
      if(init.le.0)then
        mds=1
        ndo=1
        do 11 j=1,ndim
          xi(1,j)=1.
11      continue
      endif
      if (init.le.1)then
        si=0.
        swgt=0.
        schi=0.
      endif
      if (init.le.2)then
        nd=NDMX
        ng=1
        if(mds.ne.0)then
          ng=(ncall/2.+0.25)**(1./ndim)
          mds=1
          if((2*ng-NDMX).ge.0)then
            mds=-1
            npg=ng/NDMX+1
            nd=ng/npg
            ng=npg*nd
          endif
        endif
        k=ng**ndim
        npg=max(ncall/k,2)
        calls=float(npg)*float(k)
        dxg=1./ng
        dv2g=(calls*dxg**ndim)**2/npg/npg/(npg-1.)
        xnd=nd
        dxg=dxg*xnd
        xjac=1./calls
        do 12 j=1,ndim
          dx(j)=region(j+ndim)-region(j)
          xjac=xjac*dx(j)
12      continue
        if(nd.ne.ndo)then
          do 13 i=1,max(nd,ndo)
            r(i)=1.
13        continue
          do 14 j=1,ndim
            call rebin(ndo/xnd,nd,r,xin,xi(1,j))
14        continue
          ndo=nd
        endif
        if(nprn.ge.0) write(*,200) ndim,calls,it,itmx,nprn,ALPH,mds,nd,
     *(j,region(j),j,region(j+ndim),j=1,ndim)
      endif
      do 28 it=1,itmx
        ti=0.
        tsi=0.
        do 16 j=1,ndim
          kg(j)=1
          do 15 i=1,nd
            d(i,j)=0.
            di(i,j)=0.
15        continue
16      continue
10      continue
          fb=0.
          f2b=0.
          do 19 k=1,npg
            wgt=xjac
            do 17 j=1,ndim
              xn=(kg(j)-ran2(idum))*dxg+1.
              ia(j)=max(min(int(xn),NDMX),1)
              if(ia(j).gt.1)then
                xo=xi(ia(j),j)-xi(ia(j)-1,j)
                rc=xi(ia(j)-1,j)+(xn-ia(j))*xo
              else
                xo=xi(ia(j),j)
                rc=(xn-ia(j))*xo
              endif
              x(j)=region(j)+rc*dx(j)
              wgt=wgt*xo*xnd
17          continue
            f=wgt*fxn(x,wgt)
            f2=f*f
            fb=fb+f
            f2b=f2b+f2
            do 18 j=1,ndim
              di(ia(j),j)=di(ia(j),j)+f
              if(mds.ge.0) d(ia(j),j)=d(ia(j),j)+f2
18          continue
19        continue
          f2b=sqrt(f2b*npg)
          f2b=(f2b-fb)*(f2b+fb)
          if (f2b.le.0.) f2b=TINY
          ti=ti+fb
          tsi=tsi+f2b
          if(mds.lt.0)then
            do 21 j=1,ndim
              d(ia(j),j)=d(ia(j),j)+f2b
21          continue
          endif
        do 22 k=ndim,1,-1
          kg(k)=mod(kg(k),ng)+1
          if(kg(k).ne.1) goto 10
22      continue
        tsi=tsi*dv2g
        wgt=1./tsi
        si=si+dble(wgt)*dble(ti)
        schi=schi+dble(wgt)*dble(ti)**2
        swgt=swgt+dble(wgt)
        tgral=si/swgt
        chi2a=max((schi-si*tgral)/(it-.99d0),0.d0)
        sd=sqrt(1./swgt)
        tsi=sqrt(tsi)
        if(nprn.ge.0)then
          write(*,201) it,ti,tsi,tgral,sd,chi2a
          if(nprn.ne.0)then
            do 23 j=1,ndim
              write(*,202) j,(xi(i,j),di(i,j),i=1+nprn/2,nd,nprn)
23          continue
          endif
        endif
        do 25 j=1,ndim
          xo=d(1,j)
          xn=d(2,j)
          d(1,j)=(xo+xn)/2.
          dt(j)=d(1,j)
          do 24 i=2,nd-1
            rc=xo+xn
            xo=xn
            xn=d(i+1,j)
            d(i,j)=(rc+xn)/3.
            dt(j)=dt(j)+d(i,j)
24        continue
          d(nd,j)=(xo+xn)/2.
          dt(j)=dt(j)+d(nd,j)
25      continue
        do 27 j=1,ndim
          rc=0.
          do 26 i=1,nd
            if(d(i,j).lt.TINY) d(i,j)=TINY
            r(i)=((1.-d(i,j)/dt(j))/(log(dt(j))-log(d(i,j))))**ALPH
            rc=rc+r(i)
26        continue
          call rebin(rc/xnd,nd,r,xin,xi(1,j))
27      continue
28    continue
      return
200   FORMAT(/' input parameters for vegas:  ndim=',i3,'  ncall=',
     *f8.0/28x,'  it=',i5,'  itmx=',i5/28x,'  nprn=',i3,'  alph=',
     *f5.2/28x,'  mds=',i3,'   nd=',i4/(30x,'xl(',i2,')= ',g11.4,' xu(',
     *i2,')= ',g11.4))
201   FORMAT(/' iteration no.',I3,': ','integral =',g14.7,'+/- ',g9.2/
     *' all iterations:   integral =',g14.7,'+/- ',g9.2,
     *' chi**2/it''n =',g9.2)
202   FORMAT(/' data for axis ',I2/'    X       delta i       ',
     *'   x       delta i       ','    x       delta i       ',/(1x,
     *f7.5,1x,g11.4,5x,f7.5,1x,g11.4,5x,f7.5,1x,g11.4))
      END

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/vegas.for