text:
      	subroutine centered_diff_fort(u,D,dx,N,Nsq)
  		integer :: N
  		integer, intent(in) :: Nsq
	  	double precision, intent(in) :: u(Nsq), dx
	  	double precision, intent(out) :: D(Nsq)
	  	integer :: i,j
	  	do i = 2,N-1
	      		do j = 2,N-1
	          		D((i-1)*N+j) = (u((i-1)*N+j+1) + u((i-2)*N+j)
     &			+ u((i-1)*N+j-1) + u(i*N+j) - 4*u((i-1)*N+j))/dx**2
	      		end do
	  	end do
C --- i=1,i=N
		do j=2,N-1
			D(j) = (u(j+1) + u((N-2)*N+j) + u(j-1) + u(N+j) - 4*u(j))/dx**2
			D((N-1)*N+j) = D(j)
		end do
C --- j=1,j=N
		do i=2,N-1
			D((i-1)*N+1) = (u((i-2)*N+1) + u((i-1)*N+2) + u((i-1)*N+N-1)
     &			 + u(i*N+1) - 4*u((i-1)*N+1))/dx**2
			D((i-1)*N+N) = D((i-1)*N+1)
		end do
C --- edges
		D(1) = (u(2) + u(N+1) + u(N-1) + u((N-2)*N+1) - 4*u(1))/dx**2
		D(N) = D(1)
		D((N-1)*N+1) = D(1)
		D(N*N) = D(1)
      	end subroutine centered_diff_fort

        SUBROUTINE fnbruss(Y,t,F,N,Nsq)
C --- RIGHT-HAND SIDE OF N-body problem
		integer, intent(in) :: Nsq, N
		double precision, intent(in) :: t, Y(2*Nsq)
		double precision, intent(out) :: F(2*Nsq)
		double precision :: Du(Nsq), Dv(Nsq), U(Nsq), V(Nsq)
		double precision :: uval, vval, finhom, dx
		alpha=0.1
		dx = 1.D0/(N-1)
		U = Y(1:Nsq)
		V = Y((Nsq+1):2*Nsq)
		CALL centered_diff_fort(U,Du,dx,N, Nsq)
		CALL centered_diff_fort(V,Dv,dx,N, Nsq)
		do i=1,Nsq
			uval = U(i)
			vval = V(i)
			finhom = BRUSS2DInhom(i,t,N,dx)
			F(i) = 1+(uval**2)*vval-4.4*uval+alpha*Du(i)+finhom
			F(Nsq+i) = 3.4*uval-(uval**2)*vval+alpha*Dv(i)
        	end do
        end subroutine fnbruss

	function BRUSS2DInhom(i,t,N,dx)
		integer, intent(in) :: i, N
		integer :: row
		double precision, intent(in) :: t, dx
		BRUSS2DInhom = 0.D0
		row = (i-1)/N
		y = dx*row
		x = dx*(i-1-row*N)
   		if((x-0.3)**2+(y-0.6)**2<=0.01 .AND. t>=1.1) then
	   		BRUSS2DInhom = 5.D0
		end if

	end function BRUSS2DInhom

ext: f

max_stars_repo_name: numerical-mathematics/extrapolation

max_issues_repo_name: numerical-mathematics/extrapolation

max_forks_repo_name: numerical-mathematics/extrapolation


text: iszero (pred (succ (succ 0)));

max_stars_repo_name: brunoflores/arith

max_issues_repo_name: brunoflores/arith

max_forks_repo_name: brunoflores/arith


