! Chemfiles, an efficient IO library for chemistry file formats
! Copyright (C) 2015-2019 Guillaume Fraux -- BSD licence
!
! =========================================================================== !
! !!!! AUTO-GENERATED FILE !!!! Do not edit. See bindgen repository for the
! generating code (https://github.com/chemfiles/bindgen).
! This file contains Fortran 2003 ISO C Binding interface to the C API
! =========================================================================== !

interface
! Function "chfl_selection", at selection.h:24
function c_chfl_selection(selection) bind(C, name="chfl_selection")
    use iso_c_binding

    implicit none
    type(c_ptr) :: c_chfl_selection
    character(len=1, kind=c_char), dimension(*), intent(in) :: selection
end function

! Function "chfl_selection_copy", at selection.h:37
function c_chfl_selection_copy(selection) bind(C, name="chfl_selection_copy")
    use iso_c_binding

    implicit none
    type(c_ptr) :: c_chfl_selection_copy
    type(c_ptr), value, intent(in) :: selection
end function

! Function "chfl_selection_size", at selection.h:49
function c_chfl_selection_size(selection, size) bind(C, name="chfl_selection_size")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_selection_size
    type(c_ptr), value, intent(in) :: selection
    integer(kind=c_int64_t), intent(inout) :: size
end function

! Function "chfl_selection_string", at selection.h:62
function c_chfl_selection_string(selection, string, buffsize) bind(C, name="chfl_selection_string")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_selection_string
    type(c_ptr), value, intent(in) :: selection
    character(len=1, kind=c_char), dimension(*) :: string
    integer(kind=c_int64_t), value :: buffsize
end function

! Function "chfl_selection_evaluate", at selection.h:75
function c_chfl_selection_evaluate(selection, frame, n_matches) bind(C, name="chfl_selection_evaluate")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_selection_evaluate
    type(c_ptr), value :: selection
    type(c_ptr), value, intent(in) :: frame
    integer(kind=c_int64_t), intent(inout) :: n_matches
end function

! Function "chfl_selection_matches", at selection.h:87
function c_chfl_selection_matches(selection, matches, n_matches) bind(C, name="chfl_selection_matches")
    use iso_c_binding
    import chfl_status

    implicit none
    integer(kind=chfl_status) :: c_chfl_selection_matches
    type(c_ptr), value, intent(in) :: selection
    type(c_ptr), value :: matches
    integer(kind=c_int64_t), value :: n_matches
end function

end interface


text:
!-----subroutine---------------------------------------------g.guidati--
!
     subroutine SETMAT
!
!-----------------------------------------------------------------------
!
!     franco guidati    IAG
!                       Uni Stuttgart
!
!
!     scope             set up aerodynamic influence
!                       coefficients matrix
!
!.......................................................................
!               declarations
!.......................................................................
      USE TINoiseGeneric
      USE TINoiseDDD
      USE TINoiseGeo
      IMPLICIT                        NONE

   ! Local variables.
INTEGER (4)                  :: i,j,ishift,jshift,ng,info

REAL(DbKi)                   :: rhelp,solidangle,summ,herm1,herm2,herm3,herm4
REAL(DbKi)                   :: HERM(n+1,n+1,4),dipok,dv1,dv2,dv1d1,dv1d2
REAL(DbKi)                   :: dv2d1,dv2d2

!.......................................................................
!
!.......................................................................
      ng = 4
      rhelp = acos(-tc1(1)*tc1(n+1) - tc2(1)*tc2(n+1))
      solidangle = (1.0 - 0.5*(2.0*rhelp/pi2))

      do i=1,n+2
        do j=1,n+2
	  Kern(i,j) = 0.0
	end do
      end do

      do i=1,N
!        write(*,*) i
        x1 = yc1(i)
        x2 = yc2(i)
!        write(55,*) real(x1),real(x2),real(nc1(i)),real(nc2(i)),
!     >              real(tc1(i)),real(tc2(i))

        Kern(i,i) = 0.5d0
        if(i.eq.1) Kern(i,i) = solidangle

        summ = 0.0d0
        do j=1,N
          s1 = swork(j)   + 0.0000001d0
          s2 = swork(j+1) - 0.0000001d0
          call CDI0   (swork, yc1, yc2, d2yc1, d2yc2,n, ng, td, Ad, pi2i,herm1, herm2, herm3, herm4 )

          HERM(i,j,1) = herm1
          HERM(i,j,2) = herm2
          HERM(i,j,3) = herm3
          HERM(i,j,4) = herm4

          summ = summ + HERM(i,j,1) + HERM(i,j,2)

!          write(*,*) HERM(i,j,3),HERM(i,j,4)
          do jshift=0,1

            Kern(i,j+jshift) = Kern(i,j+jshift) - HERM(i,j,jshift+1)

            ishift=0
            if(j+jshift.eq.1     ) ishift= 2
            if(j+jshift.eq.2     ) ishift= 1
            if(j+jshift.eq.n     ) ishift=-1
            if(j+jshift.eq.n+1   ) ishift=-2
  	    Kern(i,j-2+ishift+jshift) = Kern(i,j-2+ishift+jshift) -dst(-2,ishift)*HERM(i,j,jshift+3)*.5

            Kern(i,j-1+ishift+jshift) = Kern(i,j-1+ishift+jshift) -dst(-1,ishift)*HERM(i,j,jshift+3)*.5

            Kern(i,j  +ishift+jshift) = Kern(i,j  +ishift+jshift) -dst( 0,ishift)*HERM(i,j,jshift+3)*.5

            Kern(i,j+1+ishift+jshift) = Kern(i,j+1+ishift+jshift) -dst( 1,ishift)*HERM(i,j,jshift+3)*.5

	    Kern(i,j+2+ishift+jshift) = Kern(i,j+2+ishift+jshift) -dst( 2,ishift)*HERM(i,j,jshift+3)*.5

          enddo

        enddo
!	if(i.eq.1) then
!  	  write(*,*) i,summ+1.0
!	  write(*,*) rhelp,solidangle
!	endif

        call calll(x1,x2,yc1(1),yc2(1),ywinf1,ywinf2,ywn1,ywn2,dipok,dv1,dv2,dv1d1,dv1d2,dv2d1,dv2d2)

        Kern(i,n+2) = -dipok

      enddo

      if(.false.) then
      Kern(n+1,1)      = -25.0d0/12.0
      Kern(n+1,2)      =  48.0d0/12.0
      Kern(n+1,3)      = -36.0d0/12.0
      Kern(n+1,4)      =  16.0d0/12.0
      Kern(n+1,5)      =  -3.0d0/12.0
      Kern(n+1,n+1)    =  25.0d0/12.0
      Kern(n+1,n  )    = -48.0d0/12.0
      Kern(n+1,n-1)    =  36.0d0/12.0
      Kern(n+1,n-2)    = -16.0d0/12.0
      Kern(n+1,n-3)    =   3.0d0/12.0
      else
      Kern(n+1,1)      =  2.0
      Kern(n+1,2)      =  -3.0
      Kern(n+1,3)      =  1.0
      Kern(n+1,n+1)    =  -2.0
      Kern(n+1,n  )    =  3.0
      Kern(n+1,n-1)    =  -1.0
      endif

      Kern(n+2,1)   =  1.0
      Kern(n+2,n+1) = -1.0
      Kern(n+2,n+2) = -1.0
      Kern(1,n+2) = -solidangle/2.0

      call dgetrf (n+2,n+2,Kern,n+2,ipiv,info)

!.......................................................................
!               end of subroutine
!.......................................................................
     return
!***********************************************************************
      end

max_stars_repo_name: CarloSucameli/Openfast-noise

max_issues_repo_name: CarloSucameli/Openfast-noise

max_forks_repo_name: CarloSucameli/Openfast-noise


text:
      subroutine g1dat(ier,noutpt,nttyo,udastr,var)
c
c     This subroutine reads a number (var) from a string (udastr)
c
c     This subroutine is called by:
c
c       EQPT/gnenb.f
c       EQPT/rdpca.f
c       EQPT/rdpth.f
c       EQPT/rdpni.f
c       EQPT/rdpn2.f
c       EQPT/rdpnn.f
c
c-----------------------------------------------------------------------
c
c     Principal input:
c
c       udastr = a string containing one numerical data field
c
c     Principal output:
c
c       ier    = error flag
c       var    = the number contained in that field
c
c-----------------------------------------------------------------------
c
      implicit none
c
c-----------------------------------------------------------------------
c
c     Calling sequence variable declarations.
c
      integer noutpt,nttyo
c
      integer ier
c
      character(len=80) udastr
c
      real(8) var
c
c-----------------------------------------------------------------------
c
c     Local variable declarations.
c
      integer jj,jlen
c
      integer ilnobl
c
c-----------------------------------------------------------------------
c
      var = 0.
      ier = 0
c
      call lejust(udastr)
      jj = index(udastr,' ')
      if (jj .eq. 0) jj = 81
      if (jj .lt. 80) udastr(jj:80) = ' '
      jlen = ilnobl(udastr)
c
      if (jlen .gt. 25) then
        write (noutpt,1000) udastr(1:jlen)
        write (nttyo,1000) udastr(1:jlen)
 1000   format(/' * Error - (EQPT/g1dat) Have found a data',
     $  ' field in the string:',/7x,'"',a,'"',/7x,'that appears',
     $  ' to exceed the allowed 25 characters.')
        ier = 1
      else
        read (udastr,1010,err=995) var
 1010   format(e25.18)
      endif
      go to 999
c
  995 ier = 1
c
  999 continue
      end

ext: f


text:
! BSD 2-Clause License
!
! Copyright (c) 2019, Oscar Riveros - www.peqnp.science.
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

module utils
    implicit none
    public :: sum_matrix
contains

    function sum_matrix(matrix, image, num_images) result(sum)
        integer :: i, j, image, num_images
        real :: sum
        real, allocatable :: matrix(:, :)[:]

        sum = 0.0
        do i = image, size(matrix, 1), num_images
            do j = 1, size(matrix, 2)
                sum = sum + matrix(i, j)[image]
            end do
        end do

    end function sum_matrix

end module utils

max_stars_repo_name: maxtuno/Distributed-Algorithms

max_issues_repo_name: maxtuno/Distributed-Algorithms

max_forks_repo_name: maxtuno/Distributed-Algorithms


text:
      double precision function  fdisc(blockno,xc,yc)
      implicit none

      double precision xc,yc, xp, yp, zp
      integer blockno
      integer*8 cont, get_context

      double precision pi, pi2
      common /compi/ pi, pi2

      double precision beta, theta(2)
      common /annulus_comm/ beta, theta

      double precision init_radius
      common /initradius_comm/ init_radius

      double precision x0, y0, r0, r, ravg, th

      cont = get_context()

      call fclaw2d_map_c2m(cont,
     &      blockno,xc,yc,xp,yp,zp)

      ravg = (1 + beta)/2.d0
      th = pi2*(0.25 + 1.d0/32.d0)
      x0 = ravg*cos(th)
      y0 = ravg*sin(th)

      r = sqrt((xp - x0)**2 + (yp-y0)**2)

      r0 = init_radius
      fdisc = r - r0

      end

max_stars_repo_name: scivision/forestclaw

max_issues_repo_name: scivision/forestclaw

max_forks_repo_name: scivision/forestclaw


text:
      SUBROUTINE THERMO(A,B,C,LINEAR,SYM,WT,VIBS,NVIBS,ESCF)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION VIBS(*)
      LOGICAL LINEAR
      CHARACTER*80 KEYWRD, KOMENT, TITLE, TMPKEY*80
      COMMON /KEYWRD/ KEYWRD
      COMMON /TITLES/ KOMENT,TITLE
C
C
C   THERMO CALCULATES THE VARIOUS THERMODYNAMIC QUANTITIES FOR A
C   SPECIFIED TEMPERATURE GIVEN THE VIBRATIONAL FREQUENCIES, MOMENTS OF
C   INERTIA, MOLECULAR WEIGHT AND SYMMETRY NUMBER.
C
C   REFERENCE: G.HERZBERG MOLECULAR SPECTRA AND MOLECULAR STRUCTURE
C              VOL 2, CHAP. 5
C
C   ----    TABLE OF SYMMETRY NUMBERS    ----
C
C        C1 CI CS     1      D2 D2D D2H  4       C(INF)V   1
C        C2 C2V C2H   2      D3 D3D D3H  6       D(INF)H   2
C        C3 C3V C3H   3      D4 D4D D4H  8       T TD     12
C        C4 C4V C4H   4      D6 D6D D6H  12      OH       24
C        C6 C6V C6H   6      S6          3
C
C
C   PROGRAM LIMITATIONS:  THE EQUATIONS USED ARE APPROPRIATE TO THE
C   HIGH TEMPERATURE LIMIT AND WILL BEGIN TO BE INADEQUATE AT TEMPERA-
C   TURES BELOW ABOUT 100 K.  SECONDLY THIS PROGRAM IS ONLY APPROPRIATE
C   IN THE CASE OF MOLECULES IN WHICH THERE IS NO FREE ROTATION
C
C
C
C
*******************************************************************
*
*  THE FOLLOWING CONSTANTS ARE NOW DEFINED:
*          PI  = CIRCUMFERENCE TO DIAMETER OF A CIRCLE
*          R   = GAS CONSTANT IN CALORIES/MOLE
*          H   = PLANCK'S CONSTANT IN ERG-SECONDS
*          AK  = BOLTZMANN CONSTANT IN ERG/DEGREE
*          AC  = SPEED OF LIGHT IN CM/SEC
*******************************************************************
      DIMENSION TRANGE(300)
      DATA PI /3.14159D0 /
      DATA R/1.98726D0/
      DATA H/6.626D-27/
      DATA AK/1.3807D-16/
      DATA AC/2.99776D+10/
*******************************************************************
      IT1=200
      IT2=400
      ISTEP=10
      TMPKEY=KEYWRD
      I=INDEX(TMPKEY,'THERMO(')
      IF(I.NE.0) THEN
C
C   ERASE ALL TEXT FROM TMPKEY EXCEPT THERMO DATA
C
         TMPKEY(:I)=' '
         TMPKEY(INDEX(TMPKEY,')'):)=' '
         IT1=READA(TMPKEY,I)
         IF(IT1.LT.100) THEN
            WRITE(6,'(//10X,''TEMPERATURE RANGE STARTS TOO LOW,'',
     1'' LOWER BOUND IS RESET TO 30K'')')
            IT1=100
         ENDIF
         I=INDEX(TMPKEY,',')
         IF(I.NE.0) THEN
            TMPKEY(I:I)=' '
            IT2=READA(TMPKEY,I)
            IF(IT2.LT.IT1) THEN
               IT2=IT1+200
               ISTEP=10
               GOTO 10
            ENDIF
            I=INDEX(TMPKEY,',')
            IF(I.NE.0) THEN
               TMPKEY(I:I)=' '
               ISTEP=READA(TMPKEY,I)
               IF(ISTEP.LT.1)ISTEP=1
            ELSE
               ISTEP=(IT2-IT1)/20
               IF(ISTEP.EQ.0)ISTEP=1
               IF(ISTEP.GE.2.AND. ISTEP.LT.5)ISTEP=2
               IF(ISTEP.GE.5.AND. ISTEP.LT.10)ISTEP=5
               IF(ISTEP.GE.10.AND. ISTEP.LT.20)ISTEP=10
               IF(ISTEP.GT.20.AND. ISTEP.LT.50)ISTEP=20
               IF(ISTEP.GT.50.AND. ISTEP.LT.100)ISTEP=50
               IF(ISTEP.GT.100)ISTEP=100
            ENDIF
         ELSE
            IT2=IT1+200
         ENDIF
      ENDIF
   10 CONTINUE
      WRITE(6,'(//,A)')TITLE
      WRITE(6,'(A)')KOMENT
      IF(LINEAR) THEN
         WRITE(6,'(//10X,''MOLECULE IS LINEAR'')')
      ELSE
         WRITE(6,'(//10X,''MOLECULE IS NOT LINEAR'')')
      ENDIF
      WRITE(6,'(/10X,''THERE ARE'',I3,'' GENUINE VIBRATIONS IN THIS '',
     1''SYSTEM'')')NVIBS
      WRITE(6,20)
   20 FORMAT(10X,'THIS THERMODYNAMICS CALCULATION IS LIMITED TO',/
     110X,'MOLECULES WHICH HAVE NO INTERNAL ROTATIONS'//)
      WRITE(6,'(//20X,''CALCULATED THERMODYNAMIC PROPERTIES'')')
      WRITE(6,'(42X,''*'')')
      WRITE(6,'(''   TEMP. (K)   PARTITION FUNCTION   H.O.F.'',
     1''    ENTHALPY   HEAT CAPACITY  ENTROPY'')')
      WRITE(6,'(  ''                                    KCAL/MOL'',
     1''   CAL/MOLE    CAL/K/MOL   CAL/K/MOL'',/)')
      DO 30 I=1,NVIBS
   30 VIBS(I)=ABS(VIBS(I))
      ILIM=1
      DO 40 ITEMP=IT1,IT2,ISTEP
         ILIM=ILIM+1
   40 TRANGE(ILIM)=ITEMP
      TRANGE(1)=298.D0
      DO 80 IR=1,ILIM
         ITEMP=TRANGE(IR)
         T=ITEMP
C   ***   INITIALISE SOME VARIABLES   ***
         C1=H*AC/AK/T
         QV=1.0D0
         HV=0.0D0
         E0=0.0D0
         CPV=0.0D0
         SV1=0.0D0
         SV2=0.0D0
C   ***   CONSTRUCT THE FREQUENCY DEPENDENT PARTS OF PARTITION FUNCTION
         DO 50 I=1,NVIBS
            WI=VIBS(I)
            EWJ=EXP(-WI*C1)
            QV=QV/(1-EWJ)
            HV=HV+WI*EWJ/(1-EWJ)
            E0=E0+WI
            CPV=CPV+WI*WI*EWJ/(1-EWJ)/(1-EWJ)
            SV1=SV1+LOG(1.0D0-EWJ)
   50    SV2=SV2+WI*EWJ/(1-EWJ)
C   ***   FINISH CALCULATION OF VIBRATIONAL PARTS   ***
         HV=HV*R*H*AC/AK
         E0=E0*1.4295D0
         CPV=CPV*R*C1*C1
         SV=SV2*R*C1-R*SV1
C   ***   NOW CALCULATE THE ROTATIONAL PARTS  (FIRST LINEAR MOLECULES
         IF(.NOT.LINEAR) GOTO 60
         QR=1/(C1*A*SYM)
         HR=R*T
         CPR=R
         SR=R*(LOG(T*AK/(H*AC*A*SYM)))+R
         GOTO 70
   60    QR=SQRT(PI/(A*B*C*C1*C1*C1))/SYM
         HR=3.0D0*R*T/2.0D0
         CPR=3.0D0*R/2.0D0
         SR=0.5D0*R*(3.D0*LOG(T*AK/(H*AC))
     1-2.D0*LOG(SYM)+LOG(PI/(A*B*C))+3.D0)
   70    CONTINUE
C   ***   CALCULATE INTERNAL CONTRIBUTIONS   ***
         QINT=QV*QR
         HINT=HV+HR
         CPINT=CPV+CPR
         SINT=SV+SR
C   ***   CONSTRUCT TRANSLATION CONTRIBUTIONS   ***
         QTR=(SQRT(2.D0*PI*WT*T*AK*1.6606D-24)/H)**3
         HTR=5.0D0*R*T/2.0D0
         CPTR=5.0D0*R/2.0D0
         STR=2.2868D0*(5.0D0*LOG10(T)+3.0D0*LOG10(WT))-2.3135D0
C   ***   CONSTRUCT TOTALS   ***
         CPTOT=CPTR+CPINT
         STOT=STR+SINT
         HTOT=HTR+HINT
C   ***   OUTPUT SECTION   ***
         IF(IR.EQ.1)THEN
            H298=HTOT
         ELSE
            WRITE(6,'(/,I7,''  VIB.'',G18.4
     1           ,13X,3F11.5        )')ITEMP,QV,  HV,  CPV,  SV
            WRITE(6,'(7X,''  ROT.'',G13.3
     1           ,16X,3F11.3        )')      QR,  HR,  CPR,  SR
            WRITE(6,'(7X,''  INT.'',G13.3
     1           ,16X,3F11.3        )')      QINT,HINT,CPINT,SINT
            WRITE(6,'(7X,''  TRA.'',G13.3
     1           ,16X,3F11.3)')
     2                                      QTR, HTR, CPTR, STR
            WRITE(6,'(7X,''  TOT.'',13X,F17.3,F11.4,2F11.4)')
     1                     ESCF+(HTOT-H298)/1000.D0,HTOT,CPTOT,STOT
         ENDIF
   80 CONTINUE
      WRITE(6,'(/3X,'' * NOTE: HEATS OF FORMATION ARE RELATIVE TO THE'',
     1/12X,'' ELEMENTS IN THEIR STANDARD STATE AT 298K'')')
      END

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


! {dg-do run }
!
! Tests the fix for PR82978 in which all the parameterized string
! lengths with the same value of parameter 'k' had the same value
! regardless of the value of 'l'. In this testcase, the length for
! 'l' = 5 was taken.
!
! Contributed by Fritz Reese  <foreese@gcc.gnu.org>
!
  implicit none

  type :: pdt_t(k, l)
    integer, kind :: k
    integer, len :: l
    character(kind=k,len=l) :: chr
    integer :: i(l)
  end type

  type(pdt_t(1, 4))   :: x1
  type(pdt_t(1, 5))   :: x2
  type(pdt_t(4, 5))   :: x3

  call test (x1, 4)
  call test (x2, 5)

! Kind tests appear because of problem identified in comment #!
! due to Dominque d'Humieres  <dominiq@lps.ens.fr>

  if (kind (x2%chr) .ne. 1) STOP 1
  if (kind (x3%chr) .ne. 4) STOP 2

contains

  subroutine test (x, i)
    type(pdt_t(1, *)) :: x
    integer :: i

    if (x%l .ne. i) STOP 3
    if (len(x%chr) .ne. i) STOP 4
    if (size(x%i,1) .ne. i) STOP 5
  end subroutine

end


text:
      SUBROUTINE HIGROW(I,IG,IM,ECC,SEMI,EMAX,EMIN,TG,EDAV,ZI,IQ)
*
*
*       Induced change of hierarchical binary.
*       --------------------------------------
*
      INCLUDE 'common6.h'
      COMMON/BINARY/  CM(4,MMAX),YREL(3,MMAX),ZREL(3,MMAX),
     &                HM(MMAX),UM(4,MMAX),UMDOT(4,MMAX),TMDIS(MMAX),
     &                NAMEM(MMAX),NAMEG(MMAX),KSTARM(MMAX),IFLAGM(MMAX)
      COMMON/MODES/  EB0(NTMAX),ZJ0(NTMAX),ECRIT(NTMAX),AR(NTMAX),
     &               BR(NTMAX),EOSC(4,NTMAX),EDEC(NTMAX),TOSC(NTMAX),
     &               RP(NTMAX),ES(NTMAX),ZM(2,NTMAX),IOSC(NTMAX),
     &               NAMEC(NTMAX)
      COMMON/SLOW0/  RANGE,ISLOW(10)
      common/tidal/  cq(2),ct(2),cgr,dedt
      common/rksave/  coeff,HOhat(3),e0,a0,hh,xmb
      REAL*8  A1(3),A2(3),XREL(3),VREL(3),EI(3),HI(3),HO(3),BHAT(3),
     &        UI(4),V(4),EVEC(3),XR0(3),VR0(3),EI0(3),WW(6),W(4)
      REAL*8  BODYI(2),WG(2)
      LOGICAL ICOLL
      DATA  WW  /2.119,3.113,8.175,3.742,4.953,9.413/
      SAVE  ICALL,ITIME,ITRY,NAMEI,DTPREV,TCHECK,ZFAC,PMIN1
      DATA  ICALL,ITRY,NAMEI,TCHECK  /0,0,0,0.0D0/
*
*
*       Copy KS variables to local scalars.
      DO 1 K = 1,4
          UI(K) = UM(K,IM)
          V(K) = UMDOT(K,IM)
    1 CONTINUE
*
*       Transform to physical variables and multiply by 4 (momentum formula).
      CALL KSPHYS(UI,V,XREL,VREL)
      DO 2 K = 1,3
          VREL(K) = 4.0*VREL(K)
    2 CONTINUE
*
*       Specify index for outer body (second KS component) & pair index.
      JCOMP = I + 1
      IPAIR = KVEC(I)
      DTSUM = 0.0
      ITIME = 0
      ICALL = 0
      ECC0 = ECC
      SEMI0 = SEMI
      a0 = SEMI
      HM0 = HM(IM)
      ZM3 = BODY(JCOMP)
      BODYI(1) = CM(1,IM)
      BODYI(2) = CM(2,IM)
*
*       Obtain inner & outer angular momentum by cyclic notation.
    4 A12 = 0.0
      A22 = 0.0
      A1A2 = 0.0
      RI2 = 0.0
      VI2 = 0.0
      RVI = 0.0
      DO 5 K = 1,3
          K1 = K + 1
          IF (K1.GT.3) K1 = 1
          K2 = K1 + 1
          IF (K2.GT.3) K2 = 1
          A1(K) = XREL(K1)*VREL(K2) - XREL(K2)*VREL(K1)
          A2(K) = (X(K1,JCOMP) - X(K1,I))*(XDOT(K2,JCOMP) - XDOT(K2,I))
     &          - (X(K2,JCOMP) - X(K2,I))*(XDOT(K1,JCOMP) - XDOT(K1,I))
          A12 = A12 + A1(K)**2
          A22 = A22 + A2(K)**2
          A1A2 = A1A2 + A1(K)*A2(K)
          RI2 = RI2 + XREL(K)**2
          VI2 = VI2 + VREL(K)**2
          RVI = RVI + XREL(K)*VREL(K)
    5 CONTINUE
*
*       Evaluate orbital parameters for outer orbit from KS elements.
      ZMB = BODY(I) + BODY(JCOMP)
      SEMI1 = -0.5*ZMB/H(IPAIR)
      ECC2 = (1.0 - R(IPAIR)/SEMI1)**2 + TDOT2(IPAIR)**2/(SEMI1*ZMB)
      ECC1 = SQRT(ECC2)
*       Determine inclination in radians.
      FAC = A1A2/SQRT(A12*A22)
      IF(FAC.LT.-1.D0.OR.FAC.GT.1.D0)THEN
         IF(FAC.LT.-1.05D0.OR.FAC.GT.1.05D0)THEN
            WRITE (6,9)  IPAIR, I, JCOMP, ECC1, SEMI1*SU, FAC
   9        FORMAT (' WARNING! HIGROW FAC: IP I J E A FAC ',
     &                                   3I6,F7.3,F8.1,F9.5)
         ENDIF
         FAC = MAX(FAC,-1.D0)
         FAC = MIN(FAC,1.D0)
      ENDIF
      ZI = ACOS(FAC)
      CC = (1.0 - ECC**2)*COS(ZI)**2
*
*       Construct the Runge-Lenz vector (Heggie & Rasio 1995, Eq.(5)).
      EI2 = 0.0
      DO 10 K = 1,3
          EI(K) = (VI2*XREL(K) - RVI*VREL(K))/BODY(I) -
     &                                                 XREL(K)/SQRT(RI2)
          EI2 = EI2 + EI(K)**2
          EVEC(K) = EI(K)
   10 CONTINUE
      EI2 = MIN(EI2,0.9999d0)
*
*       Define unit vectors for inner eccentricity and angular momenta.
      COSJ = 0.0
      SJSG = 0.0
      DO 15 K = 1,3
          EI(K) = EI(K)/SQRT(EI2)
          HI(K) = A1(K)/SQRT(A12)
          HO(K) = A2(K)/SQRT(A22)
          COSJ = COSJ + HI(K)*HO(K)
          SJSG = SJSG + EI(K)*HO(K)
   15 CONTINUE
*
*       Form unit vector BHAT and scalars AH & BH (Douglas Heggie, 10/9/96).
      AH = 0.0
      BH = 0.0
      DO 16 K = 1,3
          K1 = K + 1
          IF (K1.GT.3) K1 = 1
          K2 = K1 + 1
          IF (K2.GT.3) K2 = 1
          BHAT(K) = HI(K1)*EI(K2) - HI(K2)*EI(K1)
          AH = AH + EI(K)*HO(K)
          BH = BH + BHAT(K)*HO(K)
   16 CONTINUE
*
*       Evaluate the expressions A & Z.
      A = COSJ*SQRT(1.0 - EI2)
      Z = (1.0 - EI2)*(2.0 - COSJ**2) + 5.0*EI2*SJSG**2
*
*       Obtain maximum inner eccentricity (Douglas Heggie, Sept. 1995).
      Z2 = Z**2 + 25.0 + 16.0*A**4 - 10.0*Z - 20.0*A**2 - 8.0*A**2*Z
      EMAX = ONE6*(Z + 1.0 - 4.0*A**2 + SQRT(Z2))
      EMAX = MAX(EMAX,0.0001d0)
      EMAX = SQRT(EMAX)
*
*       Form minimum eccentricity (Douglas Heggie, Sept. 1996).
      AZ = A**2 + Z - 2.0
      IF (AZ.GE.0.0) THEN
          AZ1 = 1.0 + Z - 4.0*A**2
          EMIN2 = ONE6*(AZ1 - SQRT(AZ1**2 - 12.0*AZ))
      ELSE
          EMIN2 = 1.0 - 0.5*(A**2 + Z)
      END IF
      EMIN2 = MAX(EMIN2,0.0001D0)
      EMIN = SQRT(EMIN2)
*
*       Estimate eccentricity growth time-scale (N-body units).
      TK = TWOPI*SEMI*SQRT(SEMI/BODY(I))
      TK1 = TWOPI*ABS(SEMI1)*SQRT(ABS(SEMI1)/ZMB)
      TG = TK1**2*ZMB*(1.0 - ECC1**2)**1.5/(BODY(JCOMP)*TK)
*
*       Evaluate numerical precession factor (involves elliptic integral).
      CONST = PFAC(A,Z)
      CONST = CONST*4.0/(1.5*TWOPI*SQRT(6.0))
*
*       Convert growth time to units of 10**6 yrs.
      TG = CONST*TG*TSTAR
*
*       Form doubly averaged eccentricity derivative (Douglas Heggie 9/96).
      YFAC = 15.0*BODY(JCOMP)/(4.0*ZMB)*TWOPI*TK/TK1**2
      YFAC = YFAC*ECC*SQRT(1.0 - ECC**2)/(1.0 - ECC1**2)**(1.5)
      EDAV = YFAC*AH*BH
*
*       Skip small EMAX or large pericentre distance (TMDIS increased).
      QPERI = SEMI*(1.0 - ECC)
      PMIN = SEMI*(1.0 - EMAX)
      RM = MAX(RADIUS(I),RADIUS(IG),1.0D-20)
      IF (KSTARM(IM).GE.0.AND.EMAX.LT.0.9) THEN
          IQ = 1
          GO TO 40
      END IF
*
*       Delay Kozai cycle for EMAX > 0.9 and TC > 2000.
      IF (EMAX.GE.0.9) THEN
          CALL HICIRC(PMIN,EMAX,I,IG,BODYI,TG,TC,EC,EDT,W)
          IF (TC.GT.2000.OR.(ECC.LT.0.1.AND.EMAX.LT.0.95)) THEN
              IQ = 2
              GO TO 40
          END IF
      END IF
*
*       Check termination time for marginal stability criterion.
      IF (NAME(I).EQ.NAMEI.AND.(TIME+TOFF).GT.TCHECK) THEN
          WRITE (6,17)  NAME(I), ECC, EMAX, ALPH, ZFAC, ECC1, PMIN1,
     &                  PCRIT
   17     FORMAT (' ECCMOD UNSTAB    NM E EX IN YF E1 PM PC ',
     &                               I6,2F8.4,F7.1,F6.2,F7.3,1P,2E10.2)
          IQ = -4
          GO TO 40
      END IF
*
*       Ensure evaluation of EDT for large eccentricity and EDAV > 0.
      IF (ECC.GT.0.9.AND.EDAV.GT.0.0) THEN
          CALL HICIRC(QPERI,ECC,I,IG,BODYI,TG,TC,EC,EDT,W)
          IF (TC.LT.2000.0) THEN
              ICOLL = .TRUE.
          ELSE
              ICOLL = .FALSE.
              EDT = EDAV
          END IF
      ELSE
          ICOLL = .FALSE.
          EDT = EDAV
      END IF
*
*       Activate indicator on small QPERI to prevent collision before EMAX.
      IF (QPERI.LT.2.0*RM) ICOLL = .TRUE.
*
*       Check circularization time near turning point of normal binary.
      IF ((KSTARM(IM).GE.0.AND.EDAV.GT.0.0.AND.
     &   (DEDT.LT.0.0.OR.ECC.GT.ABS(EMAX-0.00005))).OR.ICOLL) THEN
          IF (.NOT.ICOLL) THEN
              CALL HICIRC(QPERI,ECC,I,IG,BODYI,TG,TC,EC,EDT,W)
          END IF
          IF (TC.LT.7.0D+04.AND.QPERI.LT.6.0*RM.AND.ITRY.LT.4.OR.
     &        TC.LT.1.0D+04.AND.QPERI.LT.5.0*RM.AND.ITRY.LT.8) THEN
              WRITE (6,18)  ECC, EMAX, EMIN, SEMI*SU, TG, TC, EDAV,
     &                      EDT, QPERI/RM
   18         FORMAT (' HICIRC TRY:    E EX EM A TG TC EDA EDT QP/R ',
     &                                 2F9.5,F7.3,F7.1,1P,5E9.1)
              ITRY = ITRY + 1
          END IF
          IF (ITRY.GT.8) ITRY = 0
*
*       Check termination for TC < 2000 to be consistent with routine TCIRC.
          IF (TC.GT.2000.0) GO TO 25
*
*       Evaluate chaos boundary parameters.
          EB = -0.5*CM(1,IM)*CM(2,IM)/a0
          CJ = CM(1,IM)*CM(2,IM)/BODY(I)*SQRT(BODY(I))
          ZJ = CJ*SQRT(QPERI*(1.0 + ECC))
          IC = NCHAOS + 1
          IDIS = 0
          ITRY = 0
*
*       Define oscillation period (dimensionless time).
          DO 20 K = 1,2
              IF (K.EQ.1) THEN
                  IK = I
              ELSE
                  IK = IG
              END IF
*       Specify polytropic index for each star (n = 3, 2 or 3/2).
              IF (KSTAR(IK).EQ.3.OR.KSTAR(IK).EQ.5) THEN
                  BODI = CM(K,IM)
                  CALL GIANT3(IK,BODI,WG,QG,ZN,QL)
                  W(K) = WG(1)
              ELSE
                  IP = 3
                  IF (KSTAR(IK).GE.3) IP = 2
                  IF (KSTAR(IK).EQ.4.OR.KSTAR(IK).EQ.6) IP = 3
                  IF (KSTAR(IK).EQ.0) IP = 1
                  W(K) = WW(IP)
              END IF
   20     CONTINUE
*
          CALL CHAOS0(QPERI,ECC,EB,ZJ,CM(1,IM),CM(2,IM),RADIUS(I),
     &                      RADIUS(IG),W,ECRIT(IC),AR(IC),BR(IC),IDIS)
*
*       Begin chaos/spiral stage if chaos boundary has been crossed.
          IF (IDIS.EQ.0) THEN
              WRITE (6,22)  TTOT, NAME(I), IC, ECC, EMAX, SEMI, QPERI,
     &                      EDAV, QPERI/RM
   22         FORMAT (' ECCMOD CHAOS    T NAM IC E EX A QP EDAV QP/R ',
     &                                  F9.2,I6,I4,2F9.5,1P,4E10.2)
              CALL FLUSH(6)
              IQ = -1
*       Activate chaos indicator for calling KSTIDE from RESET.
              KSTARM(IM) = -1
              GO TO 40
	  ELSE IF (IDIS.EQ.-1) THEN
              WRITE (6,23)  TTOT, NAME(I), IC, ECC, EMAX, SEMI, QPERI,
     &                      EDAV, QPERI/RM
   23         FORMAT (' ECCMOD SPIRAL    T NAM IC E EX A QP EDAV QP/R ',
     &                                   F9.2,I6,I4,2F9.5,1P,4E10.2)
              CALL FLUSH(6)
              IQ = -2
*       Note that KSTARM(IM) = -2 would be treated as existing spiral.
              KSTARM(IM) = -1
              GO TO 40
          ELSE IF (IDIS.EQ.1) THEN
*       Check collision condition to be sure.
              IF (QPERI.LT.RADIUS(I) + RADIUS(IG)) THEN
              WRITE (6,24)  TTOT, NAME(I), IC, ECC, SEMI, QPERI, EDAV,
     &                      QPERI/RM
   24         FORMAT (' ECCMOD DISRUPT    T NAM IC E A QP EDAV QP/R  ',
     &                                    F9.2,I6,I4,F8.4,1P,4E10.2)
                  IQ = -3
                  GO TO 40
              END IF
          END IF
      END IF
*
*       Form quantities for the outer orbit.
   25 RR = 0.0
      RV0 = 0.0
      V20 = 0.0
      DO 30 K = 1,3
          XR0(K) = X(K,JCOMP) - X(K,I)
          VR0(K) = XDOT(K,JCOMP) - XDOT(K,I)
          RV0 = RV0 + XR0(K)*VR0(K)
          V20 = V20 + VR0(K)**2
          RR = RR + XR0(K)**2
   30 CONTINUE
*
*       Construct the outer Runge-Lenz vector.
      DO 35 K = 1,3
          EI0(K) = (V20*XR0(K) - RV0*VR0(K))/ZMB - XR0(K)/SQRT(RR)
   35 CONTINUE
*
*       Specify total time interval (unperturbed or perturbed case).
      IF (LIST(1,I).EQ.0) THEN
          DT0 = TIME - T0(I)
*         DT = MIN(DT0,0.1*(1.0 - ECC)/(ABS(EDAV) + 1.0D-15))
      ELSE
*       Ensure that interval extends to next apocentre (including slow-down).
          IMOD = KSLOW(IPAIR)
          DT0 = 0.98*FLOAT(ISLOW(IMOD))*TK1
*         DT = MIN(DT0,0.1*(1.0 - ECC)*TG/TSTAR)
      END IF
*
*       Set partial interval and check remaining time.
      DT1 = MIN(DT0,0.1*TG*SQRT(1.0 - ECC**2)/TSTAR)
      IF (DTSUM + DT1.GT.DT0) DT1 = DT0 - DTSUM + 1.0D-12
*
*       Obtain quadrupole and tidal constants for each new binary.
      DTM = TIME - MAX(TEV0(I),TEV0(IG))
      IF (NAME(I).NE.NAMEI.OR.DTM.LE.DT0) THEN
          CALL QTIDES(I,IG,IM,SEMI,ECC)
          NAMEI = NAME(I)
          DEDT = 0.0
          DTPREV = 1.0
          NST = NSTAB(SEMI,SEMI1,ECC,ECC1,ZI,BODYI(1),BODYI(2),ZM3)
          IF (NST.EQ.0) THEN
              PCRIT = 0.98*QPERI*(1.0 - PERT)
              PCR = stability(BODYI(1),BODYI(2),ZM3,ECC,ECC1,ZI)
              PCR = PCR*SEMI
*       Specify reduced peri if old criterion < PCRIT/2 (avoids switching).
              IF (PCR.LT.0.5*PCRIT) THEN
                  PCRIT = 0.75*PCRIT
              END IF
              IF (PCRIT.LT.PCR.AND.PERT.LT.0.01.AND.
     &            ITIME.LT.20) THEN
                  ALPH = 360.0*ZI/TWOPI
                  FAIL = QPERI*(1-PERT) - PCR
                  WRITE (6,42)  TTOT, ALPH, ECC, ECC1, QPERI, FAIL, PERT
   42             FORMAT (' NEWSTAB    T INC EI EO QP FAIL PERT ',
     &                                 F7.1,F7.2,2F8.4,1P,3E10.2)
              END IF
          ELSE
              PCRIT = 1.01*QPERI
          END IF
          PMIN1 = SEMI1*(1.0 - ECC1)
          IF (PMIN1.LT.PCRIT) THEN
              NK = 1 + 10.0*ECC1/(1.0 - ECC1)
              TCHECK = TIME + TOFF + NK*TK1
          ELSE
              TCHECK = 1.0D+10
          END IF
      END IF
*
*       Determine time-step for integration from several criteria.
      ge=30*ECC+45*ECC**3+3.75*ECC**5
*       Apply factor of 2 correction for agreement with classical result.
      ge = 0.5*ge
      slr=SEMI*(1.0 - ECC**2)
      zq=ge/(slr**5*SEMI*sqrt(SEMI))
*       Correct for wrong eccentricity dependence.
      zq = zq*(1.0 - ECC**2)
      EDQ = zq*(cq(1) + cq(2))
      TQ = ECC/EDQ
      zgr = cgr/(slr*SEMI*SQRT(SEMI))
*       Note cgr is angular velocity (hence multiply by 2*pi; 12/03).
      TGR = TWOPI*ECC/zgr
*
*       Delay first time if apsidal motion or GR are important (ECC < 0.9).
      IF (ECC.LT.0.90.AND.ITIME.EQ.0) THEN
          IF (MIN(TQ,TGR)*TSTAR.LT.TG) THEN
              IQ = 2
              GO TO 40
          END IF
      END IF
*
*       Include check on tidal dissipation for large eccentricity.
      IF (KSTARM(IM).GE.0.AND.ECC.GT.0.9) THEN
          TT = ECC/ABS(EDT)
          TQ = MIN(TT,TQ)
*       Note that EDT from HICIRC agrees well with actual value from RKINT.
      END IF
*
*       Adopt harmonic mean of inner period and growth/quadrupole time.
      DT = 0.5*SQRT(TK*MIN(TG/TSTAR,TQ))
      IF (ECC.GT.0.98) DT = 0.1*DT
      DT = MIN(DT,0.001/(ABS(EDAV) + 1.0D-20))
*       Ensure conservative step first time and limit increase to factor 2.
      IF (ITIME.EQ.0) THEN
          DT = MIN(DT,1.0D-03*TG/TSTAR)
      ELSE
          DT = MIN(DT,2.0D0*DTPREV)
      END IF
*
*       Include safety check to avoid problems.
      ETRY = ECC + EDAV*DT1
      IF (ECC.GT.0.96.AND.ETRY.GT.0.98) THEN
          ETRY = MIN(ETRY,EMAX)
          QPTRY = SEMI*(1.0 - ETRY)
          CALL HICIRC(QPTRY,ETRY,I,IG,BODYI,TG,TC,EC,EDT,W)
          IF (TC.LT.10.0) THEN
              DT1 = 0.5*DT1
              DT = 0.5*DT
          END IF
      END IF
*
*       Ensure frequent termination check inside 5*RM to avoid overshooting.
      IF (QPERI.LT.5.0*RM.AND.EDAV.GT.0.0.AND.KSTARM(IM).GE.0) THEN
*       Restrict HIMOD interval so E + EDAV*DT1 gives new peri outside 4*RM.
          DTT = (QPERI - 4.0*RM)/(SEMI*EDAV)
          DT1 = MIN(DTT,DT1)
          DT1 = MAX(DT,DT1)
      END IF
*
*       Modify the inner orbit using averaging of de/dt & dh/dt.
      DTPREV = DT
      CALL HIMOD(EVEC,A1,ZM3,ZMB,EI0,A2,ICALL,DT1,DT,ECC,XREL,VREL)
*
      IF (e0.ge.1.0) THEN
          IQ = 1
          GO TO 40
      END IF
      SEMI = a0
      ECC = e0
      ITIME = ITIME + 1
      DTSUM = DTSUM + DT1
      IF (DTSUM.LT.DT0) GO TO 4
*
*       Update step counter and next time for checking.
      NEINT = NEINT + ITIME
      TMDIS(IM) = TIME + DTSUM
*
*       Transform back to KS variables after integration.
   40 IF (ITIME.GT.0) THEN
          CALL PHYSKS(XREL,VREL,UI,V)
          DO 45 K = 1,4
              UM(K,IM) = UI(K)
              UMDOT(K,IM) = 0.25*V(K)
              IF (K.LT.4) THEN
                  YREL(K,IM) = XREL(K)
                  ZREL(K,IM) = VREL(K)
              END IF
   45     CONTINUE
          HM(IM) = -0.5D0*BODY(I)/a0
          CALL HIRECT(IM)
          ZMU = CM(1,IM)*CM(2,IM)/BODY(I)
          DECORR = ZMU*(HM0 - HM(IM))
          EMERGE = EMERGE - DECORR
          ECOLL = ECOLL + DECORR
          EGRAV = EGRAV + DECORR
*       Define circularization if e0 < 0.002 (end of Kozai cycle).
          IF (e0.LT.0.002) THEN
              KM = KSTARM(IM)
              KSTARM(IM) = 10
              IF (KM.LT.0) NCIRC = NCIRC + 1
              TMDIS(IM) = 1.0D+10
              TEV(N+IPAIR) = MIN(TEV(I),TEV(IG)) - 2.0*STEPX
              TB = DAYS*SEMI*SQRT(SEMI/BODY(I))
              ALPH = ZI*360.0/TWOPI
              WRITE (6,48)  NAME(I), NAME(IG), KM, e0, a0*SU, TB, ALPH
   48         FORMAT (' END KOZAI    NM K* E A P INC ',
     &                               2I6,I4,F7.3,F6.1,2F7.1)
          END IF
*       Include rectification of active SPIRAL (suppressed).
          IC = 0
          DO 50 K = 1,NCHAOS
              IF (NAMEC(K).EQ.NZERO - NAMEM(IM)) IC = K
*             IF (NAMEC(K).EQ.NZERO + NAME(I)) IC = K
   50     CONTINUE
          IF (IC.EQ.-3) THEN
              RP(IC) = a0*(1.0 - ECC)
              ES(IC) = ECC
              TOSC(IC) = TIME
              DH = HM0 - HM(IM)
              WRITE (6,60)  TTOT, TMDIS(IM), ECC, RP(IC), DH
   60         FORMAT (' UPDATE:    T TMD E RP DH ',
     &                             2F10.2,F8.4,1P,2E10.2)
          END IF
      END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!    Implementation of the IEEE_FEATURES standard intrinsic module
!    Copyright (C) 2013-2017 Free Software Foundation, Inc.
!    Contributed by Francois-Xavier Coudert <fxcoudert@gcc.gnu.org>
!
! This file is part of the GNU Fortran runtime library (libgfortran).
!
! Libgfortran is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public
! License as published by the Free Software Foundation; either
! version 3 of the License, or (at your option) any later version.
!
! Libgfortran is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
!
! Under Section 7 of GPL version 3, you are granted additional
! permissions described in the GCC Runtime Library Exception, version
! 3.1, as published by the Free Software Foundation.
!
! You should have received a copy of the GNU General Public License and
! a copy of the GCC Runtime Library Exception along with this program;
! see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
! <http://www.gnu.org/licenses/>.  */

module IEEE_FEATURES

  implicit none
  private

  type, public :: IEEE_FEATURES_TYPE
    private
    integer :: hidden
  end type

  type(IEEE_FEATURES_TYPE), parameter, public :: &
    IEEE_DATATYPE       = IEEE_FEATURES_TYPE(0), &
    IEEE_DENORMAL       = IEEE_FEATURES_TYPE(1), &
    IEEE_DIVIDE         = IEEE_FEATURES_TYPE(2), &
    IEEE_HALTING        = IEEE_FEATURES_TYPE(3), &
    IEEE_INEXACT_FLAG   = IEEE_FEATURES_TYPE(4), &
    IEEE_INF            = IEEE_FEATURES_TYPE(5), &
    IEEE_INVALID_FLAG   = IEEE_FEATURES_TYPE(6), &
    IEEE_NAN            = IEEE_FEATURES_TYPE(7), &
    IEEE_ROUNDING       = IEEE_FEATURES_TYPE(8), &
    IEEE_SQRT           = IEEE_FEATURES_TYPE(9), &
    IEEE_UNDERFLOW_FLAG = IEEE_FEATURES_TYPE(10)

end module IEEE_FEATURES

max_stars_repo_name: peclatj/claw-compiler

max_issues_repo_name: peclatj/claw-compiler

max_forks_repo_name: peclatj/claw-compiler


text:
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c      written by the UFO converter
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      SUBROUTINE COUP2()

      IMPLICIT NONE

      DOUBLE PRECISION PI
      PARAMETER  (PI=3.141592653589793D0)
      INCLUDE 'input.inc'
      INCLUDE 'coupl.inc'

      INCLUDE 'model_functions.inc'
      UV_GQQB = MDL_COMPLEXI*(COND(DCMPLX(MDL_MB),DCMPLX(0.000000D
     $ +00),DCMPLX(-((MDL_G__EXP__2)/(2.000000D+00*4.800000D+01*PI*
     $ *2))*4.000000D+00*MDL_TF*REGLOG(DCMPLX(MDL_MB__EXP__2/MDL_MU_R__
     $ EXP__2)))))*G
      UV_GQQT = MDL_COMPLEXI*(COND(DCMPLX(MDL_MT),DCMPLX(0.000000D
     $ +00),DCMPLX(-((MDL_G__EXP__2)/(2.000000D+00*4.800000D+01*PI*
     $ *2))*4.000000D+00*MDL_TF*REGLOG(DCMPLX(MDL_MT__EXP__2/MDL_MU_R__
     $ EXP__2)))))*G
      UVWFCT_G_2 = COND(DCMPLX(MDL_MT),DCMPLX(0.000000D+00),DCMPLX(((MD
     $ L_G__EXP__2)/(2.000000D+00*4.800000D+01*PI**2))*4.000000D
     $ +00*MDL_TF*REGLOG(DCMPLX(MDL_MT__EXP__2/MDL_MU_R__EXP__2))))
      UVWFCT_G_1 = COND(DCMPLX(MDL_MB),DCMPLX(0.000000D+00),DCMPLX(((MD
     $ L_G__EXP__2)/(2.000000D+00*4.800000D+01*PI**2))*4.000000D
     $ +00*MDL_TF*REGLOG(DCMPLX(MDL_MB__EXP__2/MDL_MU_R__EXP__2))))
      R2_SXCW = ((MDL_CKM22*MDL_EE*MDL_COMPLEXI)/(MDL_SW*MDL_SQRT__2))
     $ *(-(MDL_G__EXP__2*(1.000000D+00+MDL_LHV)*(MDL_NCOL__EXP__2_M_1_0
     $ ))/(2.000000D+00*MDL_NCOL*1.600000D+01*PI**2))
      GC_4 = -G
      GC_5 = MDL_COMPLEXI*G
      R2_GQQ = -MDL_COMPLEXI*MDL_G__EXP__3/(1.600000D+01*PI**2)
     $ *((MDL_NCOL__EXP__2_M_1)/(2.000000D+00*MDL_NCOL))*(1.000000D+00
     $ +MDL_LHV)
      R2_QQQ = MDL_LHV*MDL_COMPLEXI*MDL_G__EXP__2*(MDL_NCOL__EXP__2_M_1
     $ )/(3.200000D+01*PI**2*MDL_NCOL)
      UV_GQQG_1EPS = MDL_COMPLEXI*(-((MDL_G__EXP__2)/(2.000000D
     $ +00*4.800000D+01*PI**2))*1.100000D+01*MDL_CA)*G
      UV_GQQQ_1EPS = MDL_COMPLEXI*(((MDL_G__EXP__2)/(2.000000D
     $ +00*4.800000D+01*PI**2))*4.000000D+00*MDL_TF)*G
      UVWFCT_G_2_1EPS = COND(DCMPLX(MDL_MT),DCMPLX(0.000000D+00)
     $ ,DCMPLX(-((MDL_G__EXP__2)/(2.000000D+00*4.800000D+01*PI**2))
     $ *4.000000D+00*MDL_TF))
      END


text:
C Copyright(C) 2009-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     * 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=======================================================================
      SUBROUTINE CMDWHE (VERB, INLINE,
     &   IFLD, INTYP, CFIELD, IFIELD, RFIELD,
     &   XN, YN, ZN, XNN, YNN, ZNN, HIDENP, XE, YE, ZE, MAPEL, MAPND, *)
C=======================================================================

C   --*** CMDWHE *** (MESH) Process locate commands
C   --   Written by Amy Gilkey - revised 02/05/88
C   --
C   --Parameters:
C   --   VERB - IN/OUT - the verbs for the SHOW command
C   --   INLINE - IN/OUT - the parsed input line for the log file
C   --   IFLD, INTYP, CFIELD, IFIELD, RFIELD - IN/OUT - the free-field
C   --      reader index and fields
C   --   XN, YN, ZN - IN - the coordinates of the original nodes
C   --   XNN, YNN, ZNN - IN - the coordinates of the displayed nodes
C   --   HIDENP(i) - IN - true iff node i is visible,
C   --      displayed mesh coordinates are defined if not HIDENP(i)
C   --   XE, YE, ZE - IN - the coordinates of the element centroids
C   --
C   --Common Variables:
C   --   Uses NUMNP of /DBNUMS/
C   --   Uses IS3DIM of /D3NUMS/

      PARAMETER (KLFT=1, KRGT=2, KBOT=3, KTOP=4, KNEA=5, KFAR=6)

      include 'params.blk'
      include 'dbnums.blk'
      include 'd3nums.blk'
      include 'rotopt.blk'

      CHARACTER*(*) VERB
      CHARACTER*(*) INLINE
      INTEGER     INTYP(*)
      CHARACTER*(*) CFIELD(*)
      INTEGER     IFIELD(*)
      REAL        RFIELD(*)
      REAL XN(NUMNP), YN(NUMNP), ZN(NUMNP)
      REAL XNN(NUMNP), YNN(NUMNP), ZNN(NUMNP)
      LOGICAL HIDENP(NUMNP)
      REAL XE(NUMEL), YE(NUMEL), ZE(NUMEL)
      INTEGER MAPEL(*), MAPND(*)

      CHARACTER*(MXSTLN) WORD
      LOGICAL ISON
      REAL RNUM(3)
      CHARACTER*20 RSTR(3)
      LOGICAL WINDOW
      LOGICAL FFMATC, MATSTR
      LOGICAL LDUM

      IF ((VERB .EQ. 'WHAT') .OR. (VERB .EQ. 'WHAT3')) THEN
         INLINE = ' '
         WINDOW = (VERB .EQ. 'WHAT')
         VERB = ' '

         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)

         IF (MATSTR (WORD, 'NODE', 1)) THEN
            ISON = FFMATC (IFLD, INTYP, CFIELD, 'CURSOR', 1)
            IF (WINDOW) THEN
               CALL PICKNP ('object location', ISON, WINDOW,
     &            NDIM, NUMNP, XNN, YNN, ZNN, HIDENP,
     &            .TRUE., IFLD, INTYP, RFIELD,
     &            INP, *100)
            ELSE
               CALL PICKNP ('object location', ISON, WINDOW,
     &            NDIM, NUMNP, XN, YN, ZN, LDUM,
     &            .TRUE., IFLD, INTYP, RFIELD,
     &            INP, *100)
            END IF

            CALL INTSTR (1, 0, MAPND(INP), RSTR(1), L)
            WRITE (*, 10000) 'Nearest node is node ', RSTR(1)(:L)
            RNUM(1) = XN(INP)
            RNUM(2) = YN(INP)
            IF (IS3DIM) RNUM(3) = ZN(INP)
            CALL NUMSTR (3, 4, RNUM(1), RSTR, LSTR)
            WRITE (*, 10000) 'Undeformed object location is',
     &         (' ', RSTR(I)(:LSTR), I=1,3)

          ELSE IF (MATSTR (WORD, 'ELEMENT', 1)) THEN
            IF (FFMATC (IFLD, INTYP, CFIELD, 'CURSOR', 1)) THEN
              CALL PRTERR ('WARNING',
     *          'CURSOR option not supported for elements')
              GO TO 100
            ENDIF
            ISON = .FALSE.
            CALL PICKNP ('object location', ISON, .FALSE.,
     &        NDIM, NUMEL, XE, YE, ZE, LDUM,
     &        .TRUE., IFLD, INTYP, RFIELD,
     &        INP, *100)

            CALL INTSTR (1, 0, MAPEL(INP), RSTR(1), L)
            WRITE (*, 10000) 'Nearest object is element ', RSTR(1)(:L)
            RNUM(1) = XE(INP)
            RNUM(2) = YE(INP)
            IF (IS3DIM) RNUM(3) = ZE(INP)
            CALL NUMSTR (3, 4, RNUM(1), RSTR, LSTR)
            WRITE (*, 10000) 'Undeformed object location is',
     &         (' ', RSTR(I)(:LSTR), I=1,3)

         ELSE
            CALL PRTERR ('CMDERR', 'Expected "NODE" or "ELEMENT"')
            GOTO 100
         END IF

      ELSE IF (VERB .EQ. 'WHERE') THEN
         INLINE = ' '
         VERB = ' '
         CALL FFCHAR (IFLD, INTYP, CFIELD, ' ', WORD)

         IF (MATSTR (WORD, 'NODE', 1)) THEN
            CALL FFINTG (IFLD, INTYP, IFIELD,
     &         'node number', 0, NOD, *100)
            INP = locint(NOD, NUMNP, MAPND)
            IF ((INP .LT. 1) .OR. (INP .GT. NUMNP)) THEN
               CALL PRTERR ('CMDERR', 'Invalid node id')
               GOTO 100
            END IF

            RNUM(1) = XN(INP)
            RNUM(2) = YN(INP)
            IF (IS3DIM) RNUM(3) = ZN(INP)
            CALL NUMSTR (3, 4, RNUM(1), RSTR, LSTR)
            WRITE (*, 10000) 'Undeformed object location is',
     &         (' ', RSTR(I)(:LSTR), I=1,NDIM)
            IF (IS3DIM .AND. (.NOT. HIDENP(INP))) THEN
               RNUM(1) = XNN(INP)
               RNUM(2) = YNN(INP)
               RNUM(3) = ZNN(INP)
               CALL NUMSTR (2, 4, RNUM(1), RSTR, LSTR)
               WRITE (*, 10000) 'Displayed window location is',
     &            (' ', RSTR(I)(:LSTR), I=1,2)
            END IF

          ELSE IF (MATSTR (WORD, 'ELEMENT', 1)) THEN
            CALL FFINTG (IFLD, INTYP, IFIELD,
     &         'element number', 0, IEL, *100)
            INP = locint(IEL, NUMEL, MAPEL)
            IF ((INP .LT. 1) .OR. (INP .GT. NUMEL)) THEN
               CALL PRTERR ('CMDERR', 'Invalid element number')
               GOTO 100
            END IF

            RNUM(1) = XE(INP)
            RNUM(2) = YE(INP)
            IF (IS3DIM) RNUM(3) = ZE(INP)
            CALL NUMSTR (3, 4, RNUM(1), RSTR, LSTR)
            WRITE (*, 10000) 'Undeformed object location is',
     &         (' ', RSTR(I)(:LSTR), I=1,NDIM)

         ELSE
            CALL PRTERR ('CMDERR', 'Expected "NODE" or "ELEMENT"')
            GOTO 100
         END IF

      END IF

      RETURN

  100 CONTINUE
      RETURN 1
10000  FORMAT (1X, 10A)
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
module gsi_thinning
!$$$ subprogram documenation block
!
! abstract:  This module contains code which may be used to selectively
!            thin satellite data.
!
! program history log:
!   2006-10-28 Jianjun Xu, developed from GSI
!   2007-03-30 Zhiquan Liu, modify and add comments
!
!
! Subroutines Included:
!   sub makegrids      - set up thinning grids
!   sub map2tgrids     - map observation to location on thinning grid
!   sub destroygrids   - deallocate thinning grid arrays
!   sub destroy_sfc    - deallocate full horizontal surface arrays
!
! Variable Definitions:
!   def mlat           - number of latitudes in thinning grid
!   def mlon           - number of longitudes in thinning grid
!   def maxthin        - maximum number of obs to retain in thinning grid box
!   def itxmax         - total number of points in thinning grid
!   def istart_val     - starting location on thinning grid for superobs (not used)
!   def icount         - observation counter
!   def ibest_obs      - index of "best" quality obs in thinning grid box
!   def glat           - latitudes on thinning grid
!   def glon           - longitudes on thinning grid
!   def hll            - (i,j) index of point on thinning grid
!   def score_crit     - "best" quality obs score in thinning grid box
!
!$$$ end documentation block
  use da_control, only: kms, kme
  use da_reporting, only : message, da_error
  use gsi_kinds, only: r_kind,i_kind
  use gsi_constants, only: deg2rad,rearth_equator,zero,two,pi,half,one,quarter,&
       rad2deg
  implicit none

  real(r_kind),parameter:: r90   = 90.0_r_kind
  real(r_kind),parameter:: r360  = 360.0_r_kind
  real(r_kind),parameter:: r999  = 999.0_r_kind
  real(r_kind),parameter:: r1000 = 1000.0_r_kind

  ! lat/lon range inside tile
  real(r_kind) rlat_min,rlat_max,rlon_min,rlon_max,dlat_grid,dlon_grid

  type thinning_type
! mlat: lat #, mlonx: max lon #, itxmax: grid box #
    integer(i_kind) mlat,maxthin,itxmax,dthin,mlonx,mlony
    integer(i_kind),dimension(0:51):: istart_val

! mlon(mlat): lon # in each lat
    integer(i_kind),allocatable,dimension(:):: mlon,icount,ibest_obs
    integer(i_kind),allocatable,dimension(:,:):: isli

! glat(mlat): lat #, glon(mlat,mlonx), hll(mlat,mlonx)
    integer(i_kind),allocatable,dimension(:,:) :: hll
    integer(i_kind),allocatable,dimension(:,:,:) :: hll_3d
    real(r_kind),allocatable,dimension(:)   :: glat
    real(r_kind),allocatable,dimension(:,:) :: glon,sli,sno
    real(r_kind),allocatable,dimension(:)   :: score_crit
  end type thinning_type

  type(thinning_type), allocatable  :: thinning_grid(:,:)
  type(thinning_type), allocatable  :: thinning_grid_conv(:)

contains

  subroutine makegrids (n,rmesh,ifgat)
! compute dimention of thinning box
! output (mlat,mlonx,istart_val)
    implicit none

    integer(i_kind), intent(in) :: n  ! sensor index
    real(r_kind), intent(in) :: rmesh ! thinning box size
    integer(i_kind), intent(in) :: ifgat

    logical odd
    integer(i_kind) i,ii,j,k,nlat,nlon
    integer(i_kind) icnt,mlonj
    real(r_kind) delonx,delat,dgv,dx,dy
    real(r_kind) twopi,halfpi,dlon_g,dlat_g,dlon_e,dlat_e
    real(r_kind) factor,factors,delon
    real(r_kind) rkm2dg,glatm,glatx

!   Initialize variables, set constants
      thinning_grid(n,ifgat)%dthin = 1
      thinning_grid(n,ifgat)%maxthin=thinning_grid(n,ifgat)%dthin

      thinning_grid(n,ifgat)%istart_val=0
      twopi  = two*pi
      halfpi = pi*half
      rkm2dg = r360/(twopi*rearth_equator)*r1000

       dx    = rmesh*rkm2dg
       dy    = dx
       thinning_grid(n,ifgat)%mlat  = dlat_grid/dy + half
       thinning_grid(n,ifgat)%mlonx = dlon_grid/dx + half
       delat = dlat_grid/thinning_grid(n,ifgat)%mlat
       delonx= dlon_grid/thinning_grid(n,ifgat)%mlonx
       dgv   = delat*half

       thinning_grid(n,ifgat)%mlat=max(2,thinning_grid(n,ifgat)%mlat)
       thinning_grid(n,ifgat)%mlonx=max(2,thinning_grid(n,ifgat)%mlonx)

      do ii=1,thinning_grid(n,ifgat)%maxthin
       thinning_grid(n,ifgat)%istart_val(ii+1)=thinning_grid(n,ifgat)%istart_val(ii)
          icnt=0
          do j = 1,thinning_grid(n,ifgat)%mlat
             glatx = rlat_min + (j-1)*delat
             glatx = glatx*deg2rad
             glatm = glatx + dgv*deg2rad
             factor = abs(cos(abs(glatm)))
             mlonj = nint(thinning_grid(n,ifgat)%mlonx*factor)
             mlonj = max(2,mlonj)
             do i = 1,mlonj
                icnt=icnt+1
                thinning_grid(n,ifgat)%istart_val(ii+1)=thinning_grid(n,ifgat)%istart_val(ii+1)+1
             enddo
          enddo
      end do

! making thinning box
! output: mlon(mlat),glat(mlat),glon(mlonx,mlat),hll(mlonx,mlat)

    allocate(thinning_grid(n,ifgat)%mlon(thinning_grid(n,ifgat)%mlat), &
             thinning_grid(n,ifgat)%glat(thinning_grid(n,ifgat)%mlat), &
             thinning_grid(n,ifgat)%glon(thinning_grid(n,ifgat)%mlonx,thinning_grid(n,ifgat)%mlat), &
             thinning_grid(n,ifgat)%hll(thinning_grid(n,ifgat)%mlonx,thinning_grid(n,ifgat)%mlat))

!   Set up thinning grid lon & lat.  The lon & lat represent the location of the
!   lower left corner of the thinning grid box.

       thinning_grid(n,ifgat)%itxmax=0
      do j = 1,thinning_grid(n,ifgat)%mlat
       thinning_grid(n,ifgat)%glat(j) = rlat_min + (j-1)*delat
       thinning_grid(n,ifgat)%glat(j) = thinning_grid(n,ifgat)%glat(j)*deg2rad
       glatm = thinning_grid(n,ifgat)%glat(j) + dgv*deg2rad

       factor = abs(cos(abs(glatm)))
       mlonj  = nint(thinning_grid(n,ifgat)%mlonx*factor)
       thinning_grid(n,ifgat)%mlon(j) = max(2,mlonj)
       delon = dlon_grid/thinning_grid(n,ifgat)%mlon(j)

       thinning_grid(n,ifgat)%glat(j) = min(max(-halfpi,thinning_grid(n,ifgat)%glat(j)),halfpi)
       do i = 1,thinning_grid(n,ifgat)%mlon(j)
          thinning_grid(n,ifgat)%itxmax=thinning_grid(n,ifgat)%itxmax+1
          thinning_grid(n,ifgat)%hll(i,j)=thinning_grid(n,ifgat)%itxmax
          thinning_grid(n,ifgat)%glon(i,j) = rlon_min + (i-1)*delon
          thinning_grid(n,ifgat)%glon(i,j) = thinning_grid(n,ifgat)%glon(i,j)*deg2rad
          thinning_grid(n,ifgat)%glon(i,j) = min(max(zero,thinning_grid(n,ifgat)%glon(i,j)),twopi)
       enddo
       !write(6,'(f10.5,i8,2i10)') glat(j)*rad2deg, mlon(j),hll(1,j),hll(mlon(j),j)
       !write(6,'(10f8.3)')   (glon(i,j)*rad2deg,i=1,mlon(j))

    end do

!   Allocate  and initialize arrays
    allocate(thinning_grid(n,ifgat)%icount(thinning_grid(n,ifgat)%itxmax))
    allocate(thinning_grid(n,ifgat)%ibest_obs(thinning_grid(n,ifgat)%itxmax))
    allocate(thinning_grid(n,ifgat)%score_crit(thinning_grid(n,ifgat)%itxmax))

    do j=1,thinning_grid(n,ifgat)%itxmax
       thinning_grid(n,ifgat)%icount(j)     = 0
       thinning_grid(n,ifgat)%ibest_obs(j)  = 0
       thinning_grid(n,ifgat)%score_crit(j) = 9.99e6_r_kind
    end do

    return
  end subroutine makegrids

  subroutine make3grids (n,rmesh, thin_3d)
! compute dimention of thinning box
! output (mlat,mlonx,istart_val)
    implicit none

    integer(i_kind), intent(in) :: n  ! sensor index
    real(r_kind), intent(in) :: rmesh ! thinning box size
    logical, optional, intent(in) :: thin_3d

    integer(i_kind) i,ii,j,k,nlat,nlon
    integer(i_kind) icnt,mlonj
    real(r_kind) delonx,delat,dgv,dx,dy
    real(r_kind) twopi,halfpi,dlon_g,dlat_g,dlon_e,dlat_e
    real(r_kind) factor,factors,delon
    real(r_kind) rkm2dg,glatm,glatx

!   Initialize variables, set constants
      thinning_grid_conv(n)%dthin = 1
      thinning_grid_conv(n)%maxthin=thinning_grid_conv(n)%dthin

      thinning_grid_conv(n)%istart_val=0
      twopi  = two*pi
      halfpi = pi*half
      rkm2dg = r360/(twopi*rearth_equator)*r1000

       dx    = rmesh*rkm2dg
       dy    = dx
       thinning_grid_conv(n)%mlat  = dlat_grid/dy + half
       thinning_grid_conv(n)%mlonx = dlon_grid/dx + half
       delat = dlat_grid/thinning_grid_conv(n)%mlat
       delonx= dlon_grid/thinning_grid_conv(n)%mlonx
       dgv   = delat*half

       thinning_grid_conv(n)%mlat=max(2,thinning_grid_conv(n)%mlat)
       thinning_grid_conv(n)%mlonx=max(2,thinning_grid_conv(n)%mlonx)

      do ii=1,thinning_grid_conv(n)%maxthin
       thinning_grid_conv(n)%istart_val(ii+1)=thinning_grid_conv(n)%istart_val(ii)
          icnt=0
          do j = 1,thinning_grid_conv(n)%mlat
             glatx = rlat_min + (j-1)*delat
             glatx = glatx*deg2rad
             glatm = glatx + dgv*deg2rad
             factor = abs(cos(abs(glatm)))
             mlonj = nint(thinning_grid_conv(n)%mlonx*factor)
             mlonj = max(2,mlonj)
             do i = 1,mlonj
                icnt=icnt+1
                thinning_grid_conv(n)%istart_val(ii+1)=thinning_grid_conv(n)%istart_val(ii+1)+1
             enddo
          enddo
      end do

! making thinning box
! output: mlon(mlat),glat(mlat),glon(mlonx,mlat),hll(mlonx,mlat)

    allocate(thinning_grid_conv(n)%mlon(thinning_grid_conv(n)%mlat), &
             thinning_grid_conv(n)%glat(thinning_grid_conv(n)%mlat), &
             thinning_grid_conv(n)%glon(thinning_grid_conv(n)%mlonx,thinning_grid_conv(n)%mlat))
     if( .not. present(thin_3d)) then
          allocate(thinning_grid_conv(n)%hll(thinning_grid_conv(n)%mlonx,thinning_grid_conv(n)%mlat))
     else
       allocate(thinning_grid_conv(n)%hll_3d(thinning_grid_conv(n)%mlonx,thinning_grid_conv(n)%mlat, kms:kme))
     end if

!   Set up thinning grid lon & lat.  The lon & lat represent the location of the
!   lower left corner of the thinning grid box.

       thinning_grid_conv(n)%itxmax=0

      do j = 1,thinning_grid_conv(n)%mlat
       thinning_grid_conv(n)%glat(j) = rlat_min + (j-1)*delat
       thinning_grid_conv(n)%glat(j) = thinning_grid_conv(n)%glat(j)*deg2rad
       glatm = thinning_grid_conv(n)%glat(j) + dgv*deg2rad

       factor = abs(cos(abs(glatm)))
       mlonj  = nint(thinning_grid_conv(n)%mlonx*factor)
       thinning_grid_conv(n)%mlon(j) = max(2,mlonj)
       delon = dlon_grid/thinning_grid_conv(n)%mlon(j)

       thinning_grid_conv(n)%glat(j) = min(max(-halfpi,thinning_grid_conv(n)%glat(j)),halfpi)
       do i = 1,thinning_grid_conv(n)%mlon(j)
          thinning_grid_conv(n)%glon(i,j) = rlon_min + (i-1)*delon
          thinning_grid_conv(n)%glon(i,j) = thinning_grid_conv(n)%glon(i,j)*deg2rad
          thinning_grid_conv(n)%glon(i,j) = min(max(zero,thinning_grid_conv(n)%glon(i,j)),twopi)
          if( .not. present(thin_3d)) then
             thinning_grid_conv(n)%itxmax=thinning_grid_conv(n)%itxmax+1
             thinning_grid_conv(n)%hll(i,j)=thinning_grid_conv(n)%itxmax
          else
             do k=kms, kme
                thinning_grid_conv(n)%itxmax=thinning_grid_conv(n)%itxmax+1
                thinning_grid_conv(n)%hll_3d(i,j,k)=thinning_grid_conv(n)%itxmax
             end do
          end if
       enddo
       !write(6,'(f10.5,i8,2i10)') glat(j)*rad2deg, mlon(j),hll(1,j),hll(mlon(j),j)
       !write(6,'(10f8.3)')   (glon(i,j)*rad2deg,i=1,mlon(j))

    end do

!   Allocate  and initialize arrays
    allocate(thinning_grid_conv(n)%icount(thinning_grid_conv(n)%itxmax))
    allocate(thinning_grid_conv(n)%ibest_obs(thinning_grid_conv(n)%itxmax))
    allocate(thinning_grid_conv(n)%score_crit(thinning_grid_conv(n)%itxmax))

    do j=1,thinning_grid_conv(n)%itxmax
       thinning_grid_conv(n)%icount(j)     = 0
       thinning_grid_conv(n)%ibest_obs(j)  = 0
       thinning_grid_conv(n)%score_crit(j) = 9.99e6_r_kind
    end do

    return
  end subroutine make3grids

  subroutine map2grids(n,ifgat,dlat_earth,dlon_earth,crit1,iobs,itx,ithin,itt,iobsout,iuse)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    map2grids
!     prgmmr:    treadon     org: np23                date: 2002-10-17
!
! abstract:  This routine maps observations to the thinning grid.
!
! program history log:
!   2002-10-17  treadon
!   2004-06-22  treadon - update documentation
!   2004-07-23  derber - modify code to thin obs as read in
!   2004-12-08  li, xu - fix bug --> set iuse=.true. when use_all=.true.
!   2005-10-14  treadon - variable name change (dlat0,dlon0) --> d*_earth
!   2006-03-25  kistler - define iobsout for the case use_all=.true.
!
!   input argument list:
!         n      - sensor index
!     ifgat      - number of time slots for FGAT and 4DVAR
!     dlat_earth - earth relative observation latitude (radians)
!     dlon_earth - earth relative observation longitude (radians)
!     crit1      - quality indicator for observation (smaller = better)
!     ithin      - number of obs to retain per thinning grid box
!
!   output argument list:
!     iobs  - observation counter
!     itx   - combined (i,j) index of observation on thinning grid
!     itt   - superobs thinning counter
!     iobsout- location for observation to be put
!     iuse  - .true. if observation should be used
!
    implicit none
    logical, intent(out) :: iuse
    integer(i_kind), intent(out) :: itt,itx
    integer(i_kind), intent(in)  :: ithin,n,ifgat
    integer(i_kind), intent(inout) :: iobs, iobsout
    real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1

    integer(i_kind) :: ix,iy
    real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy
    real(r_kind) dist1,crit

!   Compute (i,j) indices of coarse mesh grid (grid number 1) which
!   contains the current observation.
    dlat1=dlat_earth
    dlon1=dlon_earth

    call grdcrd(dlat1,1,thinning_grid(n,ifgat)%glat,thinning_grid(n,ifgat)%mlat,1)
    iy=int(dlat1)
    dy=dlat1-iy
    iy=max(1,min(iy,thinning_grid(n,ifgat)%mlat))

    call grdcrd(dlon1,1,thinning_grid(n,ifgat)%glon(1,iy),thinning_grid(n,ifgat)%mlon(iy),1)
    ix=int(dlon1)
    dx=dlon1-ix
    ix=max(1,min(ix,thinning_grid(n,ifgat)%mlon(iy)))

    dxx=half-min(dx,one-dx)
    dyy=half-min(dy,one-dy)
    dist1=dxx*dxx+dyy*dyy+half
    itx=thinning_grid(n,ifgat)%hll(ix,iy)
    itt=thinning_grid(n,ifgat)%istart_val(ithin)+itx
    if(ithin == 0) itt=0

!   Increment obs counter on coarse mesh grid.  Also accumulate observation
!   score and distance functions

    thinning_grid(n,ifgat)%icount(itx)=thinning_grid(n,ifgat)%icount(itx)+1
!   dist1=one - quarter*(dista + distb)  !dist1 is min at grid box center and
                                    !ranges from 1 (at corners)to
                                    !.5 (at center of box)
    crit=crit1*dist1
    iuse=.false.

    if(thinning_grid(n,ifgat)%icount(itx) == 1)then

!   Increment obs counter

      iuse=.true.
      iobs=iobs+1
      thinning_grid(n,ifgat)%score_crit(itx)= crit
      thinning_grid(n,ifgat)%ibest_obs(itx) = iobs
      iobsout=iobs

    end if
    if(crit < thinning_grid(n,ifgat)%score_crit(itx) .and. thinning_grid(n,ifgat)%icount(itx) > 1)then
      iuse=.true.
      iobs=iobs+1
      thinning_grid(n,ifgat)%score_crit(itx)= crit
      thinning_grid(n,ifgat)%ibest_obs(itx)=iobs
      iobsout=iobs
    end if

    return
  end subroutine map2grids

  subroutine map2grids_conv(n,dlat_earth,dlon_earth,crit1,iobs,itx,ithin,itt,iobsout,iuse, zk, thin_3d)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    map2grids
!     prgmmr:    treadon     org: np23                date: 2002-10-17
!
! abstract:  This routine maps observations to the thinning grid.
!
! program history log:
!   2002-10-17  treadon
!   2004-06-22  treadon - update documentation
!   2004-07-23  derber - modify code to thin obs as read in
!   2004-12-08  li, xu - fix bug --> set iuse=.true. when use_all=.true.
!   2005-10-14  treadon - variable name change (dlat0,dlon0) --> d*_earth
!   2006-03-25  kistler - define iobsout for the case use_all=.true.
!
!   input argument list:
!         n      - sensor index
!     dlat_earth - earth relative observation latitude (radians)
!     dlon_earth - earth relative observation longitude (radians)
!     crit1      - quality indicator for observation (smaller = better)
!     ithin      - number of obs to retain per thinning grid box
!
!   output argument list:
!     iobs  - observation counter
!     itx   - combined (i,j) index of observation on thinning grid
!     itt   - superobs thinning counter
!     iobsout- location for observation to be put
!     iuse  - .true. if observation should be used
!
    implicit none
    logical, intent(out) :: iuse
    integer(i_kind), intent(out) :: itt,itx
    integer(i_kind), intent(in)  :: ithin,n
    integer(i_kind), intent(inout) :: iobs, iobsout
    real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1
    real(r_kind), optional, intent(in) :: zk
    logical, optional, intent(in) :: thin_3d

    integer(i_kind) :: ix,iy,iz
    real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy,dz,dzz
    real(r_kind) dist1,crit,dist2
    integer(i_kind) :: model_lev

    model_lev=kme-kms+1

!   Compute (i,j) indices of coarse mesh grid (grid number 1) which
!   contains the current observation.
    dlat1=dlat_earth
    dlon1=dlon_earth

    call grdcrd(dlat1,1,thinning_grid_conv(n)%glat,thinning_grid_conv(n)%mlat,1)
    iy=int(dlat1)
    dy=dlat1-iy
    iy=max(1,min(iy,thinning_grid_conv(n)%mlat))

    call grdcrd(dlon1,1,thinning_grid_conv(n)%glon(1,iy),thinning_grid_conv(n)%mlon(iy),1)
    ix=int(dlon1)
    dx=dlon1-ix
    ix=max(1,min(ix,thinning_grid_conv(n)%mlon(iy)))

    dxx=half-min(dx,one-dx)
    dyy=half-min(dy,one-dy)
    dist1=dxx*dxx+dyy*dyy+half
    dist2=1.
    if( present (thin_3d)) then
       if( present(zk)) then
          iz=int(zk)
          dz=zk-iz
          iz=max(1,min(iz,model_lev))
          dzz=half-min(dz,one-dz)
          dist2=dzz*dzz+half
          itx=thinning_grid_conv(n)%hll_3d(ix,iy,iz)
       else
          write(unit=message(1),fmt='(A)')' For 3D thinning zk is required'
          call da_error(__FILE__,__LINE__,message(1:1))
       end if
    else
       itx=thinning_grid_conv(n)%hll(ix,iy)
    end if
    itt=thinning_grid_conv(n)%istart_val(ithin)+itx
    if(ithin == 0) itt=0

!   Increment obs counter on coarse mesh grid.  Also accumulate observation
!   score and distance functions

    thinning_grid_conv(n)%icount(itx)=thinning_grid_conv(n)%icount(itx)+1
!   dist1=one - quarter*(dista + distb)  !dist1 is min at grid box center and
                                    !ranges from 1 (at corners)to
                                    !.5 (at center of box)
    crit=crit1*dist1*dist2
    iuse=.false.

    if(thinning_grid_conv(n)%icount(itx) == 1)then

!   Increment obs counter

      iuse=.true.
      iobs=iobs+1
      thinning_grid_conv(n)%score_crit(itx)= crit
      thinning_grid_conv(n)%ibest_obs(itx) = iobs
      iobsout=iobs

    end if
    if(crit < thinning_grid_conv(n)%score_crit(itx) .and. thinning_grid_conv(n)%icount(itx) > 1)then
      iuse=.true.
      thinning_grid_conv(n)%score_crit(itx)= crit
      iobsout = thinning_grid_conv(n)%ibest_obs(itx)
    end if

    return
  end subroutine map2grids_conv

  subroutine map2tgrid(n,ifgat,dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse)
!   input argument list:
!     dlat_earth - earth relative observation latitude (radians)
!     dlon_earth - earth relative observation longitude (radians)
!     crit1      - quality indicator for observation (smaller = better)
!     ithin      - number of obs to retain per thinning grid box
!     ifgat      - number of time slots for FGAT and 4DVAR
!   output argument list:
!     iobs  - observation counter
!     itx   - combined (i,j) index of observation on thinning grid
!     itt   - superobs thinning counter
!     iobsout- location for observation to be put
!     iuse  - .true. if observation should be used
!
    implicit none
    logical,intent(out):: iuse
    integer(i_kind),intent(in):: ithin,n,ifgat
    integer(i_kind),intent(out):: itt,itx
    real(r_kind),intent(in):: dlat_earth,dlon_earth,crit1
    real(r_kind),intent(out):: dist1

    integer(i_kind) ix,iy
    real(r_kind) dlat1,dlon1,dx,dy,dxx,dyy

!   Compute (i,j) indices of coarse mesh grid (grid number 1) which
!   contains the current observation.
    dlat1=dlat_earth
    dlon1=dlon_earth

    call grdcrd(dlat1,1,thinning_grid(n,ifgat)%glat,thinning_grid(n,ifgat)%mlat,1)
    iy=int(dlat1)
    dy=dlat1-iy
    iy=max(1,min(iy,thinning_grid(n,ifgat)%mlat))

    call grdcrd(dlon1,1,thinning_grid(n,ifgat)%glon(1,iy),thinning_grid(n,ifgat)%mlon(iy),1)
    ix=int(dlon1)
    dx=dlon1-ix
    ix=max(1,min(ix,thinning_grid(n,ifgat)%mlon(iy)))

    dxx=half-min(dx,one-dx)
    dyy=half-min(dy,one-dy)
    dist1=dxx*dxx+dyy*dyy+half
    itx=thinning_grid(n,ifgat)%hll(ix,iy)
    itt=thinning_grid(n,ifgat)%istart_val(ithin)+itx
    if(ithin == 0) itt=0
    iuse=.true.
    if(dist1*crit1 > thinning_grid(n,ifgat)%score_crit(itx) .or. thinning_grid(n,ifgat)%icount(itx) == 0)iuse=.false.

    !write(6,'(a,3f10.3)') 'dlat_earth dlon_earth crit1 ',dlat_earth*rad2deg,dlon_earth*rad2deg,crit1
    !write(6,'(a,2i5,3f10.3,i10,e12.5,2x,L)') 'ix iy',ix,iy,dx,dy,dist1,itx,score_crit(itx),iuse
    return
  end subroutine map2tgrid

  subroutine grdcrd(d,nd,x,nx,flg)
  implicit none
  integer(i_kind),intent(in) :: nd,nx,flg
  real(r_kind),intent(inout) :: d
  real(r_kind),dimension(nx), intent(in) :: x

  integer(i_kind) :: id,ix
! Treat "normal" case in which nx>1
  if(nx>1) then
        if (flg.eq.1) then

!          Case in which x is in increasing order
           if(d<=x(1)) then
              ix=1
           else
              ix=isrchf(nx-1,x,d,flg)-1
           end if
           if(ix==nx) ix=ix-1

        else if (flg.eq.(-1)) then

!          Case in which x is in decreasing order
           if(d>=x(1)) then
              ix=1
           else
              ix=isrchf(nx-1,x,d,flg)-1
           end if
        end if
        d=float(ix)+(d-x(ix))/(x(ix+1)-x(ix))

! Treat special case of nx=1
  elseif (nx==1) then
        d = one
  endif

  return
end subroutine grdcrd

  subroutine checkob(n,ifgat,dist1,crit1,itx,iuse)
!
!   input argument list:
!     ifgat  - number of time slots for FGAT and 4DVAR
!     dist1  - quality indicator for distance (smaller = better)
!     crit1      - quality indicator for observation (smaller = better)
!     itx   - combined (i,j) index of observation on thinning grid
!     iuse  - .true. if observation should be used
!
!   output argument list:
!     iuse  - .true. if observation should be used
!
    implicit none
    logical,intent(inout):: iuse
    integer(i_kind),intent(in):: n,itx,ifgat
    real(r_kind),intent(in):: dist1,crit1

!   If data (no thinning), simply return to calling routine
    if(.not. iuse .or. thinning_grid(n,ifgat)%icount(itx)==0)return
    if(crit1*dist1 > thinning_grid(n,ifgat)%score_crit(itx))iuse=.false.

    return
  end subroutine checkob

  subroutine finalcheck(n,ifgat,dist1,crit1,iobs,itx,iobsout,iuse,sis)
!
!   input argument list:
!     ifgat  - number of time slots for FGAT and 4DVAR
!     dist1  - quality indicator for distance (smaller = better)
!     crit1  - quality indicator for observation (smaller = better)
!     itx    - combined (i,j) index of observation on thinning grid
!     iobs   - observation counter
!     iuse   - .true. if observation should be used
!     sis    - sensor/instrument/satellite
!
!   output argument list:
!     iobs   - observation counter
!     iobsout- location for observation to be put
!     iuse   - .true. if observation should be used
!
    implicit none
    logical,intent(inout):: iuse
    integer(i_kind),intent(inout):: iobs,iobsout
    integer(i_kind),intent(in):: n,itx,ifgat
    real(r_kind),intent(in):: dist1,crit1
    character(20),intent(in):: sis

    real(r_kind) crit

    if(.not. iuse)return

!   If using all data (no thinning), simply return to calling routine

    crit=crit1*dist1

    if(thinning_grid(n,ifgat)%icount(itx) == 0)then

!   Increment obs counter

      if(iobs < thinning_grid(n,ifgat)%itxmax)then
       iobs=iobs+1
       iobsout=iobs
       thinning_grid(n,ifgat)%score_crit(itx)= crit
       thinning_grid(n,ifgat)%ibest_obs(itx) = iobs
       thinning_grid(n,ifgat)%icount(itx)=thinning_grid(n,ifgat)%icount(itx)+1
      else
       iuse = .false.
       write(6,*)' ndata > maxobs when reading data for ',sis,thinning_grid(n,ifgat)%itxmax
      end if

    else if(crit < thinning_grid(n,ifgat)%score_crit(itx))then
      thinning_grid(n,ifgat)%score_crit(itx)= crit
      iobsout=thinning_grid(n,ifgat)%ibest_obs(itx)
      thinning_grid(n,ifgat)%icount(itx)=thinning_grid(n,ifgat)%icount(itx)+1
    else
      iuse = .false.
    end if

    return
  end subroutine finalcheck

  subroutine destroygrids(n,ifgat)
    implicit none
    integer(i_kind), intent(in) :: n,ifgat
    deallocate(thinning_grid(n,ifgat)%mlon,thinning_grid(n,ifgat)%glat, &
               thinning_grid(n,ifgat)%glon,thinning_grid(n,ifgat)%hll)
    deallocate(thinning_grid(n,ifgat)%icount)
    deallocate(thinning_grid(n,ifgat)%ibest_obs)
    deallocate(thinning_grid(n,ifgat)%score_crit)
    return
  end subroutine destroygrids

  subroutine destroygrids_conv(n, thin_3d)
    implicit none
    integer(i_kind), intent(in) :: n
    logical, optional, intent(in) :: thin_3d
    deallocate(thinning_grid_conv(n)%mlon,thinning_grid_conv(n)%glat,thinning_grid_conv(n)%glon)
    deallocate(thinning_grid_conv(n)%icount)
    deallocate(thinning_grid_conv(n)%ibest_obs)
    deallocate(thinning_grid_conv(n)%score_crit)
    if( present (thin_3d) ) then
       deallocate(thinning_grid_conv(n)%hll_3d)
    else
       deallocate(thinning_grid_conv(n)%hll)
    end if
    return
  end subroutine destroygrids_conv

  subroutine cleangrids_conv(n)
    implicit none
    integer(i_kind), intent(in) :: n
    thinning_grid_conv(n)%icount(:)     = 0
    thinning_grid_conv(n)%ibest_obs(:)  = 0
    thinning_grid_conv(n)%score_crit(:) = 9.99e6_r_kind
    return
  end subroutine cleangrids_conv

  subroutine destroy_sfc(n,ifgat)
    implicit none
    integer(i_kind), intent(in) :: n,ifgat
    deallocate(thinning_grid(n,ifgat)%sli,thinning_grid(n,ifgat)%sno,thinning_grid(n,ifgat)%isli)
    return
  end subroutine destroy_sfc

function isrchf(nx1,x,y,flg)
!                .      .    .                                       .
! subprogram:    isrchf
!   prgmmr: parrish          org: np22                date: 1990-10-11
!
! abstract: get grid coordinates from monotonically increasing or
!           decreasing points
!
! program history log:
!   2005-03-07  treadon - add doc block
!
!   input argument list:
!     nx1    - number of input points
!     x      - grid values
!     y      - target value
!     flg    - marks order of values in x
!              (1=increasing, -1=decreasing)
!
!   output argument list:
!     isrchf  - array index of input grid value near target value
!
  implicit none
  integer(i_kind):: isrchf
  integer(i_kind),intent(in):: nx1
  integer(i_kind),intent(in):: flg
  real(r_kind),intent(in):: y
  real(r_kind),dimension(nx1),intent(in):: x

  integer(i_kind) k

  if(flg.eq.1) then
    do k=1,nx1
      if(y<=x(k)) then
        isrchf=k

        go to 100
      end if
    end do
  else
    do k=1,nx1
      if(y>=x(k)) then
         isrchf=k
        go to 100
      end if
    end do
  end if

  isrchf=nx1+1
  if(nx1<=0) isrchf=0

100 continue
  return
end  function isrchf

end module gsi_thinning

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
program main_stats
    implicit none
    integer, parameter :: nr_data = 21
    type :: statistics
        real :: mean, stddev
        integer :: n
    end type statistics
    integer :: i
    type(statistics) :: stats
    real, dimension(nr_data) :: data

    data = [ (100.0 + 0.3*i, i=-size(data)/2, size(data)/2) ]
    stats = compute_stats(data)
    print *, stats%mean, stats%stddev, stats%n

contains

    function compute_stats(data) result(stats)
        implicit none
        real, dimension(:), intent(in) :: data
        type(statistics) :: stats
        real :: data_sum, data_sum2

        stats%n = size(data)
        data_sum = sum(data)
        data_sum2 = sum(data**2)
        associate(n => stats%n)
            stats%mean = data_sum/n
            stats%stddev = sqrt((data_sum2 - data_sum**2/n)/(n - 1))
        end associate
    end function compute_stats

end program main_stats

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
type intpointer
  integer, pointer :: p
end type intpointer

!...
  type(intpointer), dimension(100)  :: parray

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


module rotflip_rank
!! this demos assumed rank operations
!! generally this requires a compiler released in year 2020 or later including:
!!
!! * Gfortran >= 10
!! * Intel Fortran >= 20.0
implicit none (type, external)

contains

subroutine flipud(A)

integer, intent(inout) :: A(..)

select rank(A)
  rank (0)
    ! do nothing
  rank (1)
    A = A(ubound(A,1):lbound(A,1):-1)
  rank (2)
    A = A(ubound(A,1):lbound(A,1):-1, :)
  rank default
    error stop 'flipud: only rank 0..2 is handled for now.'
end select

end subroutine flipud

end module rotflip_rank

program assumed_rank

use rotflip_rank, only: flipud

implicit none (type, external)

integer :: i, A0 = 0, A1(3) = [0,1,2]
integer :: A2(3,3) = reshape( &
[0, 1, 2, &
 3, 4, 5, &
 6, 7, 8], &
 shape(A2), order=[2,1])
integer :: A3(3,3,3)

A3(:,:,1) = A2
A3(:,:,2) = A2+9
A3(:,:,3) = A2+18

call flipud(A0)

print '(A20,3I3)', '1D vector', A1
call flipud(A1)
print '(A20,3I3)', '1D vector reverse', A1

print *, '2D array'
do i = 1,size(A2,1)
  print '(3I3)', A2(i,:)
enddo
call flipud(A2)
print *, '2D array flipud'
do i = 1,size(A2,1)
  print '(3I3)', A2(i,:)
enddo

end program assumed_rank


text:
C     ------- ------
      PROGRAM AUTCON
C
C This program converts a data file into a labeled AUTO solution
C
C
      IMPLICIT NONE
      include 'fcon.h'
C
      INTEGER NDIM,NOLD,NTST,NCOL,ISW,IPS,NPAR,I,J,ios
      DOUBLE PRECISION TEMP,T,PERIOD
      DOUBLE PRECISION, ALLOCATABLE :: U(:),TM(:),UPS(:,:)
      DOUBLE PRECISION, ALLOCATABLE :: TMR(:),UPSR(:,:),PAR(:)
C
       OPEN(2,FILE='fort.2',STATUS='old',ACCESS='sequential',IOSTAT=ios)
       IF(ios/=0)THEN
          WRITE(6,'(A,A)')'The constants file (fort.2 or c. file) ',
     *         'could not be found.'
          STOP
       ENDIF
       OPEN(3,FILE='fort.3',STATUS='unknown',ACCESS='sequential')
       OPEN(8,FILE='fort.8',STATUS='unknown',ACCESS='sequential')
C
        CALL INIT(NDIM,IPS,NTST,NCOL,ISW,NPAR)
        ALLOCATE(U(NDIM),PAR(NPAR))
        PAR(:)=0d0
C
        NOLD=0
        DO
          NOLD=NOLD+1
          READ(3,*,END=2)TEMP,(TEMP,I=1,NDIM)
        ENDDO
C
 2      NOLD=NOLD-2
        ALLOCATE(TM(0:NTST),TMR(0:NOLD))
        ALLOCATE(UPSR(NDIM,0:NOLD),UPS(NDIM,0:NTST*NCOL))
C
        REWIND 3
        DO J=0,NOLD
          READ(3,*)TMR(J),UPSR(1:NDIM,J)
        ENDDO
C
        PERIOD=TMR(NOLD)-TMR(0)
        DO I=NOLD,0,-1
          TMR(I)=(TMR(I)-TMR(0))/PERIOD
        ENDDO
        CALL ADAPT(NOLD,NDIM,NTST,NCOL,NDIM,TMR,UPSR,TM,UPS,
     &       ((IPS.EQ.2.OR.IPS.EQ.12) .AND. ABS(ISW).LE.1))
C
        PAR(11)=PERIOD
        CALL STPNT(NDIM,U,PAR,T)
        CALL WRTBV8(NDIM,NTST,NCOL,ISW,PAR,UPS,TM,NPAR)
C
        DEALLOCATE(U,TM,TMR,UPS,UPSR)
C
      STOP
      CONTAINS

C     ---------- ----
      SUBROUTINE INIT(NDIM,IPS,NTST,NCOL,ISW,NPAR)
C
      IMPLICIT NONE
      INTEGER, INTENT(OUT) :: NDIM,IPS,NTST,NCOL,ISW,NPAR
      CHARACTER(LEN=2048) :: STR
C
C Reads the file of continuation constants
C
      INTEGER KEYEND,POS,NPOS
      INTEGER LINE
C
C Defaults
C
      NDIM=2
      IPS=1
      NTST=20
      NCOL=4
      ISW=1
      NPAR=NPARX

      LINE=0
      NPOS=1
      scanloop: DO
         IF(NPOS==1)THEN
            LINE=LINE+1
            READ(2,'(A)',END=1) STR
         ELSE
            STR=STR(NPOS:)
         ENDIF
         STR=ADJUSTL(STR)
         IF(LEN_TRIM(STR)==0)CYCLE
         DO I=1,LEN_TRIM(STR)
            ! comment on line
            IF(STR(I:I)=='#'.OR.STR(I:I)=='!')THEN
               NPOS=1
               CYCLE scanloop
            ENDIF
            ! keyword detected
            IF((LGE(STR(I:I),'A').AND.LLE(STR(I:I),'Z')).OR.
     &         (LGE(STR(I:I),'a').AND.LLE(STR(I:I),'z')))THEN
               STR=STR(I:)
               EXIT
            ELSE
               EXIT scanloop
            ENDIF
            IF(I==LEN_TRIM(STR))THEN
               NPOS=1
               CYCLE scanloop
            ENDIF
         ENDDO
         ! look for = after keyword
         KEYEND=SCAN(STR,'= ')-1
         IF(KEYEND==-1)THEN
            LINE=LINE-1
            EXIT scanloop
         ENDIF
         POS=SCAN(STR,'=')+1
         STR(POS:)=ADJUSTL(STR(POS:))
         NPOS=SCANVALUE(STR(POS:))
         IF(NPOS/=1)THEN
            NPOS=NPOS+POS-1
         ENDIF
         SELECT CASE(STR(1:KEYEND))
         CASE('NDIM')
            READ(STR(POS:),*,ERR=3)NDIM
         CASE('IPS')
            READ(STR(POS:),*,ERR=3)IPS
         CASE('NTST')
            READ(STR(POS:),*,ERR=3)NTST
         CASE('NCOL')
            READ(STR(POS:),*,ERR=3)NCOL
         CASE('ISW')
            READ(STR(POS:),*,ERR=3)ISW
         CASE('NPAR')
            READ(STR(POS:),*,ERR=3)NPAR
         ! ignore all others
         END SELECT
      ENDDO scanloop

      BACKSPACE 2
      READ(2,*,ERR=3,END=4) NDIM,IPS
      LINE=LINE+1
      READ(2,*,ERR=3,END=4)
      LINE=LINE+1
      READ(2,*,ERR=3,END=4) NTST,NCOL,ISW,ISW,ISW

 1    CONTINUE
      RETURN
 3    WRITE(6,"(A,I2,A)")
     *     " Error in fort.2 or c. file: bad value on line ",
     *     LINE,"."
      STOP
 4    WRITE(6,"(A,I2,A)")
     *     " Error in fort.2 or c. file: ends prematurely on line ",
     *     LINE,"."
      STOP
      END SUBROUTINE INIT

C     ------- -------- ---------
      INTEGER FUNCTION SCANVALUE(STR)
      IMPLICIT NONE
C
C     Scans STR(:) for a value
C     NPOS points to the next keyword on the same line,
C       or is set to 1 if there is none
C
      CHARACTER(*), INTENT(INOUT) :: STR

      INTEGER NPOS,I,LEVEL,LENSTR,ios
      CHARACTER(1) C,PREV,QUOTE
      LOGICAL QUOTEESC
      LEVEL=0
      QUOTE=' '
      QUOTEESC=.FALSE.
      PREV=' '

      NPOS=1
      LENSTR=LEN_TRIM(STR)
      I=1
      DO
         IF(I>LENSTR)THEN
            IF(LEVEL==0)EXIT
            LENSTR=LEN_TRIM(STR)
            READ(2,'(A)',IOSTAT=ios) STR(LENSTR+1:)
            IF(ios/=0)EXIT
            LENSTR=LEN_TRIM(STR)
         ENDIF
         NPOS=I
         C=STR(I:I)
         IF(QUOTE==' ')THEN
            SELECT CASE(C)
            CASE(',',' ')
               IF(LEVEL==0)EXIT
               IF(PREV==':')C=PREV !eat ',' and ' ' after ':'
            CASE(':')
            CASE(']','}')
               LEVEL=LEVEL-1
            CASE DEFAULT
               SELECT CASE(C)
               CASE('[','{')
                  LEVEL=LEVEL+1
               CASE('"',"'")
                  QUOTE=C
               END SELECT
            END SELECT
         ELSEIF(C==QUOTE)THEN
            ! ignore "" and ''
            IF(STR(I+1:I+1)==C.OR.QUOTEESC)THEN
               QUOTEESC=.NOT.QUOTEESC
            ELSE
               QUOTE=' '
            ENDIF
         ENDIF
         PREV=C
         I=I+1
      ENDDO
      NPOS=NPOS+VERIFY(STR(NPOS:)," ,")-1
      IF(NPOS>=LEN_TRIM(STR))NPOS=1
      SCANVALUE=NPOS
      END FUNCTION SCANVALUE
C
C     ---------- -----
      SUBROUTINE ADAPT(NOLD,NDOLD,NNEW,NCNEW,NDIM,
     &     TMR,UPSR,TM,UPS,IPER)

C Adapts the distribution of the mesh points so that the increase of the
C monotone function EQDF becomes approximately equidistributed over the
C intervals. The function UPS is interpolated on new mesh.

      INTEGER, INTENT(IN) :: NOLD,NDOLD,NNEW,NCNEW,NDIM
      LOGICAL, INTENT(IN) :: IPER
      DOUBLE PRECISION, INTENT(IN) :: UPSR(NDOLD,0:*), TMR(0:NOLD)
      DOUBLE PRECISION, INTENT(OUT) :: UPS(NDIM,0:*), TM(0:*)

C Generate the new mesh :

       CALL NEWMSH(NDIM,NDOLD,UPSR,NOLD,TMR,NNEW,TM,IPER)

C Replace UPS by its interpolant on the new mesh :

       CALL INTERP(NDIM,NOLD,TMR,UPSR,NDOLD,NNEW,NCNEW,TM,UPS)

      END SUBROUTINE ADAPT

C     ---------- ------
      SUBROUTINE INTERP(NDIM,N,TM,UPS,NDOLD,N1,NC1,TM1,UPS1)

C Finds interpolant (TM(.) , UPS(.) ) on new mesh TM1.

      INTEGER, INTENT(IN) :: NDIM,N,N1,NC1,NDOLD
      DOUBLE PRECISION, INTENT(IN) :: UPS(NDOLD,0:N)
      DOUBLE PRECISION, INTENT(IN) :: TM(0:N),TM1(0:N1)
      DOUBLE PRECISION, INTENT(OUT) :: UPS1(NDIM,0:N1*NC1)
C Local
      INTEGER I,J,J1,K
      DOUBLE PRECISION W0,W1,Z

       J=1
       DO J1=0,N1-1
         DO I=0,NC1-1
           Z=TM1(J1)+I*( TM1(J1+1)-TM1(J1) )/NC1
           DO
             IF(J>N)EXIT
             IF(TM(J)>Z)EXIT
             J=J+1
           ENDDO
           J=J-1
           W0=( Z-TM(J+1) )/( TM(J)-TM(J+1) )
           W1=( Z-TM(J) )/( TM(J+1)-TM(J) )
           DO K=1,NDIM
             UPS1(K,J1*NC1+I)=W0*UPS(K,J)+W1*UPS(K,J+1)
           ENDDO
         ENDDO
       ENDDO

       DO I=1,NDIM
         UPS1(I,N1*NC1)=UPS(I,N)
       ENDDO

      END SUBROUTINE INTERP

C     ---------- ------
      SUBROUTINE NEWMSH(NDIM,NDOLD,UPS,NOLD,TMOLD,NNEW,TMNEW,IPER)

C Redistributes the mesh according to the function EQDF.

      INTEGER, INTENT(IN) :: NDIM,NDOLD,NOLD,NNEW
      LOGICAL, INTENT(IN) :: IPER
      DOUBLE PRECISION, INTENT(IN) :: UPS(NDOLD,0:*),TMOLD(0:NOLD)
      DOUBLE PRECISION, INTENT(OUT) :: TMNEW(0:NNEW)
C Local
      INTEGER J,J1
      DOUBLE PRECISION X,DAL,UNEQ
      DOUBLE PRECISION, ALLOCATABLE :: EQF(:)
      ALLOCATE(EQF(0:NOLD))

C Put the values of the monotonely increasing function EQDF in EQF.

       CALL EQDF(NOLD,NDIM,TMOLD,UPS,NDOLD,EQF,IPER)

C Uniformly divide the range of EQDF :

       DAL=EQF(NOLD)/NNEW

C Generate the new mesh in TMNEW :

       J=1
       DO J1=0,NNEW-1

C EQF is an ascending array with values in [0,1]. Get the
C value of the index of the TM-interval in which EQF(i) lies.
        UNEQ=J1*DAL
        DO
          IF(J>NOLD)EXIT
          IF(EQF(J)>UNEQ)EXIT
          J=J+1
        ENDDO
        J=J-1
        X=(UNEQ-EQF(J))/(EQF(J+1)-EQF(J))
        TMNEW(J1)=(1.d0-X)*TMOLD(J)+X*TMOLD(J+1)
      ENDDO

C Assign TMNEW(NNEW) explicitly because of loss of precision
C problems when EQF(NOLD) and EQF(NOLD-1) are very close

      TMNEW(NNEW)=TMOLD(NOLD)

      DEALLOCATE(EQF)
      END SUBROUTINE NEWMSH

C     ---------- ----
      SUBROUTINE EQDF(NTST,NDIM,TM,UPS,NDOLD,EQF,IPER)

      DOUBLE PRECISION, PARAMETER :: HMACH=1.0d-7

      INTEGER, INTENT(IN) :: NTST,NDIM,NDOLD
      LOGICAL, INTENT(IN) :: IPER
      DOUBLE PRECISION, INTENT(IN) :: UPS(NDOLD,0:NTST),TM(0:NTST)
      DOUBLE PRECISION, INTENT(OUT) :: EQF(0:NTST)

C Local
      DOUBLE PRECISION E,DTAV,ND,DT
      DOUBLE PRECISION, ALLOCATABLE :: HD(:)
      LOGICAL SMALL
      INTEGER I,J

       ALLOCATE(HD(NDIM))

C Compute approximation to NCOL-th and (NCOL+1)-st derivative
C and define the equidistribution function :

       SMALL=.TRUE.
       EQF(0)=0.d0
       DO J=1,NTST
         E=0.d0
         DT=TM(J)-TM(J-1)
         DO I=1,NDIM
           ND=(UPS(I,J)-UPS(I,J-1))/DT
           IF(J>1)THEN
             E=E+SQRT(ABS( 2*( ND-HD(I) )/(TM(J)-TM(J-2)) ))
           ENDIF
           IF(ABS(ND)>HMACH)SMALL=.FALSE.
           HD(I)=ND
         ENDDO
         IF(J>1)EQF(J-1)=EQF(J-2)+(TM(J-1)-TM(J-2))*E
       ENDDO

       E=0.d0
       J=NTST-1
       DTAV=TM(NTST)-TM(NTST-1)
       IF(IPER)THEN
C      *Extend by periodicity :
         J=1
         DTAV=(DTAV+TM(1)-TM(0))/2
       ENDIF
C      *else extend by extrapolation :
       DO I=1,NDIM
         ND=(UPS(I,J)-UPS(I,J-1))/(TM(J)-TM(J-1))
         E=E+SQRT(ABS( (ND-HD(I))/DTAV ))
         IF(ABS(ND)>HMACH)SMALL=.FALSE.
       ENDDO
       EQF(NTST)=EQF(NTST-1)+(TM(NTST)-TM(NTST-1))*E

C Take care of "small derivative" case.

       IF(SMALL)THEN
         DO I=0,NTST
            EQF(I)=I
         ENDDO
       ENDIF

       DEALLOCATE(HD)
      END SUBROUTINE EQDF

C     ---------- ------
      SUBROUTINE WRTBV8(NDIM,NTST,NCOL,ISW,PAR,UPS,TM,NPAR)

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: NDIM,NTST,NCOL,ISW,NPAR
      DOUBLE PRECISION, INTENT(IN) :: UPS(NDIM,0:*),TM(0:*),PAR(*)

      INTEGER IBR,ITP,NFPR,NTOT,LAB,NTPL,NAR,NROWPR,MTOT

       ITP=9
       NFPR=1
       IBR=1
       NTOT=1
       LAB=1

C Write information identifying the solution :

       NTPL=NCOL*NTST+1
       NAR=NDIM+1
       NROWPR=((NDIM+6)/7)*NTPL + (NPAR+6)/7
       MTOT=MOD(NTOT-1,9999)+1
       WRITE(8,101)IBR,MTOT,ITP,LAB,NFPR,ISW,NTPL,NAR,NROWPR,
     &      NTST,NCOL,NPAR

C Write the entire solution on unit 8 :

       DO J=0,NTST*NCOL
         T=TM(J/NCOL)+MOD(J,NCOL)*(TM(J/NCOL+1)-TM(J/NCOL))/NCOL
         WRITE(8,102)T,UPS(:,J)
       ENDDO

C Write the parameter values.

       WRITE(8,102)(PAR(I),I=1,NPAR)

 101   FORMAT(6I6,I8,I6,I8,3I5)
 102   FORMAT(4X,7ES19.10)

      END SUBROUTINE WRTBV8

      END PROGRAM AUTCON

C
      DOUBLE PRECISION FUNCTION GETP(CODE,IC,UPS)
C Dummy routine (in case GETP is called in the user routine PVLS)
      GETP=0.0d0
      END FUNCTION GETP


! test_bobyqa.f90 --
!
!     Test problem for BOBYQA, the objective function being the sum of
!     the reciprocals of all pairwise distances between the points P_I,
!     I=1,2,...,M in two dimensions, where M=N/2 and where the components
!     of P_I are X(2*I-1) and X(2*I). Thus each vector X of N variables
!     defines the M points P_I. The initial X gives equally spaced points
!     on a circle. Four different choices of the pairs (N,NPT) are tried,
!     namely (10,16), (10,21), (20,26) and (20,41). Convergence to a local
!     minimum that is not global occurs in both the N=10 cases. The details
!     of the results are highly sensitive to computer rounding errors. The
!     choice IPRINT=2 provides the current X and optimal F so far whenever
!     RHO is reduced. The bound constraints of the problem require every
!     component of X to be in the interval [-1,1].
!
module bobyqa_functions
    use bobyqa_minimisation, only: wp
    implicit none

contains
subroutine calfun (x,f)
    real(kind=wp), dimension(:) :: x
    real(kind=wp)               :: f

    integer                     :: i, j, n
    real(kind=wp)               :: temp

    n = size(x)
    f = 0.0d0
    do i = 4,n,2
        do j=2,i-2,2
            temp = (x(i-1)-x(j-1))**2+(x(i)-x(j))**2
            temp = max( temp, 1.0d-6 )
            f    = f + 1.0d0 / sqrt(temp)
        enddo
    enddo
end subroutine calfun
end module bobyqa_functions

program test_bobyqa
    use bobyqa_minimisation
    use bobyqa_functions

    implicit none
    real(kind=wp), dimension(100) :: x, xl, xu
    real(kind=wp)                 :: twopi, bdl, bdu, rhobeg, rhoend, temp
    integer                       :: iprint, maxfun
    integer                       :: m, n, j, jcase, npt

    twopi  =  8.0d0 * atan(1.0d0)
    bdl    = -1.0d0
    bdu    =  1.0d0
    iprint =  2
    maxfun =  500000
    rhobeg =  1.0d-1
    rhoend =  1.0d-6

    m = 5
    do while ( m <= 10 )
        n  = 2 * m
        xl = bdl
        xu = bdu

        do jcase = 1,2
            npt = n + 6
            if ( jcase == 2 ) npt = 2 * n + 1
            print 30, m,n,npt
   30       format (//5x,'2d output with m =',i4,',  n =',i4, &
                         '  and  npt =',i4)
            do j = 1,m
                temp = real(j,wp) * twopi / real(m,wp)
                x(2*j-1) = cos(temp)
                x(2*j)   = sin(temp)
            enddo

            call bobyqa( npt, x(1:n), xl(1:n), xu(1:n), rhobeg, rhoend, iprint, maxfun, calfun)
        enddo
        m = 2 * m
    enddo
end program


text:
      subroutine lint(x,y,xx,yy,n1,n2,kt,ikop)
C...subroutine lint 01/09/78
      common/gfcell/i1,j1,i2,j2,d1,d2,kent,qq,const
       common/gfint/ is
       dimension x(n1),y(n2),xx(kt),yy(kt)
      go to(531,532,531,532),kent
  531 rd=x(i1)
      tt=y(j1)+(y(j2)-y(j1))/(1.-d2/d1)
       go to 203
532   rd=x(i1)+(x(i2)-x(i1))/(1.-d2/d1)
      tt=y(j1)
  203 xx(is)=rd
      yy(is)=tt
      if(ikop.ne.1) goto 2
       xx(is)=xpol(rd,tt)
      yy(is)=ypol(rd,tt)
    2 return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
test03/test_sem.cpp
test03/main.cpp


text:
AmericanWest Bank is a banks bank that is headquartered in Spokane, Washington.

This branch location was previously a Premier West Bank. Before that it was a Wachovia branch until Wells Fargo sold the branch to PremierWest Bank. Originally this location was a World Savings Bank branch until October 11, 2007.

AmericanWest Bank took over PremierWest Bank and rebranded the bank branches in 2013.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search