text:
C
C $Id: mapsav.f,v 1.19 2008-09-18 00:42:16 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MAPSAV (IFNO)
        INTEGER IFNO
        IF (ICFELL('MAPSAV - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
        CALL MDPSAV (IFNO)
        IF (ICFELL('MAPSAV',2).NE.0) RETURN
        RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/work/rep/PARIO/src/parutilio.f,v 1.2 2012/04/11 13:47:45 sjr Exp $

C what(1) key, module and SID; SCCS file; date and time of last delta:
C %W% %P% %G% %U%

C....................................................................
C  CONTAINS:  Interface for PAR I/O procedures
C
C  REVISION HISTORY:
C       Original version 02/11 by Shawn Roselle
C       Modified 12/09/2015 by David Wong
C          -- Added an optional argument in PIO_RE_INIT routine
C       Modified 01/26/2016 by Jeff Young
C          -- Change WFLG to an optional argument, eliminate PIO_RE_INIT, and make
C          PIO_INIT the only routine with a new optional IO_PE_INCLUSIVE argument;
C          replace WRSUBMAP with the new WRSUBDMAP
C       Modified 02/01/2019 by David Wong
C          removed unnecessary function declaration
C....................................................................

      MODULE PARUTILIO

         IMPLICIT NONE

C...Declare PWRITE3 outside the interface block:  with check interface
C...  option, compilers complained about type mismatches for buffer
C...  array passed to function PWRITE3 (when included in interface).

         LOGICAL, EXTERNAL :: PWRITE3

C...Interface

         INTERFACE

!           SUBROUTINE BOUNDARY ( GNBNDY, NBNDY, NLAYS, BTHICK, NEIGHBOR,
!    &                            MY_COL1, MY_ROW1, NCOL, NROW,
!    &                            GNCOLS, GNROWS, GARRAY, LARRAY )
!              INTEGER :: NLAYS
!              INTEGER :: NBNDY
!              INTEGER :: GNBNDY
!              INTEGER :: BTHICK
!              INTEGER :: NEIGHBOR(8)
!              INTEGER :: MY_COL1
!              INTEGER :: MY_ROW1
!              INTEGER :: NCOL
!              INTEGER :: NROW
!              INTEGER :: GNCOLS
!              INTEGER :: GNROWS
!              REAL    :: GARRAY( GNBNDY,NLAYS )
!              REAL    :: LARRAY( NBNDY,NLAYS )
!           END SUBROUTINE BOUNDARY

            SUBROUTINE GET_WRITE_MAP ( NP, NPC, NPR, NCOLS3D, NROWS3D, NLAYS3D )
               INTEGER :: NP
               INTEGER :: NPC
               INTEGER :: NPR
               INTEGER :: NCOLS3D
               INTEGER :: NROWS3D
               INTEGER :: NLAYS3D
            END SUBROUTINE GET_WRITE_MAP

!           LOGICAL FUNCTION GROWBUF( PTR_ARRAY, END )
!              USE PINTERPB_MODULE
!              TYPE (MEM_TYPE) :: PTR_ARRAY
!              INTEGER, INTENT(IN) :: END
!           END FUNCTION GROWBUF

            LOGICAL FUNCTION GTNDXHDV( FILE, VAR, JDATE, JTIME, VARSIZE, NBVS,
     &                                 ENDBUF, VX, NEWVAR )
               CHARACTER(16) :: FILE
               CHARACTER(16) :: VAR
               INTEGER :: JDATE
               INTEGER :: JTIME
               INTEGER :: VARSIZE
               INTEGER :: NBVS
               INTEGER :: ENDBUF
               INTEGER :: VX
               LOGICAL :: NEWVAR
            END FUNCTION GTNDXHDV

            LOGICAL FUNCTION INTERPOL( DATE, TIME, DATE1, TIME1, DATE2, TIME2,
     &                                 NUMVALS, VALSIN1, VALSIN2, VALSOUT, S_IND, E_IND )
               INTEGER, INTENT(IN) :: NUMVALS
               INTEGER, INTENT(IN) :: DATE
               INTEGER, INTENT(IN) :: TIME
               INTEGER, INTENT(IN) :: DATE1
               INTEGER, INTENT(IN) :: TIME1
               INTEGER, INTENT(IN) :: DATE2
               INTEGER, INTENT(IN) :: TIME2
               REAL, INTENT(IN)    :: VALSIN1( NUMVALS )
               REAL, INTENT(IN)    :: VALSIN2( NUMVALS )
               REAL, INTENT(OUT)   :: VALSOUT( NUMVALS )
               INTEGER, INTENT(IN), OPTIONAL :: S_IND, E_IND
            END FUNCTION INTERPOL

!           LOGICAL FUNCTION PINTERPB( FILNAME, VARNAME, CALLER, JDATE, JTIME,
!    &                                  VSIZE, VARRAY, LVL )
!              INTEGER, INTENT(IN)       :: VSIZE
!              CHARACTER(16), INTENT(IN) :: FILNAME
!              CHARACTER(*), INTENT(IN)  :: VARNAME
!              CHARACTER(*), INTENT(IN)  :: CALLER
!              INTEGER, INTENT(IN)       :: JDATE
!              INTEGER, INTENT(IN)       :: JTIME
!              REAL, INTENT(OUT)         :: VARRAY( VSIZE )
!              INTEGER, INTENT(IN), OPTIONAL :: LVL
!           END FUNCTION PINTERPB

            LOGICAL FUNCTION PIO_INIT( colrow, gl_ncols, gl_nrows, nlays,
     &                                 nthik, ncols, nrows, npcol, nprow,
     &                                 nprocs, mype, wflg, io_pe_inclusive )
               CHARACTER(2), INTENT(INOUT) :: colrow
               INTEGER, INTENT(IN) :: gl_ncols
               INTEGER, INTENT(IN) :: gl_nrows
               INTEGER, INTENT(IN) :: nlays
               INTEGER, INTENT(IN) :: nthik
               INTEGER, INTENT(IN) :: ncols
               INTEGER, INTENT(IN) :: nrows
               INTEGER, INTENT(IN) :: npcol
               INTEGER, INTENT(IN) :: nprow
               INTEGER, INTENT(IN) :: nprocs
               INTEGER, INTENT(IN) :: mype
               LOGICAL, INTENT(IN), OPTIONAL :: wflg
               LOGICAL, INTENT(IN), OPTIONAL :: io_pe_inclusive
            END FUNCTION PIO_INIT
C Note: more than one optional argument requires that the caller calls by name

            SUBROUTINE PM3ERR ( CALLER, JDATE, JTIME, ERRTXT, FATAL )
               CHARACTER(*) :: CALLER
               INTEGER      :: JDATE
               INTEGER      :: JTIME
               CHARACTER(*) :: ERRTXT
               LOGICAL      :: FATAL
            END SUBROUTINE PM3ERR

            SUBROUTINE PM3EXIT ( CALLER, JDATE, JTIME, MSGTXT, EXITSTAT )
               CHARACTER(*) :: CALLER
               INTEGER      :: JDATE
               INTEGER      :: JTIME
               CHARACTER(*) :: MSGTXT
               INTEGER      :: EXITSTAT
            END SUBROUTINE PM3EXIT

            SUBROUTINE PM3WARN ( CALLER, JDATE, JTIME, ERRTXT )
               CHARACTER(*) :: CALLER
               INTEGER      :: JDATE
               INTEGER      :: JTIME
               CHARACTER(*) :: ERRTXT
            END SUBROUTINE PM3WARN

            LOGICAL FUNCTION PSHUT3 ()
            END FUNCTION PSHUT3

            LOGICAL FUNCTION PTRWRITE3( FNAME, VNAME, JDATE, JTIME, BUFFER )
               use HDMod
               CHARACTER(*)  :: FNAME
               CHARACTER(*)  :: VNAME
               INTEGER       :: JDATE
               INTEGER       :: JTIME
               TYPE(hyperdual), POINTER :: BUFFER( :,:,:,: )
            END FUNCTION PTRWRITE3

            LOGICAL FUNCTION PWRGRDD( FILNAME, VARNAME, DATE, TIME, BUFFER,
     &                                NCOLS3D, NROWS3D, NLAYS3D, NCOLS, NROWS,
     &                                NP )
               INTEGER      :: NROWS
               INTEGER      :: NCOLS
               INTEGER      :: NLAYS3D
               CHARACTER(*) :: FILNAME
               CHARACTER(*) :: VARNAME
               INTEGER      :: DATE
               INTEGER      :: TIME
               REAL         :: BUFFER( NCOLS,NROWS,NLAYS3D )
               INTEGER      :: NCOLS3D
               INTEGER      :: NROWS3D
               INTEGER      :: NP
            END FUNCTION PWRGRDD

            LOGICAL FUNCTION READBNDY( FILNAME, VARNAME, VX, NBNDY, NLAYS,
     &                                 DATE, TIME, RFLAG, FLIP )
               CHARACTER(16) :: FILNAME
               CHARACTER(16) :: VARNAME
               INTEGER :: VX
               INTEGER :: NBNDY
               INTEGER :: NLAYS
               INTEGER :: DATE(2)
               INTEGER :: TIME(2)
               INTEGER :: RFLAG
               INTEGER :: FLIP
            END FUNCTION READBNDY

            SUBROUTINE SUBDMAP ( NUMPROCS, NCOLS, NROWS, NLAYS, NPCOL, NPROW,
     &                           NCOLS_PE, NROWS_PE, COLSX_PE, ROWSX_PE,
     &                           MAXCELLS, IERR )
               INTEGER :: NUMPROCS
               INTEGER :: NCOLS
               INTEGER :: NROWS
               INTEGER :: NLAYS
               INTEGER :: NPCOL
               INTEGER :: NPROW
               INTEGER :: NCOLS_PE( NUMPROCS )
               INTEGER :: NROWS_PE( NUMPROCS )
               INTEGER :: COLSX_PE( 2,NUMPROCS )
               INTEGER :: ROWSX_PE( 2,NUMPROCS )
               INTEGER :: MAXCELLS
               INTEGER :: IERR
            END SUBROUTINE SUBDMAP

            SUBROUTINE WRSUBDMAP ( NUMPROCS, NCOLS_PE, NROWS_PE, COLSX_PE,
     &                             ROWSX_PE )
               INTEGER :: NUMPROCS
               INTEGER :: NCOLS_PE(*)
               INTEGER :: NROWS_PE(*)
               INTEGER :: COLSX_PE(2,*)
               INTEGER :: ROWSX_PE(2,*)
            END SUBROUTINE WRSUBDMAP

         END INTERFACE

      END MODULE PARUTILIO

max_stars_repo_name: echen-git/echenhydCMAQ

max_issues_repo_name: echen-git/echenhydCMAQ

max_forks_repo_name: echen-git/echenhydCMAQ


! { dg-do run }
!
! PR 41829: [OOP] Runtime error with dynamic dispatching.  Tests
! dynamic dispatch in a case where the caller knows nothing about
! the dynamic type at compile time.
!
! Contributed by Salvatore Filippone <sfilippone@uniroma2.it>
!
module foo_mod
  type foo
    integer :: i
  contains
    procedure, pass(a) :: doit
    procedure, pass(a) :: getit
  end type foo

  private doit,getit
contains
  subroutine  doit(a)
    class(foo) :: a

    a%i = 1
!    write(*,*) 'FOO%DOIT base version'
  end subroutine doit
  function getit(a) result(res)
    class(foo) :: a
    integer :: res

    res = a%i
  end function getit

end module foo_mod
module foo2_mod
  use foo_mod

  type, extends(foo) :: foo2
    integer :: j
  contains
    procedure, pass(a) :: doit  => doit2
    procedure, pass(a) :: getit => getit2
  end type foo2

  private doit2, getit2

contains

  subroutine  doit2(a)
    class(foo2) :: a

    a%i = 2
    a%j = 3
!    write(*,*) 'FOO2%DOIT derived version'
  end subroutine doit2
  function getit2(a) result(res)
    class(foo2) :: a
    integer :: res

    res = a%j
  end function getit2

end module foo2_mod

module bar_mod
  use foo_mod
  type bar
    class(foo), allocatable :: a
  contains
    procedure, pass(a) :: doit
    procedure, pass(a) :: getit
  end type bar
  private doit,getit

contains
  subroutine doit(a)
    class(bar) :: a

    call a%a%doit()
  end subroutine doit
  function getit(a) result(res)
    class(bar) :: a
    integer :: res

    res = a%a%getit()
  end function getit
end module bar_mod

program testd10
  use foo_mod
  use foo2_mod
  use bar_mod

  type(bar) :: a

  allocate(foo :: a%a)
  call a%doit()
!  write(*,*) 'Getit value : ', a%getit()
  if (a%getit() .ne. 1) call abort
  deallocate(a%a)
  allocate(foo2 :: a%a)
  call a%doit()
!  write(*,*) 'Getit value : ', a%getit()
  if (a%getit() .ne. 3) call abort

end program testd10


text:
      SUBROUTINE TPERI(SEMI,Q,UPR,MB,DT)
*
*
*       Pericentre time for KS motion.
*       ------------------------------
*
*       Routine for calculating the pericentre time
*       from the KS coordinates Q and derivatives UPR.
*       SEMI = semi-major axis of dominant motion.
*       Q = KS coordinates of relative position vector.
*       UPR = derivatives of Q in standard KS formulation.
*       MB = combined mass of the two bodies.
*       DT = pericentre time interval (DT < 0 before peri).
*       Algorithm due to Seppo Mikkola (1991).
*
      IMPLICIT REAL*8  (A-H,M,O-Z)
      REAL*8  Q(4),UPR(4)
*
*
*       Set scalar KS distance, velocity and radial velocity.
      R = Q(1)**2 + Q(2)**2 + Q(3)**2 + Q(4)**2
      V2 = UPR(1)**2 + UPR(2)**2 + UPR(3)**2 + UPR(4)**2
      VR = 2.0D0*(Q(1)*UPR(1) + Q(2)*UPR(2) + Q(3)*UPR(3) + Q(4)*UPR(4))
*
*       Define auxiliary variables.
      PSI = VR/SQRT(MB)
      ZETA = 4.0D0*V2/MB - 1.0D0
*
*       Form semi-major axis & eccentricity.
*     ALPHA = 2.0D0*(1.0D0 - 2.0D0*V2/MB)/R
*     SEMI = 1.0D0/ALPHA
*
*       Employ regular value of semi-major axis determined by routine EREL.
      ALPHA = 1.0D0/SEMI
      ECC = SQRT(ZETA*ZETA + ALPHA*PSI*PSI)
*
*       Distinguish between near-collision, elliptic & hyperbolic case.
      Q1 = (PSI/ECC)**2
      IF (ZETA.GT.0.0D0.AND.ABS(Q1/SEMI).LT.0.1) THEN
*       Expansion of THETA/(SIN(THETA)*COS(THETA)) - 1 for Q = SIN(THETA)**2.
          SN1 = 1.0
          S = 0.0D0
          Q1 = Q1/SEMI
          DO 1 IT = 1,10
              SN1 = Q1*SN1*FLOAT(2*IT)/FLOAT(2*IT + 1)
              S = S + SN1
    1     CONTINUE
          S = S + SN1*Q1*0.9D0/(1.0 - Q1)
          DT = (R*ZETA - PSI**2 + ZETA*S*SEMI)*PSI/(ECC**2*SQRT(MB))
      ELSE IF (SEMI.GT.0.0D0) THEN
*       Determine the eccentric anomaly with respect to pericentre (-PI,PI).
          THETA = DATAN2(PSI/SQRT(SEMI),ZETA)
*       Obtain pericentre time interval from Kepler's equation.
          DT = SEMI*SQRT(SEMI/MB)*(THETA - PSI/SQRT(SEMI))
      ELSE IF (SEMI.LT.0.0D0) THEN
*       Hyperbolic case.
          A1 = PSI/(ECC*SQRT(ABS(SEMI)))
*       Use EXP(F) = SINH(F) + COSH(F) to obtain the eccentric anomaly THETA.
          A2 = ABS(A1) + SQRT(A1**2 + 1.0D0)
          THETA = LOG(A2)
          IF (A1.LT.0.0D0) THETA = -THETA
          A0 = ABS(SEMI)
          DT = A0*SQRT(A0/MB)*(PSI/SQRT(A0) - THETA)
      END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
c---------------------------------------------------------------------
c---------------------------------------------------------------------

       subroutine  tzetar

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   block-diagonal matrix-vector multiplication
c---------------------------------------------------------------------

       include 'header.h'

       integer i, j, k
       double precision  t1, t2, t3, ac, xvel, yvel, zvel, r1, r2, r3,
     >                   r4, r5, btuz, ac2u, uzik1

       if (timeron) call timer_start(t_tzetar)
       do    k = 1, nz2
          do    j = 1, ny2
             do    i = 1, nx2

                xvel = us(i,j,k)
                yvel = vs(i,j,k)
                zvel = ws(i,j,k)
                ac   = speed(i,j,k)

                ac2u = ac*ac

                r1 = rhs(1,i,j,k)
                r2 = rhs(2,i,j,k)
                r3 = rhs(3,i,j,k)
                r4 = rhs(4,i,j,k)
                r5 = rhs(5,i,j,k)

                uzik1 = u(1,i,j,k)
                btuz  = bt * uzik1

                t1 = btuz/ac * (r4 + r5)
                t2 = r3 + t1
                t3 = btuz * (r4 - r5)

                rhs(1,i,j,k) = t2
                rhs(2,i,j,k) = -uzik1*r2 + xvel*t2
                rhs(3,i,j,k) =  uzik1*r1 + yvel*t2
                rhs(4,i,j,k) =  zvel*t2  + t3
                rhs(5,i,j,k) =  uzik1*(-xvel*r2 + yvel*r1) +
     >                    qs(i,j,k)*t2 + c2iv*ac2u*t1 + zvel*t3

             end do
          end do
       end do
       if (timeron) call timer_stop(t_tzetar)

       return
       end

max_stars_repo_name: jonarbo/KUBE

max_issues_repo_name: jonarbo/KUBE

max_forks_repo_name: jonarbo/KUBE


text:
      SUBROUTINE EOSALL(S,T,P0,THET,SIGTHE,ALFNEW,BETNEW,GAMNEW,SOUNDV)
C *********************************************************************
C THIS PROGRAMME WRITTEN 8 JULY 1985 BY TREVOR J. McDOUGALL
C EOSALL STANDS FOR "EQUATION OF STATE ALL"
C THIS SUBROUTINE USES THE FOLLOWING FUNCTIONS WRITEN BY BOB MILLARD,
C  - THETA(S,T,P0,PR) ; SVAN(S,T,P0,SIGTHE) ; EOS8D(S,T,P0,DRV)
C THE NEW EXPANSION COEFFICIENT , ALFNEW , DUE TO HEAT , AND THE
C CORRESPONDING SALINE CONTRACTION COEFFICIENT ARE DEFINED IN
C TERMS OF THE TWO CONSERVATIVE PROPERTIES OF SEA WATER,
C NAMELY POTENTIAL TEMPERATURE (REFERED TO ANY REFERENCE LEVEL)
C AND SALINITY. THESE COEFFICIENTS ARE DEFINED IN GILL(1982)
C AND HE LABELS THEM WITH DASHES, SEE HIS SECTION 3.7.4
C ********************************************************************
C
      IMPLICIT DOUBLE PRECISION (A-Z)
      DOUBLE PRECISION DRV(3,8)
C ******* THE REFERENCE PRESSURE PR IS KEPT GENERAL BUT WILL BE
C ******* EQUAL TO ZERO FOR ALL PERCIEVED APPLICATIONS OF NEUTRAL
C ******* SURFACES.
      PR=0.0
      THET=THETA(S,T,P0,PR)
      EDUM=EOS8D(S,T,P0,DRV)
C ******* CALCULATE THE ORDINARY EXPANSION COEFFICIENTS, ALFOLD
C ******* BETOLD AND THE OLD COMPRESSIBILITY, GAMOLD.
      ALFOLD=-DRV(2,3)/(DRV(1,3)+1000.0)
      BETOLD=DRV(1,8)/(DRV(1,3)+1000.0)
      GAMOLD=DRV(1,7)/(DRV(1,3)+1000.0)
C ******* CALCULATE THE SPECIFIC VOLUME ANOMALY, SVAN, AND THE
C ******* SIGMA THETA, SIGTHE, BY THE FUNCTION SVAN.
      SDUM=SVAN(S,THET,PR,SIGTHE)
      CALL DERTHE(S,T,P0,DTHEDT,DTHEDS,DTHEDP)
      ALFNEW=ALFOLD/DTHEDT
      BETNEW=BETOLD+ALFNEW*DTHEDS
      GAMNEW=GAMOLD+ALFNEW*DTHEDP
C ******* CHECK VALUES OF THESE NEW 'EXPANSION COEFFICIENTS'
C ******* AT S=40.0,T=40.0,THET=36.8907,P0=10000.0 DBARS,
C ******* ALFNEW=4395.6E-7 ; (ALFOLD=4181.1E-7)
C ******* BETNEW=6646.9E-7 ; (BETOLD=6653.1E-7)
C ******* GAMNEW=31.4E-7   ; (GAMOLD=32.7E-7)
C *******
C ******* NOW FOR FUN WE CALCULATE THE SPEED OF SOUND.
C ******* THE IN SITU DENSITY IS (DRV(1,3)+1000.0)
      SOUNDV=SQRT(ABS(1.0E+4/(GAMNEW*(DRV(1,3)+1000.0))))
C ******* CHECK VALUE OF THE SPEED OF SOUND IS
C ******* AT S=40.0,T=40.0,THET=36.8907,P0=10000.0 DBARS,
C ******* IS SOUNDV=1734.8 M/S.
      RETURN
      END


text:
Describe Users/DrThomas here.

20130712 14:17:20 nbsp Dr ThomasIts generally frowned upon to delete or alter other peoples comments. You might want to read the Welcome to the wiki/business owner Welcome for Business Owners.  Users/PeteB

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! written by jxzou at 20210111: 3rd order MRPT interfaces

! do CASCI/CASSCF-CASPT3 or CASSCF-NEVPT3 when npair<=7
subroutine do_mrpt3()
 use print_id, only: iout
 use mr_keyword, only: dmrgci, dmrgscf, CIonly, caspt3, nevpt3, casnofch, &
  casscf_prog, casci_prog, bgchg, chgname, mem, nproc, molpro_path, bdf_path,&
  eist
 use mol, only: caspt2_e, nevpt2_e, caspt3_e, nevpt3_e, ptchg_e, nuc_pt_e
 implicit none
 integer :: i, mem0, RENAME, system
 character(len=24) :: data_string
 character(len=240) :: pyname, mklname
 character(len=240) :: string, outname, inpname, inporb
 real(kind=8) :: ref_e, corr2_e, corr3_e
 logical :: alive(2)

 if(eist == 1) return ! excited state calculation
 alive = [caspt3, nevpt3]
 if(ALL(alive .eqv. .false.)) return

 write(iout,'(//,A)') 'Enter subroutine do_mrpt3...'
 mem0 = CEILING(DBLE(mem*125)/DBLE(nproc))

 if((dmrgci .or. dmrgscf)) then
  write(iout,'(A)') 'ERROR in subroutine do_mrpt3: CASPT3/NEVPT3 based on DMRG&
                   & reference is not supported.'
  stop
 end if

 if(.not. CIonly) then
  if(casscf_prog == 'orca') then
   write(iout,'(A)') 'Warning: ORCA is used as the CASSCF solver,&
                    & the NO coefficients in .mkl file are only 7-digits.'
   write(iout,'(A)') 'This will affect the PT3 energy up to 10^-5 a.u.&
                    & Such small error is usually not important.'
   write(iout,'(A)') 'If you care about the accuracy, please use another CASSCF solver.'
  end if

  if(caspt3) then
   string = 'CASPT3 based on CASSCF orbitals.'
  else
   string = 'NEVPT3 based on CASSCF orbitals.'
  end if

 else ! CIonly = .True.
  if(casci_prog == 'orca') then
   write(iout,'(A)') 'Warning: ORCA is used as the CASCI solver,&
                    & the NO coefficients in .mkl file are only 7-digits.'
   write(iout,'(A)') 'This will affect the PT3 energy up to 10^-5 a.u.&
                    & Such small error is usually not important.'
   write(iout,'(A)') 'If you care about the accuracy, please use another CASCI solver.'
  end if

  string = 'CASPT3 based on CASCI orbitals.'
  write(iout,'(A)') 'Warning: the CASSCF orbital optimization is strongly&
                     & recommended'
  write(iout,'(A)') 'to be performed before PT3, unless it is too time-consuming.'
 end if

 write(iout,'(A)') TRIM(string)
 write(iout,'(A)',advance='no') 'Frozen_core = F, '

 if(caspt3) then ! CASPT3
  write(iout,'(A)') 'CASPT3 using program molpro'
  call check_exe_exist(molpro_path)

  i = system('fch2com '//TRIM(casnofch))
  i = index(casnofch, '.fch', back=.true.)
  pyname = casnofch(1:i-1)//'.com'
  string = casnofch
  call convert2molpro_fname(string, '.a')
  i = index(casnofch, '_NO', back=.true.)
  inpname = casnofch(1:i-1)//'_CASPT3.com'
  outname = casnofch(1:i-1)//'_CASPT3.out'
  inporb = inpname
  call convert2molpro_fname(inporb, '.a')
  i = RENAME(TRIM(pyname), TRIM(inpname))
  i = RENAME(TRIM(string), TRIM(inporb))

  call prt_mrpt_molpro_inp(inpname, 3) ! 1/2/3 for NEVPT2/CASPT2/CASPT3
  if(bgchg) i = system('add_bgcharge_to_inp '//TRIM(chgname)//' '//TRIM(inpname))
  write(string,'(2(A,I0),A)') TRIM(molpro_path)//' -n ',nproc,' -t 1 -m ',mem0, &
                           'm '//TRIM(inpname)
  i = system(TRIM(string))

 else ! FIC-NEVPT3
  write(iout,'(A)') 'NEVPT3 using program bdf'
  call check_exe_exist(bdf_path)

  i = system('fch2bdf '//TRIM(casnofch)//' -no')
  i = index(casnofch, '.fch', back=.true.)
  inporb = casnofch(1:i-1)//'_bdf.inporb'
  string = casnofch(1:i-1)//'_bdf.inp'
  i = index(casnofch, '_NO', back=.true.)
  mklname = casnofch(1:i)//'NEVPT3.inporb'
  inpname = casnofch(1:i)//'NEVPT3.inp'
  outname = casnofch(1:i)//'NEVPT3.out'
  i = RENAME(TRIM(inporb), TRIM(mklname))
  i = RENAME(TRIM(string), TRIM(inpname))

  call prt_mrpt_bdf_inp(inpname, 3) ! 1/2/3 for SDSPT2/NEVPT2/NEVPT3
  if(bgchg) i = system('add_bgcharge_to_inp '//TRIM(chgname)//' '//TRIM(inpname))
  i = index(inpname, '.inp', back=.true.)
  string = inpname(1:i-1)
  i = system(TRIM(bdf_path)//' '//TRIM(string))
 end if

 if(i /= 0) then
  if(caspt3) then
   write(iout,'(A)') 'ERROR in subroutine do_mrpt3: CASPT3 computation failed.'
  else
   write(iout,'(A)') 'ERROR in subroutine do_mrpt3: NEVPT3 computation failed.'
  end if
  write(iout,'(A)') 'You can open file '//TRIM(outname)//' and check why.'
  stop
 end if

 if(caspt3) then ! read CASPT3 energy
  call read_caspt3_energy_from_molpro_out(outname, ref_e, corr2_e, corr3_e)
  ref_e = ref_e + ptchg_e
 else            ! read NEVPT3 energy
  call read_nevpt3_energy_from_bdf_out(outname, ref_e, corr2_e, corr3_e)
  ! here the parameter davidson_e is useless
  ref_e = ref_e + ptchg_e + nuc_pt_e
 end if

 write(iout,'(/,A,F18.8,1X,A4)')'E(ref)       = ', ref_e, 'a.u.'

 if(caspt3) then ! CASPT2
  write(iout,'(A)')             'IP-EA shift  =         0.25 (default)'
  caspt2_e = ref_e + corr2_e
  caspt3_e = ref_e + corr3_e
  write(iout,'(A,F18.8,1X,A4)') 'E(corr2)     = ', corr2_e, 'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(CASPT2)    = ', caspt2_e,'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(corr3)     = ', corr3_e, 'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(CASPT3)    = ', caspt3_e,'a.u.'
 else            ! FIC-NEVPT3
  nevpt2_e = ref_e + corr2_e
  nevpt3_e = ref_e + corr3_e
  write(iout,'(A,F18.8,1X,A4)') 'E(corr2)     = ', corr2_e, 'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(FIC-NEVPT2)= ', nevpt2_e,'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(corr3)     = ', corr3_e, 'a.u.'
  write(iout,'(A,F18.8,1X,A4)') 'E(FIC-NEVPT3)= ', nevpt3_e,'a.u.'
 end if

 call fdate(data_string)
 write(iout,'(A)') 'Leave subroutine do_mrpt3 at '//TRIM(data_string)
 return
end subroutine do_mrpt3

! read CASPT3 electronic energy from Molpro .out file
subroutine read_caspt3_energy_from_molpro_out(outname, ref_e, corr2_e, corr3_e)
 implicit none
 integer :: i, fid
 integer, parameter :: iout = 6
 real(kind=8), intent(out) :: ref_e, corr2_e, corr3_e
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 ref_e = 0d0
 corr2_e = 0d0
 corr3_e = 0d0
 open(newunit=fid,file=TRIM(outname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(2:20) == '!MCSCF STATE  1.1 E') exit
 end do ! for while

 if(i /= 0) then
  write(iout,'(A)') 'ERROR in subroutine read_caspt3_energy_from_molpro_out:'
  write(iout,'(A)') "'!MCSCF STATE  1.1 E' not found in file "//TRIM(outname)
  close(fid)
  stop
 end if

 i = index(buf,'ergy')
 read(buf(i+4:),*) ref_e

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(2:20) == '!RSPT2 STATE  1.1 E') exit
 end do ! for while

 if(i /= 0) then
  write(iout,'(A)') 'ERROR in subroutine read_caspt3_energy_from_molpro_out:'
  write(iout,'(A)') "'!RSPT2 STATE  1.1 E' not found in file "//TRIM(outname)
  close(fid)
  stop
 end if

 i = index(buf,'ergy')
 read(buf(i+4:),*) corr2_e
 corr2_e = corr2_e - ref_e

 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(2:20) == '!RSPT3 STATE  1.1 E') exit
 end do ! for while
 close(fid)

 if(i /= 0) then
  write(iout,'(A)') 'ERROR in subroutine read_caspt3_energy_from_molpro_out:'
  write(iout,'(A)') "'!RSPT3 STATE  1.1 E' not found in file "//TRIM(outname)
  stop
 end if

 i = index(buf,'ergy')
 read(buf(i+4:),*) corr3_e
 corr3_e = corr3_e - ref_e
 return
end subroutine read_caspt3_energy_from_molpro_out

! read NEVPT3 electronic energy from BDF .out file
subroutine read_nevpt3_energy_from_bdf_out(outname, ref_e, corr2_e, corr3_e)
 implicit none
 integer :: i, fid
 integer, parameter :: iout = 6
 real(kind=8) :: rtmp
 real(kind=8), intent(out) :: ref_e, corr2_e, corr3_e
 character(len=240) :: buf
 character(len=240), intent(in) :: outname

 ref_e = 0d0
 corr2_e = 0d0
 corr3_e = 0d0

 open(newunit=fid,file=TRIM(outname),status='old',position='rewind')
 do while(.true.)
  read(fid,'(A)',iostat=i) buf
  if(i /= 0) exit
  if(buf(2:6) == 'NROOT') exit
 end do ! for while

 if(i /= 0) then
  write(iout,'(A)') "ERROR in subroutine read_mrpt_nevpt3_from_bdf_out: no&
                    & 'NROOT' found in file "//TRIM(outname)
  close(fid)
  stop
 end if

 read(fid,*) i, ref_e, corr2_e, rtmp, corr3_e
 close(fid)
 corr2_e = corr2_e - ref_e
 corr3_e = corr3_e - ref_e
 return
end subroutine read_nevpt3_energy_from_bdf_out


text:
	SUBROUTINE OAGSTN  ( iflno, source, rlevl, ivert, slat, slon,
     +			     number, iret )
C************************************************************************
C* OAGSTN								*
C*									*
C* This subroutine reads in the latitudes and longitudes of stations	*
C* reporting data in a given area.					*
C*									*
C* OAGSTN  ( IFLNO, SOURCE, RLEVL, IVERT, SLAT, SLON, NUMBER, IRET )	*
C*									*
C* Input parameters:							*
C*	IFLNO		INTEGER		File number			*
C*	SOURCE		CHAR*		Source ( SN or SF )		*
C*	RLEVL		REAL		Vertical level	 		*
C*	IVERT		INTEGER		Vertical coordinates		*
C*	SLAT (NUMBER)	REAL		Latitudes 			*
C*	SLON (NUMBER)	REAL		Longitudes			*
C*	NUMBER		INTEGER		Number of stations		*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 8/85						*
C* M. desJardins/GSFC	11/88	GEMPAK 4.1				*
C* M. desJardins/GSFC	11/89	Changes to station time			*
C* A. Hardy/GSC		 3/99   Added priority parameter to PC_SSTN     *
C* A. Hardy/GSC		 3/99   Added priority parameter to SF_SNXT     *
C* A. Hardy/GSC		 3/99   Removed ispri = 0			*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	source
	REAL		slat (*), slon (*)
	LOGICAL		more
	CHARACTER	cchar*4, chdata*8
	LOGICAL		bounds
	REAL		datarr (LLMXDT)
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret   = 0
	more   = .true.
	number = 0
	IF  ( source .eq. 'SF' )  THEN
	    CALL SF_BEGS  ( iflno, ier )
	  ELSE
	    CALL SN_BEGS  ( iflno, ier )
	END IF
C
C*	Loop through each station.
C
	DO WHILE  ( more )
	    IF  ( source .eq. 'SN' )  THEN
		CALL SN_SNXT  ( iflno, cchar, id, rlat, rlon, re, ier )
		IF  ( ier .ne. 0 )  more = .false.
	      ELSE
		CALL SF_SNXT  ( iflno, cchar, id, rlat, rlon, re, ispri,
     +				ier )
		IF  ( ier .ne. 0 )  more = .false.
	    END IF
C*
	    IF  ( more )  THEN
	        CALL LC_INBN  ( rlat, rlon, bounds, ier )
	      ELSE
	        bounds = .false.
	    END IF
C*
	    IF  ( bounds )  THEN
		IF  ( source .eq. 'SN' )  THEN
		    CALL SN_RDAT  ( iflno, numlev, datarr, ihhmm, ier )
		  ELSE
		    CALL SF_RDAT  ( iflno, datarr, ihhmm, ier )
		    numlev = 1
		END IF
		IF  ( ier .ne. 0 )  bounds = .false.
	    END IF
C*
	    IF  ( bounds )  THEN
		ispri  = 0
		CALL PC_SSTN ( cchar, id, rlat, rlon, re, ispri, ihhmm,
     +				numlev, ier )
		CALL PC_CMVR  ( rlevl, ivert, datarr, data, chdata,
     +				ier )
		IF  ( ( .not. ERMISS ( data ) ) .and.
     +		      ( number .lt. LLSTFL ) )  THEN
		    number = number + 1
		    slat ( number ) = rlat
		    slon ( number ) = rlon
		END IF
	    END IF
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!-----------------------------------------------------------------------
!     Subroutine:    stellopt_fieldlines
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          11/05/2019
!     Description:   This subroutine provides and interface to the
!                    FIELDLINES code.
!-----------------------------------------------------------------------
      SUBROUTINE stellopt_fieldlines(lscreen,iflag)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stellopt_runtime, ONLY: proc_string
      USE fieldlines_grid
      USE fieldlines_lines, ONLY: R_lines, Z_lines, PHI_lines, &
                                  Rhc_lines, Zhc_lines, B_lines
      USE fieldlines_runtime, id_string_fieldlines => id_string, pi_fieldlines => pi, &
                              pi2_fieldlines => pi2, mu0_fieldlines =>mu0, &
                              lverb_fieldlines => lverb
      USE EZspline_obj
      USE EZspline
      USE mpi_params
      USE mpi_inc
      USE mpi_sharmem
      USE wall_mod, ONLY: wall_free

!-----------------------------------------------------------------------
!     Subroutine Parameters
!        iflag         Error flag
!----------------------------------------------------------------------
      IMPLICIT NONE
      LOGICAL, INTENT(in)    :: lscreen
      INTEGER, INTENT(inout) :: iflag
!-----------------------------------------------------------------------
!     Local Variables
!        ier         Error flag
!----------------------------------------------------------------------
      INTEGER ::  ier, nshar

!----------------------------------------------------------------------
!     BEGIN SUBROUTINE
!----------------------------------------------------------------------
      IF (iflag < 0) RETURN

      ! Setup MPI stuff
      myworkid = master
#if defined(MPI_OPT)
      CALL MPI_COMM_DUP( MPI_COMM_MYWORLD, MPI_COMM_FIELDLINES, ierr_mpi)
      CALL MPI_COMM_RANK(MPI_COMM_FIELDLINES, myworkid, ierr_mpi)
      CALL MPI_COMM_SIZE(MPI_COMM_FIELDLINES, nprocs_fieldlines, ierr_mpi)
      CALL MPI_COMM_SPLIT_TYPE(MPI_COMM_FIELDLINES, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, MPI_COMM_SHARMEM, ierr_mpi)
      CALL MPI_COMM_RANK(MPI_COMM_SHARMEM, myid_sharmem, ierr_mpi)
      CALL MPI_COMM_SIZE(MPI_COMM_SHARMEM, nshar, ierr_mpi)
#endif
      IF (lverb_fieldlines) THEN
         WRITE(6,'(a,f5.2)') 'FIELDLINES Version ',FIELDLINES_VERSION
         WRITE(6,'(A)')      '-----  MPI Parameters  -----'
         WRITE(6,'(A,I8)')  '   Nproc_total:  ', nprocs_fieldlines
         WRITE(6,'(A,3X,I5)')  '   Nproc_shared: ', nshar
         CALL FLUSH(6)
      END IF

      ! Initialize fieldlines
      CALL fieldlines_init_vars

      ! Setup run
      IF (myworkid == master) lverb_fieldlines = lscreen
      id_string_fieldlines = TRIM(proc_string)

      ! Initalize grids
      CALL fieldlines_init

      IF (lemc3 .or. lbfield_only .or. lafield_only) nruntype=runtype_norun

      SELECT CASE(nruntype)
         CASE(runtype_old)
            CALL fieldlines_follow
         CASE(runtype_full)
            IF (lverb_fieldlines) THEN
               WRITE(6,'(A)') '===========ROUGH GRID=========='
               WRITE(6,'(A,F8.5,A,F8.5,A)') '   EDGE_INNER[R,Z]   = [',&
                     r_start(1),',',z_start(1),']'
               WRITE(6,'(A,F8.5,A,F8.5,A)') '   EDGE_OUTER[R,Z]   = [',&
                     MAXVAL(r_start,MASK = r_start > 0),',',MAXVAL(z_start,MASK = r_start > 0),']'
            END IF
#if defined(MPI_OPT)
            CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
            IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BARRIER_ERR,'fieldlines_init',ierr_mpi)
#endif
            CALL fieldlines_follow  ! This call on field grid
#if defined(MPI_OPT)
            CALL MPI_BARRIER(MPI_COMM_FIELDLINES,ierr_mpi)
            IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_BARRIER_ERR,'fieldlines_init',ierr_mpi)
#endif
            CALL fieldlines_init_subgrid
            CALL fieldlines_follow  ! This call on subgrid grid
            !CALL fieldlines_periodic_orbits  ! This call on subgrid grid
         CASE(runtype_norun)
      END SELECT

      ! Output Data
      CALL fieldlines_write

      ! Clean up
      IF (lvessel) CALL wall_free(ier, MPI_COMM_FIELDLINES)
      IF (ASSOCIATED(raxis)) CALL mpidealloc(raxis,win_raxis)
      IF (ASSOCIATED(phiaxis)) CALL mpidealloc(phiaxis,win_phiaxis)
      IF (ASSOCIATED(zaxis)) CALL mpidealloc(zaxis,win_zaxis)
      IF (ASSOCIATED(B_R)) CALL mpidealloc(B_R,win_B_R)
      IF (ASSOCIATED(B_PHI)) CALL mpidealloc(B_PHI,win_B_PHI)
      IF (ASSOCIATED(B_Z)) CALL mpidealloc(B_Z,win_B_Z)
      IF (ASSOCIATED(BR4D)) CALL mpidealloc(BR4D,win_BR4D)
      IF (ASSOCIATED(BZ4D)) CALL mpidealloc(BZ4D,win_BZ4D)
      IF (ASSOCIATED(MU4D)) CALL mpidealloc(MU4D,win_MU4D)
      IF (ASSOCIATED(MODB4D)) CALL mpidealloc(MODB4D,win_MODB4D)
      IF (ALLOCATED(R_lines)) DEALLOCATE(R_lines)
      IF (ALLOCATED(Z_lines)) DEALLOCATE(Z_lines)
      IF (ALLOCATED(PHI_lines)) DEALLOCATE(PHI_lines)
      IF (ALLOCATED(Rhc_lines)) DEALLOCATE(Rhc_lines)
      IF (ALLOCATED(Zhc_lines)) DEALLOCATE(Zhc_lines)
      IF (ALLOCATED(B_lines)) DEALLOCATE(B_lines)
      IF (EZspline_allocated(BR_spl)) CALL EZspline_free(BR_spl,ier)
      IF (EZspline_allocated(BZ_spl)) CALL EZspline_free(BZ_spl,ier)
      IF (EZspline_allocated(MU_spl)) CALL EZspline_free(MU_spl,ier)
      IF (EZspline_allocated(MODB_spl)) CALL EZspline_free(MODB_spl,ier)

      CALL MPI_COMM_FREE(MPI_COMM_SHARMEM,ierr_mpi)
      CALL MPI_COMM_FREE(MPI_COMM_FIELDLINES,ierr_mpi)
      IF (lverb_fieldlines) THEN
         write(6,*)'==========  FIELDLINES Complete  =========='
         write(6,*)'==========================================='
         CALL FLUSH(6)
      END IF
      RETURN
!----------------------------------------------------------------------
!     END SUBROUTINE
!----------------------------------------------------------------------
      END SUBROUTINE stellopt_fieldlines

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
module qutraj_general
  !
  ! Global constants and general purpose subroutines
  !

  use qutraj_precision

  implicit none

  !
  ! Constants
  !

  ! error params
  integer, parameter :: blas_error_param    = -23

  ! imaginary unit
  complex(wp), parameter :: ii = (0._wp,1._wp)

  !
  ! Interfaces
  !

  interface new
    module procedure int_array_init
    module procedure int_array_init2
    module procedure wp_array_init
    module procedure wp_array_init2
  end interface

  interface finalize
    module procedure int_array_finalize
    module procedure wp_array_finalize
  end interface

  contains

  !
  ! Initializers and finalizers
  !

  subroutine int_array_init(this,n)
    integer, allocatable, intent(inout) :: this(:)
    integer, intent(in) :: n
    integer :: istat
    if (allocated(this)) then
      deallocate(this,stat=istat)
    endif
    allocate(this(n),stat=istat)
    if (istat.ne.0) then
      call fatal_error("int_array_init: could not allocate.",istat)
    endif
  end subroutine

  subroutine int_array_init2(this,val)
    integer, allocatable, intent(inout) :: this(:)
    integer, intent(in), dimension(:) :: val
    call int_array_init(this,size(val))
    this = val
  end subroutine

  subroutine wp_array_init(this,n)
    real(wp), allocatable, intent(inout) :: this(:)
    integer, intent(in) :: n
    integer :: istat
    if (allocated(this)) then
      deallocate(this,stat=istat)
    endif
    allocate(this(n),stat=istat)
    if (istat.ne.0) then
      call fatal_error("sp_array_init: could not allocate.",istat)
    endif
  end subroutine

  subroutine wp_array_init2(this,val)
    real(wp), allocatable, intent(inout) :: this(:)
    real(wp), intent(in), dimension(:) :: val
    call wp_array_init(this,size(val))
    this = val
  end subroutine

  subroutine int_array_finalize(this)
    integer, allocatable, intent(inout) :: this(:)
    integer :: istat=0
    if (allocated(this)) then
      deallocate(this,stat=istat)
    endif
    if (istat.ne.0) then
      call error("int_array_finalize: could not deallocate.",istat)
    endif
  end subroutine

  subroutine wp_array_finalize(this)
    real(wp), allocatable, intent(inout) :: this(:)
    integer :: istat=0
    if (allocated(this)) then
      deallocate(this,stat=istat)
    endif
    if (istat.ne.0) then
      call error("wp_array_finalize: could not deallocate.",istat)
    endif
  end subroutine

  !
  ! Error handling
  !

  subroutine error(errormsg,ierror)
    character(len=*), intent(in), optional :: errormsg
    integer, intent(in), optional :: ierror
    if (present(errormsg)) then
      write(*,*) 'error: ',errormsg
    endif
    if (present(ierror)) then
      write(*,*) 'error flag=',ierror
    endif
  end subroutine

  subroutine fatal_error(errormsg,ierror)
    character(len=*), intent(in), optional :: errormsg
    integer, intent(in), optional :: ierror
    if (present(errormsg)) then
      write(*,*) 'fatal error: ',errormsg
    endif
    if (present(ierror)) then
      write(*,*) 'error flag=',ierror
    endif
    write(*,*) 'halting'
    stop 1
  end subroutine

end module

max_stars_repo_name: kiuthed/qutip

max_issues_repo_name: trxw/qutip

max_forks_repo_name: trxw/qutip


function erfcc_s(x)

  use mo_kind, only: sp
  use mo_nrutil, only: poly

  implicit none

  real(sp), intent(in) :: x
  real(sp)             :: erfcc_s

  real(sp) :: t,z
  real(sp), dimension(10) :: coef = (/ -1.26551223_sp, 1.00002368_sp, &
       0.37409196_sp, 0.09678418_sp, -0.18628806_sp, 0.27886807_sp, &
       -1.13520398_sp, 1.48851587_sp, -0.82215223_sp, 0.17087277_sp /)

  z = abs(x)
  t = 1.0_sp / (1.0_sp+0.5_sp*z)

  erfcc_s = t * exp(-z*z + poly(t,coef))
  if (x < 0.0_sp) erfcc_s = 2.0_sp - erfcc_s

end function erfcc_s

function erfcc_v(x)

  use mo_kind,   only: sp
  use mo_nrutil, only: poly

  implicit none

  real(sp), dimension(:), intent(in) :: x
  real(sp), dimension(size(x))       :: erfcc_v

  real(sp), dimension(size(x)) :: t,z
  real(sp), dimension(10) :: coef = (/ -1.26551223_sp, 1.00002368_sp, &
       0.37409196_sp, 0.09678418_sp, -0.18628806_sp, 0.27886807_sp, &
       -1.13520398_sp, 1.48851587_sp, -0.82215223_sp, 0.17087277_sp /)

  z = abs(x)
  t = 1.0_sp / (1.0_sp+0.5_sp*z)

  erfcc_v = t * exp(-z*z + poly(t,coef))
  where (x < 0.0_sp) erfcc_v = 2.0_sp - erfcc_v

end function erfcc_v

function derfcc_s(x)

  use mo_kind, only: dp
  use mo_nrutil, only: poly

  implicit none

  real(dp), intent(in) :: x
  real(dp)             :: derfcc_s

  real(dp) :: t,z
  real(dp), dimension(10) :: coef = (/ -1.26551223_dp, 1.00002368_dp, &
       0.37409196_dp, 0.09678418_dp, -0.18628806_dp, 0.27886807_dp, &
       -1.13520398_dp, 1.48851587_dp, -0.82215223_dp, 0.17087277_dp /)

  z = abs(x)
  t = 1.0_dp / (1.0_dp+0.5_dp*z)

  derfcc_s = t * exp(-z*z + poly(t,coef))
  if (x < 0.0_dp) derfcc_s = 2.0_dp - derfcc_s

end function derfcc_s

function derfcc_v(x)

  use mo_kind,   only: dp
  use mo_nrutil, only: poly

  implicit none

  real(dp), dimension(:), intent(in) :: x
  real(dp), dimension(size(x))       :: derfcc_v

  real(dp), dimension(size(x)) :: t,z
  real(dp), dimension(10) :: coef = (/ -1.26551223_dp, 1.00002368_dp, &
       0.37409196_dp, 0.09678418_dp, -0.18628806_dp, 0.27886807_dp, &
       -1.13520398_dp, 1.48851587_dp, -0.82215223_dp, 0.17087277_dp /)

  z = abs(x)
  t = 1.0_dp / (1.0_dp+0.5_dp*z)

  derfcc_v = t * exp(-z*z + poly(t,coef))
  where (x < 0.0_dp) derfcc_v = 2.0_dp - derfcc_v

end function derfcc_v


text:
Describe Users/aconn here.

20091102 12:34:37 nbsp Hi there. Im sorry to hear about your Lost Pets lost cat. I created a new page on the wiki:woodland Woodland Wiki, and moved the information about Kiki over there. If Kiki has been found, feel free to delete the information, or move to the associated wiki:woodland Lost_Pets/Success_Stories page. Users/robinlaughlin


text:
C
C
      SUBROUTINE ANUSET(
C         inputs
     &    NEOS, URINP,
C         outputs
     &    UARR)
C
C***********************************************************************
C
C     User defined EOS setup entry point. See ANUEOS for running
C     entry point.
C
C     WARNING - depending on the computer system
C     All variables are single or double precision.
C
C     Inputs: (all units cgs-ev)
C
C     NEOS   =  User input eos number (negative)
C
C     URINP  =  input array
C
C     Output:
C
C     UARR   =  aneos data array for this material
C
C     The input array is the aneos input record - either 24 or 40
C     variables are defined depending on aneos input record format
C     that you selected. Of these the following must be defined:
C
C     urinp(1) = number of elements in this material
C     urinp(2) = -4.  (to get here)
C     urinp(3) = reference density for this material
C     urinp(4) = reference temperature for this material
C     The others can be anything that you want to input to your model.
C
C     The other input array is of length 99. Certain location in this
C     array must not be redefined by your model. Any of the remaining
C     locations can be used by your model. Note that this array is
C     written to the restart file so that you don't have to do anything
C     about restarts. * indicates values that you MUST set.
C
C     UARR location       value
C         11        reference density
C         12        reference temperature
C     *   18        melt temperature
C         22        zero
C         23        internal value
C         26        average atomic number
C         27        number of atoms per gram
C         28        number of elements
C         29        average atomic weight
C         30        internal switch
C         31        internal pointer
C     *   47        liquid density at melt point for solid form
C         64        internal value
C         65        internal value
C         71        internal value
C         72        internal value
C         73        internal value
C         74        internal value
C
C     EOS values must be defined for all positive temperatures and
C     densities. Returned values must be realistic - negative dP/dRHO
C     is not allowed.
C
C************************************** 4/90 version slt ***************
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION URINP(48), UARR(99)
C     klst is the output listing file
      COMMON /FILEOS/ KLST,KINP
      SAVE   /FILEOS/
      PARAMETER (ZERO=0.D0)
C
C     add your code here
C
C     example code for ideal gas ---------------------------------------
C     input assumed
C     aneos #elements -4  rho0   temp0   gamma   Cv
C     aneos 0
C     aneos 0
C     aneos element records
C
      UARR(1)=URINP(5)
      UARR(2)=URINP(6)
      UARR(18)=ZERO
C
C     example code for ideal gas ---------------------------------------
C
      RETURN
      END

max_stars_repo_name: chreinhardt/M-ANEOS

max_issues_repo_name: chreinhardt/M-ANEOS

max_forks_repo_name: chreinhardt/M-ANEOS


text:
      SUBROUTINE BMG2_SymStd_SETUP_MSG(&
     &                pMSG, pMSGSO, iMSG_Geom, NMSGi, pSI_MSG,&
     &                IBC, IGRD, NOG, NOGm, NProc, MyProc,  &
     &                DimX, DimY, DimXfine, DimYfine,&
     &                ProcGrid, NProcI, NProcJ, MPICOMM&
     &                ) BIND(C, NAME='BMG2_SymStd_SETUP_MSG')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     Create pointers into the integer work space describing the storage
!     of MSG communication setup.  Note that the pointer shift, pSI_MSG,
!     is changed in every call.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================
      USE ModInterface
      IMPLICIT NONE

! -----------------------------
!     Includes
!
      INCLUDE 'mpif.h'
      INCLUDE 'MSG_f90.h'

      INCLUDE 'BMG_parameters_f90.h'
      INCLUDE 'BMG_workspace_f90.h'

! ---------------------------
!    Argument Declarations:
!
      INTEGER(C_INT),  VALUE :: IBC, NOG, NOGm, NProcI, NProcj
      INTEGER(len_t), VALUE :: NMSGi
      INTEGER, VALUE :: MPICOMM
      INTEGER(len_t) :: IGRD(NOGm,NBMG_pIGRD), iMSG_Geom(NMSGi)
      INTEGER(C_INT) :: pMSG(NBMG_pMSG,NOG), pMSGSO(NBMG_pMSG,NOG)
      INTEGER(len_t) :: DimX(NProcI,NOGm), DimY(NProcJ,NOGm)
      INTEGER(len_t) :: DimXfine(NProcI), DimYfine(NProcJ)
      INTEGER(C_INT) ::  ProcGrid(NProcI,NProcJ)
      INTEGER(C_INT) :: pSI_MSG
      INTEGER(C_INT), VALUE :: NProc, MyProc

! --------------------------
!     Local Declarations:
!
      integer(len_t) :: I, J, iGs, jGs
      Integer ierror, N

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

      call MSG_set_comm_parent(MPICOMM)
      call MSG_enable(MyProc, NProc)
!     Note that the pointer shift, pSI_MSG, is changed in every call

      CALL BMG2_SymStd_SETUP_PtrMSG(&
     &     IGRD(NOG,idL_BMG_NLx), IGRD(NOG,idL_BMG_NLy), pSI_MSG,&
     &     NProc, NOG, pMSG)

      CALL BMG2_SymStd_SETUP_PtrMSGSO(&
     &     IGRD(NOG,idL_BMG_NLx), IGRD(NOG,idL_BMG_NLy), pSI_MSG,&
     &     NProc, NOG, pMSGSO)

! -------------------------------------------
!     Compute the dimensions for all local
!     grids on all processors
! -------------------------------------------

      DO I=1, NProcI
         DimX(I,NOG) = DimXfine(I)
      END DO

      DO J=1, NProcJ
         DimY(J,NOG) = DimYfine(J)
      END DO

      DO N=NOG-1, 1, -1

         iGs = 1

         DO I=1, NProcI

            if (mod(iGs,2).eq.1) then
               DimX(I,N) = (DimX(I,N+1)+1)/2
            else
               IF (mod(DimX(I,N+1),2).eq.1) THEN
                  DimX(I,N) = (DimX(I,N+1)-1)/2
               ELSE
                  DimX(I,N) = (DimX(I,N+1)+1)/2
               ENDIF
            endif

            iGs = iGs + DimX(I,N+1)

         END DO

         jGs = 1

         DO J=1, NProcJ

            if (mod(jGs,2).eq.1) then
               DimY(J,N) = (DimY(J,N+1)+1)/2
            else
               IF (mod(DimY(J,N+1),2).eq.1) THEN
                  DimY(J,N) = (DimY(J,N+1)-1)/2
               ELSE
                  DimY(J,N) = (DimY(J,N+1)+1)/2
               ENDIF
            endif

            jGs = jGs + DimY(J,N+1)

         END DO

      END DO

! -------------------------------------------
!     Create MSG grid information for U
! -------------------------------------------

      DO N=NOG, 1, -1
         !
         CALL BMG2_SymStd_SETUP_MSGGrid(&
     &        IGRD(N,idL_BMG_NGx), IGRD(N,idL_BMG_NGy), IBC,&
     &        iMSG_Geom(pMSG(ipL_MSG_LocalArraySize,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_GlobalCoordLocalData,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_GlobalCoordActData,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_ActDataStart,N)),&
     &        DimX, DimY, ProcGrid, MyProc, &
     &        NProc, NProcI, NProcJ, NOGm, N,&
     &        MPICOMM)
         !
      ENDDO

! --------------------------------------------
!     Create MSG grid information for SO
! --------------------------------------------

      DO N=NOG, 1, -1
         !
         CALL BMG2_SymStd_SETUP_MSGGridSO(&
     &        IGRD(N,idL_BMG_NGx), IGRD(N,idL_BMG_NGy), IBC,&
     &        iMSG_Geom(pMSGSO(ipL_MSG_LocalArraySize,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_GlobalCoordLocalData,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_GlobalCoordActData,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_ActDataStart,N)),&
     &        DimX, DimY, ProcGrid, NProc, NProcI, NProcJ, NOGm, N,&
     &        MPICOMM )
         !
      ENDDO

!
!     Setup MSG communications for all grids
!

      DO N=NOG, 1, -1

!         write(*,*) N, NProc, MyProc

         iMSG_Geom(pMSG(ipL_MSG_NumAdjProc, N)) = 0

!         write(*,*)MyProc,':', (iWork(pMSG(ipL_MSG_ActDataStart,N)+i),
!     &        i=0,3*NProc-1)

         CALL MSG_tp_setup (&
     &        iMSG_Geom(pMSG(ipL_MSG_LocalArraySize,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_ActDataStart,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_GlobalCoordLocalData,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_GlobalCoordActData,N)),&
     &        NProc, MyProc,&
     &        iMSG_Geom(pMSG(ipL_MSG_NumAdjProc,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_Proc,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_Ipr,N)),&
     &        iMSG_Geom(pMSG(ipL_MSG_Index,N)),&
     &        1,1,ierror)

!         write(*,*) 'after MSG_tp_setup'

!         write(*,*) N,IGRD(N,idL_BMG_ICOORD),IGRD(N,idL_BMG_JCOORD)

         iMSG_Geom(pMSGSO(ipL_MSG_NumAdjProc, N)) = 0

         CALL MSG_tp_setup (&
     &        iMSG_Geom(pMSGSO(ipL_MSG_LocalArraySize,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_ActDataStart,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_GlobalCoordLocalData,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_GlobalCoordActData,N)),&
     &        NProc, MyProc,&
     &        iMSG_Geom(pMSGSO(ipL_MSG_NumAdjProc,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_Proc,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_Ipr,N)),&
     &        iMSG_Geom(pMSGSO(ipL_MSG_Index,N)),&
     &        1,1,ierror)

      ENDDO

      return
      end

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


text:
      SUBROUTINE VRFVAL ()

C     Bug found by Corinne on 4 October 2002 (not reproduced)

      REAL*8 X,Y,I

      IF (MAX(ABS(X),ABS(Y)).GT.1.) I = 1

      Z = ABS(X)
      U = ABS(Y)
      V = AMAX1(Z, U)

      if(V.gt.1.) then
C     Let's see what information is generated:
         print *, x, y, z, u, v
C     And what's left about X and Y? Nothing!
         read *, z, u, v
         print *, x, y
      else
C     Let's see again what information is generated:
         print *, x, y, z, u, v
C     And what's left about X and Y? Something!
         read *, z, u, v
         print *, x, y
      endif

      print *, x, y, z, u, v
      END


text:
      PROGRAM xsplie2
C     driver for routine splie2
      INTEGER M,N
      PARAMETER(M=10,N=10)
      INTEGER i,j
      REAL x1x2,x1(M),x2(N),y(M,N),y2(M,N)
      do 11 i=1,M
        x1(i)=0.2*i
11    continue
      do 12 i=1,N
        x2(i)=0.2*i
12    continue
      do 14 i=1,M
        do 13 j=1,N
          x1x2=x1(i)*x2(j)
          y(i,j)=x1x2**2
13      continue
14    continue
      call splie2(x1,x2,y,M,N,y2)
      write(*,'(/1x,a)') 'Second derivatives from SPLIE2'
      write(*,'(1x,a/)') 'Natural spline assumed'
      do 15 i=1,5
        write(*,'(1x,5f12.6)') (y2(i,j),j=1,5)
15    continue
      write(*,'(/1x,a/)') 'Actual second derivatives'
      do 17 i=1,5
        do 16 j=1,5
          y2(i,j)=2.0*(x1(i)**2)
16      continue
        write(*,'(1x,5f12.6)') (y2(i,j),j=1,5)
17    continue
      END

ext: for

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


text:
c     check the value of function arguments
      subroutine foo
      integer a, b, c
      a = 1
      b = 2
      c = bar (a, b)     ! 1, 2
      d = bar (a, b)     ! Bot, Bot
      e = bar (1 + 2, c) ! 3, Bot
      return
      end

      integer function bar(x, y)
      integer x, y
      integer bar
      bar = x + y
      return
      end

max_stars_repo_name: camfort/fortran-vars

max_issues_repo_name: camfort/fortran-vars

max_forks_repo_name: camfort/fortran-vars


text:
module module_shapiro

contains

subroutine shapiro(j,k,wet,etan,eps,eta)
      integer(4), parameter :: nx=500
      integer(4), parameter :: ny=500
!local parameters
      integer, intent(InOut) :: j
      integer, intent(InOut) :: k
      integer, dimension(0:ny+1,0:nx+1), intent(In) :: wet
      real, dimension(0:ny+1,0:nx+1), intent(In) :: etan
      real, intent(In) :: eps
      real, dimension(0:ny+1,0:nx+1), intent(Out) :: eta
real :: term1
real :: term2
real :: term3
! 1-order Shapiro filter
do j = 1,ny
do k = 1,nx
if(wet(j,k)==1)then
  term1 = (1.0-0.25*eps*(wet(j,k+1)+wet(j,k-1)+ wet(j+1,k)+wet(j-1,k)))*etan(j,k)
  term2 = 0.25*eps*(wet(j,k+1)*etan(j,k+1)+wet(j,k-1)*etan(j,k-1))
  term3 = 0.25*eps*(wet(j+1,k)*etan(j+1,k)+wet(j-1,k)*etan(j-1,k))
  eta(j,k) = term1+term2+term3
else
  eta(j,k) = etan(j,k)
end if
end do
end do
end subroutine shapiro

end module module_shapiro

max_stars_repo_name: wimvanderbauwhede/AutoParallel-Fortran

max_issues_repo_name: wimvanderbauwhede/AutoParallel-Fortran

max_forks_repo_name: wimvanderbauwhede/AutoParallel-Fortran