text:
      function RP1(m_N, tau, taup, term, j_chi, coupling_Array)
      implicit none
      double precision :: RP1, m_N, j_chi
      double precision :: coupling_Array(14,2)
      integer :: tau,taup
      integer :: term
      integer :: c, q2
      double precision :: c12, c12p
      double precision :: c13, c13p
c
      c = 0
      q2 = 2
c
      c12 = coupling_Array(11,tau)
      c12p = coupling_Array(11,taup)
      c13 = coupling_Array(12,tau)
      c13p = coupling_Array(12,taup)
c
      if (term.eq.c) then
         RP1 = (j_chi*(j_chi+1))/12. * c12*c12p
      else if (term.eq.q2) then
         RP1 = 1/m_N**2 * (j_chi*(j_chi+1))/12. * c13*c13p
      else
         RP1 = 0.
      end if
c
      end function RP1

ext: f


text:
! *****************************COPYRIGHT*******************************
! (C) Crown copyright Met Office. All rights reserved.
! For further details please refer to the file COPYRIGHT.txt
! which you should have received as part of this distribution.
! *****************************COPYRIGHT*******************************
!
! Driver for 1D Kinematic Driver (KiD) model
!
! Author: Ben Shipway
!
! For version details see header_data.f90

Program KiD

  Use main

  call main_loop

end Program KiD

max_stars_repo_name: Adehill/KiD-A

max_issues_repo_name: Adehill/KiD-A

max_forks_repo_name: Adehill/KiD-A


text:
C
C $Id: trn32s.f,v 1.6 2008-07-27 00:17:26 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE TRN32S (X,Y,Z,XT,YT,ZT,IFLAG)
C
C THIS ROUTINE IMPLEMENTS THE 3-SPACE TO 2-SPACE TRANSFOR-
C MATION BY KUBER, SZABO AND GIULIERI, THE PERSPECTIVE
C REPRESENTATION OF FUNCTIONS OF TWO VARIABLES. J. ACM 15,
C 2, 193-204,1968.
C IFLAG=0 ARGUMENTS
C X,Y,Z    ARE THE 3-SPACE COORDINATES OF THE INTERSECTION
C          OF THE LINE OF SIGHT AND THE IMAGE PLANE.  THIS
C          POINT CAN BE THOUGHT OF AS THE POINT LOOKED AT.
C XT,YT,ZT ARE THE 3-SPACE COORDINATES OF THE EYE POSITION.
C
C IFLAG=1 ARGUMENTS
C X,Y,Z    ARE THE 3-SPACE COORDINATES OF A POINT TO BE
C          TRANSFORMED.
C XT,YT    THE RESULTS OF THE 3-SPACE TO 2-SPACE TRANSFOR-
C          MATION.
C          USE INT(XT) AND INT(YT) IN GPL CALLS.
C ZT       NOT USED.
C IF LL (IN COMMON) =0 XT AND YT ARE IN THE SAME SCALE AS X, Y, AND Z.
C
      COMMON /PWRZ1S/ XXMIN      ,XXMAX      ,YYMIN      ,YYMAX      ,
     1                ZZMIN      ,ZZMAX      ,DELCRT     ,EYEX       ,
     2                EYEY       ,EYEZ
      COMMON /SRFBLK/ LIMU(1024) ,LIML(1024) ,CL(41)     ,NCL        ,
     1                LL         ,FACT       ,IROT       ,NDRZ       ,
     2                NUPPER     ,NRSWT      ,BIGD       ,UMIN       ,
     3                UMAX       ,VMIN       ,VMAX       ,RZERO      ,
     4                IOFFP      ,NSPVAL     ,SPVAL      ,BIGEST
      DIMENSION       NLU(7)     ,NRU(7)     ,NBV(7)     ,NTV(7)
C
C SAVE INSERTED BY BEN DOMENICO 9/8/85 BECAUSE OF ASSUMPTION THAT
C   JUMP, JUMP2, AND JUMP3 ARE PRESERVED BETWEEN CALLS.
C   THERE MAY BE OTHER SUCH ASSUMPTIONS AS WELL.
C
      SAVE
C
C PICTURE CORNER COORDINATES FOR LL=1
C
      DATA NLU(1),NRU(1),NBV(1),NTV(1)/  10,1014,  10,1014/
C
C PICTURE CORNER COORDINATES FOR LL=2
C
      DATA NLU(2),NRU(2),NBV(2),NTV(2)/  10, 924,  50, 964/
C
C PICTURE CORNER COORDINATES FOR LL=3
C
      DATA NLU(3),NRU(3),NBV(3),NTV(3)/ 100,1014,  50, 964/
C
C PICTURE CORNER COORDINATES FOR LL=4
C
      DATA NLU(4),NRU(4),NBV(4),NTV(4)/  10,1014,  10,1014/
C
C PICTURE CORNER COORDINATES FOR LL=5
C
      DATA NLU(5),NRU(5),NBV(5),NTV(5)/  10,1014,  10,1014/
C
C PICTURE CORNER COORDINATES FOR LL=6
C
      DATA NLU(6),NRU(6),NBV(6),NTV(6)/  10, 512, 256, 758/
C
C PICTURE CORNER COORDINATES FOR LL=7
C
      DATA NLU(7),NRU(7),NBV(7),NTV(7)/ 512,1014, 256, 758/
C
C STORE THE PARAMETERS OF THE SET32 CALL FOR USE WHEN
C TRN32 IS CALLED.
C
      IF (IFLAG)  40, 10, 40
   10 CONTINUE
      JUMP3 = 1
      IF (IOFFP .EQ. 1) JUMP3 = 2
      AX = X
      AY = Y
      AZ = Z
      EX = XT
      EY = YT
      EZ = ZT
C
C AS MUCH COMPUTATION AS POSSIBLE IS DONE DURING EXECUTION
C THIS ROUTINE WHEN IFLAG=0 BECAUSE CALLS IN THAT MODE ARE INFREQUENT.
C
      DX = AX-EX
      DY = AY-EY
      DZ = AZ-EZ
      D = SQRT(DX*DX+DY*DY+DZ*DZ)
      COSAL = DX/D
      COSBE = DY/D
      COSGA = DZ/D
      SINGA = SQRT(1.-COSGA*COSGA)
      JUMP2 = 1
      IF (LL .EQ. 0) GO TO  20
      JUMP2 = 2
      DELCRT = NRU(LL)-NLU(LL)
      U0 = UMIN
      V0 = VMIN
      U1 = NLU(LL)
      V1 = NBV(LL)
      U2 = NRU(LL)-NLU(LL)
      V2 = NTV(LL)-NBV(LL)
      U3 = U2/(UMAX-UMIN)
      V3 = V2/(VMAX-VMIN)
      U4 = NRU(LL)
      V4 = NTV(LL)
      IF (NRSWT .EQ. 0) GO TO  20
      U0 = -BIGD
      V0 = -BIGD
      U3 = U2/(2.*BIGD)
      V3 = V2/(2.*BIGD)
C
C THE 3-SPACE POINT LOOKED AT IS TRANSFORMED INTO (0,0) OF
C THE 2-SPACE.  THE 3-SPACE Z AXIS IS TRANSFORMED INTO THE
C 2-SPACE Y AXIS.  IF THE LINE OF SIGHT IS CLOSE TO PARALLEL
C TO THE 3-SPACE Z AXIS, THE 3-SPACE Y AXIS IS CHOSEN (IN-
C STEAD OF THE 3-SPACE Z AXIS) TO BE TRANSFORMED INTO THE
C 2-SPACE Y AXIS.
C
   20 IF (SINGA .LT. 0.0001) GO TO  30
      R = 1./SINGA
      JUMP = 1
      RETURN
   30 SINBE = SQRT(1.-COSBE*COSBE)
      R = 1./SINBE
      JUMP = 2
      RETURN
   40 CONTINUE
      XX = X
      YY = Y
      ZZ = Z
      GO TO ( 60, 50) , JUMP3
   50 IF (ZZ .EQ. SPVAL) GO TO 110
   60 Q = D/((XX-EX)*COSAL+(YY-EY)*COSBE+(ZZ-EZ)*COSGA)
      GO TO (70,80), JUMP
   70 XX = ((EX+Q*(XX-EX)-AX)*COSBE-(EY+Q*(YY-EY)-AY)*COSAL)*R
      YY = (EZ+Q*(ZZ-EZ)-AZ)*R
      GO TO 90
   80 XX = ((EZ+Q*(ZZ-EZ)-AZ)*COSAL-(EX+Q*(XX-EX)-AX)*COSGA)*R
      YY = (EY+Q*(YY-EY)-AY)*R
   90 GO TO (120,100), JUMP2
  100 XX = MIN(U4,MAX(U1,U1+U3*(FACT*XX-U0)))
      YY = MIN(V4,MAX(V1,V1+V3*(FACT*YY-V0)))
      GO TO 120
  110 XX = NSPVAL
      YY = NSPVAL
C
  120 XT = XX
      YT = YY
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
! Check the Flang Print Function Names example plugin doesn't count/print function/subroutine calls (should only count definitions)
! This requires that the examples are built (FLANG_BUILD_EXAMPLES=ON) to access flangPrintFunctionNames.so

! REQUIRES: plugins, examples, shell

! RUN: %flang_fc1 -load %llvmshlibdir/flangPrintFunctionNames%pluginext -plugin print-fns %s 2>&1 | FileCheck %s

!-----------------------------
! EXPECTED OUTPUT: Counts == 0
!-----------------------------
! CHECK: ==== Functions: 0 ====
! CHECK-NEXT: ==== Subroutines: 0 ====

!-----------------------------
! INPUT
!-----------------------------
program main
    call subroutine1
    fn1 = function1()
    fn2 = function2()
end program main

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm


text:
*DECK C0LGMC
      COMPLEX FUNCTION C0LGMC (Z)
C***BEGIN PROLOGUE  C0LGMC
C***PURPOSE  Evaluate (Z+0.5)*LOG((Z+1.)/Z) - 1.0 with relative
C            accuracy.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7A
C***TYPE      COMPLEX (C0LGMC-C)
C***KEYWORDS  FNLIB, GAMMA FUNCTION, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Evaluate  (Z+0.5)*LOG((Z+1.0)/Z) - 1.0  with relative error accuracy
C Let Q = 1.0/Z so that
C     (Z+0.5)*LOG(1+1/Z) - 1 = (Z+0.5)*(LOG(1+Q) - Q + Q*Q/2) - Q*Q/4
C        = (Z+0.5)*Q**3*C9LN2R(Q) - Q**2/4,
C where  C9LN2R  is (LOG(1+Q) - Q + 0.5*Q**2) / Q**3.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  C9LN2R, R1MACH
C***REVISION HISTORY  (YYMMDD)
C   780401  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C***END PROLOGUE  C0LGMC
      COMPLEX Z, Q, C9LN2R
      SAVE RBIG
      DATA RBIG / 0.0 /
C***FIRST EXECUTABLE STATEMENT  C0LGMC
      IF (RBIG.EQ.0.0) RBIG = 1.0/R1MACH(3)
C
      CABSZ = ABS(Z)
      IF (CABSZ.GT.RBIG) C0LGMC = -(Z+0.5)*LOG(Z) - Z
      IF (CABSZ.GT.RBIG) RETURN
C
      Q = 1.0/Z
      IF (CABSZ.LE.1.23) C0LGMC = (Z+0.5)*LOG(1.0+Q) - 1.0
      IF (CABSZ.GT.1.23) C0LGMC = ((1.+.5*Q)*C9LN2R(Q) - .25) * Q**2
C
      RETURN
      END


! -*- f90 -*-
! file lbfgs.f90
! Approximate Hessian inversion
! version 0.1
SUBROUTINE lbfgs(n, H_0, k, s, y, d)

!*****************************************************************
!
!   subroutine lbfgs
!
!   Approximately apply the inverse Hessian to d.
!   n: dim of H_0 (nxn) and d (n)
!   d: Vector to mupltiply by. On exit, the result of the multiplication
!
!*****************************************************************
    IMPLICIT NONE

    INTEGER ::  k, n
    REAL (KIND=8) ::  H_0(n,n), s(n,k), y(n,k), d(n)

    INTEGER :: i,j
    REAL (KIND=8) ::  a(k), b(k), rho(k)

!    print *, n,k !    print *,'s',s !    print *,'y',y

    do i=1,k-1
        rho(i) = 1.0/DOT_PRODUCT(y(:,i), s(:,i))
    end do

    do i=k-1, 1,-1
        a(i) =  rho(i)*(DOT_PRODUCT(s(:,i), d))
        d = d - a(i)*y(:,i) ! might need to modify here?
    end do

    d = MATMUL(H_0, d)

    do i=1,k-1
       b(i) = rho(i)*(DOT_PRODUCT(y(:,i), d))
       !print *,i, b(i), rho(i)
       d = d + (a(i)-b(i))*s(:,i) ! might also need a tweak here
       !print *,'d',d
    end do

END SUBROUTINE lbfgs

SUBROUTINE update_storage(n, k, istep, s, y, d_x, d_jac, info)

    IMPLICIT NONE

    INTEGER :: n,k,istep

    REAL(KIND=8) :: s(n,k), y(n,k)
    REAL(KIND=8) :: d_x(n), d_jac(n)

    INTEGER :: i, info

    info = 1
    do i=1,n
        if (ABS(d_x(i)) .LE. 1.0E-015) then
            info = 0
        end if
    end do

    !print *,'s',s(1,:)
    !print *,'y',y(1,:)
    if (info .EQ. 1) then
        if (istep .LE. k) then
            !no shuffling needed
            s(:,istep) = d_x
            y(:,istep) = d_jac
        else
            do i=2,k
                s(:,i-1) = s(:,i)
                y(:,i-1) = y(:,i)
            end do
        s(:,k) = d_x
        y(:,k) = d_jac
        end if
    end if

    !print *,'s',s(1,:)
    !print *,'y',y(1,:)
!    print *,' '

END SUBROUTINE update_storage


text:
C  *********************************************************
C  *                                                       *
C  *    TEST NUMBER: 04.01.04/02                           *
C  *    TEST TITLE : Appearance of annotation text         *
C  *                 primitives                            *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /GLOBCH/ PIDENT,    GLBERR,    TSTMSG,     FUNCID,
     1                DUMCH
      CHARACTER       PIDENT*40, GLBERR*60, TSTMSG*900, FUNCID*80,
     1                DUMCH(20)*20

      COMMON /DIALOG/ DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT,
     2                DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS
      INTEGER         DOUTYP, DINTYP, DSTDNR, DSTRID, PSTRID, DTCLIM,
     1                SCRMOD, DTXCI, SPECWT
      REAL            DSIZE, EFRAC, DYXRAT, SYXRAT, MTRPDC, WCPDC, QVIS

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C text alignment horizontal
      INTEGER    PAHNOR,     PALEFT,     PACENT,     PARITE
      PARAMETER (PAHNOR = 0, PALEFT = 1, PACENT = 2, PARITE = 3)
C text alignment vertical
      INTEGER    PAVNOR,     PATOP,      PACAP,     PAHALF,
     1           PABASE,     PABOTT
      PARAMETER (PAVNOR = 0, PATOP  = 1, PACAP = 2, PAHALF = 3,
     1           PABASE = 4, PABOTT = 5)

C text path
      INTEGER    PRIGHT,     PLEFT,     PUP,     PDOWN
      PARAMETER (PRIGHT = 0, PLEFT = 1, PUP = 2, PDOWN = 3)

C text precision
      INTEGER    PSTRP,     PCHARP,     PSTRKP
      PARAMETER (PSTRP = 0, PCHARP = 1, PSTRKP = 2)

C line type
      INTEGER    PLSOLI,    PLDASH,     PLDOT,   PLDASD
      PARAMETER (PLSOLI=1,  PLDASH=2,   PLDOT=3, PLDASD=4)

C composition type
C                preconcatenate  postconcatenate  replace
      INTEGER    PCPRE,          PCPOST,          PCREPL
      PARAMETER (PCPRE = 0,      PCPOST = 1,      PCREPL = 2)

      INTEGER    PICSTR, TXCI, IX, NGSQ, PERM(10), RNDINT

      REAL       XMID,YMID,EXT,XEXT,YEXT,X0,Y0,Z0, XF(4,4), CHHT
      REAL       NPCX, NPCY, NPCPWC, TX(2), TY(2), YBOT,YTOP, NGDELT
      REAL       XREFPT,YREFPT,ZREFPT,XANNPT,YANNPT,ZANNPT, YB,YH,YC
      REAL       XSIZ,YSIZ,XWINLO(6),YWINLO(6), RDUM1, RDUM2
      REAL       RNDRL,SIDER(6,3),IDM(4,4),XFA(4,4),XFB(4,4)
      REAL       Z, H, U, PI
      PARAMETER  (Z=0., H=0.5, U=1.0, PI=3.14159265)

      CHARACTER  UPPER*26, CHSTR*5, LOWER*26, REPSTR*10, CHRABC*3

      DATA       SIDER/18*Z/

      CALL INITGL ('04.01.04/02')

C open PHIGS
      CALL XPOPPH (ERRFIL, MEMUN)
C set-up of workstation and dialogue area
      PICSTR = 101
      TXCI = 1
      CALL SETDLG (PICSTR, 801,TXCI)

      CALL POPST (PICSTR)
C  by convention, view #1 is for picture
      CALL PSVWI (1)
C  use individual attributes
      CALL SETASF (PINDIV)
C  ratio of NPC to WC
      CALL WCNPC (0., 0., NPCX, NPCY, NPCPWC)

C  set default attributes:
      CALL PSTXFN (1)
      CALL PSTXPR (PSTRKP)
      CALL PSATAL (PACENT, PAHALF)
      CALL PSATP  (PRIGHT)
      CHHT = 0.05
      CALL PSATCH (CHHT * NPCPWC)
      CALL PSANS  (1)
      CALL PSLN   (PLSOLI)
      CALL PSPLCI (3)

C  set up 102 as sub-structure
      CALL PEXST (102)
      CALL PCLST
C
C  *** *** *** *** ***  character content  *** *** *** *** ***
C
      CALL SETMSG ('7', 'A displayed annotation text primitive ' //
     1             'should contain the specified character string.')

      CALL POPST (102)

      CALL PSATCH (0.08 * NPCPWC)
      CALL PSTXFN (2)
      UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'

C  chstr = string of 5 randomly chosen upper-case letters.
      CALL RNSET (5, 26, PERM)
      DO 100 IX = 1, 5
         CHSTR(IX:IX) = UPPER(PERM(IX):PERM(IX))
100   CONTINUE

C  display <annotation text relative 3> consisting of chstr:
C    reference position = 0.5, 0.5, 0.5
C    offset = 0., 0., 0.1
      CALL PATR3 (H,H,H, Z,Z,Z, CHSTR)

      CALL DLINE ('CHARACTER CONTENT: Enter the displayed character ' //
     1            'string.', REPSTR)
      IF (REPSTR .NE. CHSTR) THEN
         CALL FAIL
         GOTO 199
      ENDIF

      CALL PEMST (102)
      CALL PSTXFN (1)
      CALL PSATCH (0.08 * NPCPWC)

C  chstr = string of 5 randomly chosen lower-case letters
      LOWER = 'abcdefghijklmnopqrstuvwxyz'
      CALL RNSET (5, 26, PERM)
      DO 110 IX = 1, 5
         CHSTR(IX:IX) = LOWER(PERM(IX):PERM(IX))
110   CONTINUE

C  display <annotation text relative> consisting of chstr:
C     reference point = 0.5, 0.5
C     offset = 0., 0.
      CALL PATR (H,H, Z,Z, CHSTR)

      CALL DLINE ('CHARACTER CONTENT: Enter the displayed character ' //
     1            'string.', REPSTR)
      CALL IFPF (REPSTR .EQ. CHSTR)

C  end_char:
199   CONTINUE
      CALL PEMST (102)

C  *** *** *** ***   effect of reference point   *** *** *** ***

      CALL PEXST (106)
      CALL PEXST (103)
      CALL PCLST

C  set up rotations for viewing from 6 sides
C  back - Y-axis
      SIDER (2,2) = PI
C  left - Y-axis
      SIDER (3,2) = 0.5*PI
C  right - Y-axis
      SIDER (4,2) = 1.5*PI
C  top - X-axis
      SIDER (5,1) = 0.5*PI
C  bottom - X-axis
      SIDER (6,1) = 1.5*PI

C  Divide screen up into 6 squares labelled areas
      CALL WIN6 (106, 2, XSIZ, YSIZ, XWINLO, YWINLO)

      CALL POPST (103)

      CHRABC = 'ABC'
      CALL IDMAT (4, IDM)
C  amount by which "incorrect" rectangle is to be shifted
      NGDELT = 0.01

C  <Inquire text extent> to get extent rectangle for CHRABC
      CALL PQTXX (SPECWT, 1, 1., 0., CHHT, PRIGHT,PACENT,PABASE,
     1            CHRABC, ERRIND, TX, TY, RDUM1, RDUM2)
      CALL CHKINQ ('pqtxx', ERRIND)
      YB = TY(1)
      CALL PQTXX (SPECWT, 1, 1., 0., CHHT, PRIGHT,PACENT,PAHALF,
     1            CHRABC, ERRIND, TX, TY, RDUM1, RDUM2)
      CALL CHKINQ ('pqtxx', ERRIND)
      YH = TY(1)
      CALL PQTXX (SPECWT, 1, 1., 0., CHHT, PRIGHT,PACENT,PACAP,
     1            CHRABC, ERRIND, TX, TY, RDUM1, RDUM2)
      CALL CHKINQ ('pqtxx', ERRIND)
      YC = TY(1)
C  horizontal size of string in MC
      XEXT = TX(2)-TX(1)
C  MC distance from baseline to halfline
      YBOT = ABS(YH-YB)
C  MC distance from capline to halfline
      YTOP = ABS(YC-YH)
C  MC vertical distance from baseline to capline
      YEXT = YTOP+YBOT

C  select reference point randomly from cube:
C    (ext,1-ext),(ext,1-ext),(ext,1-ext)
C  to make sure that annotation text fits in window

C  Distance from center of rectangle to farthest border, divided
C  by XSIZ to cancel effect of modelling transformation, which
C  multiplies by XSIZ.  Thus, when string is located, aligned at
C  its center, it won't overlap border.
      EXT = 0.5 * MAX(XEXT, YEXT) / XSIZ

C  which square is no good
      NGSQ = RNDINT (1,6)

C  reference point must be at least EXT distant from any side of the
C  unit cube in all dimensions, because modelling transformation will
C  rotate cube, and reference point along with it, but NOT annotation
C  text itself, which is always in NPC-XY plane.  Therefore, it is
C  incorrect (even though plausible) to apply XEXT in X-dimension,
C  and YEXT in Y-dimension, because rotation mixes up dimensions.
      XREFPT = RNDRL (EXT, 1-EXT)
      YREFPT = RNDRL (EXT, 1-EXT)
      ZREFPT = RNDRL (EXT, 1-EXT)

      DO 250 IX = 1, 6
         XMID = XWINLO(IX)+XSIZ/2
         YMID = YWINLO(IX)+YSIZ/2

C  scale down to fit in square, rotate, and move to square #ix
         CALL EBLTM3 (H, H, H, XMID-H, YMID-H, Z,
     2        SIDER(IX,1),SIDER(IX,2),SIDER(IX,3), XSIZ,XSIZ,XSIZ,XF)
         CALL PSLMT3 (XF, PCREPL)
C  draw the annotation text at transformed reference point
         CALL PATR3 (XREFPT,YREFPT,ZREFPT, Z,Z,Z, CHRABC)

         CALL PSLMT3 (IDM, PCREPL)
C  find expected center of rectangle
         CALL ETP3 (XREFPT,YREFPT,ZREFPT, XF, X0,Y0,Z0)

         IF (IX .EQ. NGSQ) THEN
C  distort towards middle of square - don't overlap border
            IF (X0 .LT. XMID) THEN
                X0 = X0 + NGDELT
            ELSE
                X0 = X0 - NGDELT
            ENDIF

            IF (Y0 .LT. YMID) THEN
                Y0 = Y0 + NGDELT
            ELSE
                Y0 = Y0 - NGDELT
            ENDIF
         ENDIF

C  outline the annotation text
         CALL DRWREC (X0-XEXT/2, X0+XEXT/2, Y0-YBOT, Y0+YTOP)
250   CONTINUE

      CALL SETMSG ('7 8 9', 'The location of an annotation text ' //
     1             'primitive should reflect its 3D reference '   //
     2             'point after being transformed.')
      CALL DCHPFV ('EFFECT OF REFERENCE POINT: In which box is '  //
     1             'the annotation text NOT accurately outlined ' //
     2             'by the rectangle?', 6, NGSQ)
      CALL PEMST(103)

C  *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** *** ***

C  which square will be marked wrong?
      NGSQ = RNDINT (1,6)

C  reference point:
      XREFPT = RNDRL (EXT, 1-EXT)
      YREFPT = RNDRL (EXT, 1-EXT)
      ZREFPT = Z

      DO 260 IX = 1,6
         XMID = XWINLO(IX)+XSIZ/2
         YMID = YWINLO(IX)+YSIZ/2
C  first, move from Z=0 plane to somewhere in middle
         CALL ETR3 (Z, Z, RNDRL(EXT,1-EXT), XFA)
C  scale down to fit in square, rotate, and move to square #ix
         CALL EBLTM3 (H,H,H, XMID-H, YMID-H, Z,
     2         SIDER(IX,1),SIDER(IX,2),SIDER(IX,3), XSIZ,XSIZ,XSIZ,XFB)
         CALL ECOM3 (XFB, XFA, XF)

         CALL PSLMT3 (XF, PCREPL)
         CALL PATR (XREFPT,YREFPT, Z,Z, CHRABC)

         CALL PSLMT3 (IDM, PCREPL)
C  find expected center of rectangle
         CALL ETP3 (XREFPT,YREFPT,ZREFPT, XF, X0,Y0,Z0)

         IF (IX .EQ. NGSQ) THEN
C  distort towards middle of square - don't overlap border
            IF (X0 .LT. XMID) THEN
                X0 = X0 + NGDELT
            ELSE
                X0 = X0 - NGDELT
            ENDIF

            IF (Y0 .LT. YMID) THEN
                Y0 = Y0 + NGDELT
            ELSE
                Y0 = Y0 - NGDELT
            ENDIF
         ENDIF

C  outline annotation text
         CALL DRWREC (X0-XEXT/2, X0+XEXT/2, Y0-YBOT, Y0+YTOP)
260   CONTINUE

      CALL SETMSG ('7 8 9', 'The location of an annotation text ' //
     1             'primitive should reflect its 2D reference '   //
     2             'point after being transformed.')
      CALL DCHPFV ('EFFECT OF REFERENCE POINT: In which box is '  //
     1             'the annotation text NOT accurately outlined ' //
     2             'by the rectangle?', 6, NGSQ)
      CALL PEMST(103)

C  *** *** ***   effect of annotation offset   *** *** ***

C  ngsq = no-good square containing error = random from 1 to 6
      NGSQ = RNDINT (1,6)
C  chstr = "ABC"
      CHRABC = 'ABC'

      DO 300 IX = 1, 6
C  unscaled reference point
         XREFPT = U
         YREFPT = U
         ZREFPT = U

C  center of current square
         XMID = XWINLO(IX) + XSIZ/2
         YMID = YWINLO(IX) + YSIZ/2
C  scale into center of current square
         CALL ESC3 (XMID, YMID, H, XF)
         CALL PSLMT3 (XF, PCREPL)

C  pick NPC offset randomly - but keep within square
         XANNPT = NPCPWC * RNDRL ((XEXT-XSIZ)/2, (XSIZ-XEXT)/2)
         YANNPT = NPCPWC * RNDRL (YBOT - YSIZ/2, YSIZ/2 - YTOP)
         ZANNPT = NPCPWC * RNDRL (-0.499, 0.499)

         IF (IX .LE. 3) THEN
            CALL PATR3 (XREFPT,YREFPT,ZREFPT, XANNPT,YANNPT,ZANNPT,
     1                  CHRABC)
         ELSE
            CALL PATR  (XREFPT,YREFPT, XANNPT,YANNPT, CHRABC)
         ENDIF

         CALL PSLMT3 (IDM, PCREPL)
C  expected center of rectangle, in MC
         X0 = XMID + XANNPT/NPCPWC
         Y0 = YMID + YANNPT/NPCPWC

         IF (IX .EQ. NGSQ) THEN
C  distort towards middle of square - don't overlap border
            IF (X0 .LT. XMID) THEN
                X0 = X0 + NGDELT
            ELSE
                X0 = X0 - NGDELT
            ENDIF

            IF (Y0 .LT. YMID) THEN
                Y0 = Y0 + NGDELT
            ELSE
                Y0 = Y0 - NGDELT
            ENDIF
         ENDIF

         CALL DRWREC (X0-XEXT/2, X0+XEXT/2, Y0-YBOT, Y0+YTOP)
300   CONTINUE

      CALL SETMSG ('7 8 9', 'An annotation text primitive should be ' //
     1             'displaced from the reference point as specified ' //
     2             'by its annotation offset values.')
      CALL DCHPFV ('EFFECT OF ANNOTATION OFFSET: In which box is ' //
     1             'the annotation text NOT accurately outlined '  //
     2             'by the rectangle?', 6, NGSQ)

      CALL PEMST (103)

666   CONTINUE
C  wrap it up.
      CALL ENDIT
      END


text:
*> \brief \b SSXT1
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       REAL             FUNCTION SSXT1( IJOB, D1, N1, D2, N2, ABSTOL,
*                        ULP, UNFL )
*
*       .. Scalar Arguments ..
*       INTEGER            IJOB, N1, N2
*       REAL               ABSTOL, ULP, UNFL
*       ..
*       .. Array Arguments ..
*       REAL               D1( * ), D2( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SSXT1  computes the difference between a set of eigenvalues.
*> The result is returned as the function value.
*>
*> IJOB = 1:   Computes   max { min | D1(i)-D2(j) | }
*>                         i     j
*>
*> IJOB = 2:   Computes   max { min | D1(i)-D2(j) | /
*>                         i     j
*>                              ( ABSTOL + |D1(i)|*ULP ) }
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] IJOB
*> \verbatim
*>          IJOB is INTEGER
*>          Specifies the type of tests to be performed.  (See above.)
*> \endverbatim
*>
*> \param[in] D1
*> \verbatim
*>          D1 is REAL array, dimension (N1)
*>          The first array.  D1 should be in increasing order, i.e.,
*>          D1(j) <= D1(j+1).
*> \endverbatim
*>
*> \param[in] N1
*> \verbatim
*>          N1 is INTEGER
*>          The length of D1.
*> \endverbatim
*>
*> \param[in] D2
*> \verbatim
*>          D2 is REAL array, dimension (N2)
*>          The second array.  D2 should be in increasing order, i.e.,
*>          D2(j) <= D2(j+1).
*> \endverbatim
*>
*> \param[in] N2
*> \verbatim
*>          N2 is INTEGER
*>          The length of D2.
*> \endverbatim
*>
*> \param[in] ABSTOL
*> \verbatim
*>          ABSTOL is REAL
*>          The absolute tolerance, used as a measure of the error.
*> \endverbatim
*>
*> \param[in] ULP
*> \verbatim
*>          ULP is REAL
*>          Machine precision.
*> \endverbatim
*>
*> \param[in] UNFL
*> \verbatim
*>          UNFL is REAL
*>          The smallest positive number whose reciprocal does not
*>          overflow.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date December 2016
*
*> \ingroup single_eig
*
*  =====================================================================
      REAL             FUNCTION SSXT1( IJOB, D1, N1, D2, N2, ABSTOL,
     $                 ULP, UNFL )
*
*  -- LAPACK test routine (version 3.7.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     December 2016
*
*     .. Scalar Arguments ..
      INTEGER            IJOB, N1, N2
      REAL               ABSTOL, ULP, UNFL
*     ..
*     .. Array Arguments ..
      REAL               D1( * ), D2( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J
      REAL               TEMP1, TEMP2
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN
*     ..
*     .. Executable Statements ..
*
      TEMP1 = ZERO
*
      J = 1
      DO 20 I = 1, N1
   10    CONTINUE
         IF( D2( J ).LT.D1( I ) .AND. J.LT.N2 ) THEN
            J = J + 1
            GO TO 10
         END IF
         IF( J.EQ.1 ) THEN
            TEMP2 = ABS( D2( J )-D1( I ) )
            IF( IJOB.EQ.2 )
     $         TEMP2 = TEMP2 / MAX( UNFL, ABSTOL+ULP*ABS( D1( I ) ) )
         ELSE
            TEMP2 = MIN( ABS( D2( J )-D1( I ) ),
     $              ABS( D1( I )-D2( J-1 ) ) )
            IF( IJOB.EQ.2 )
     $         TEMP2 = TEMP2 / MAX( UNFL, ABSTOL+ULP*ABS( D1( I ) ) )
         END IF
         TEMP1 = MAX( TEMP1, TEMP2 )
   20 CONTINUE
*
      SSXT1 = TEMP1
      RETURN
*
*     End of SSXT1
*
      END

max_stars_repo_name: dnoan/OpenBLAS

max_issues_repo_name: dnoan/OpenBLAS

max_forks_repo_name: dnoan/OpenBLAS


text:
      program mpi_mod
        use mpi
      end program mpi_mod

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
      double precision function fdisc(blockno, x,y)
      implicit none

      integer blockno
      double precision x,y

      integer*8 cont, get_context
      logical fclaw2d_map_is_used

      double precision x0,y0,r0
      integer idisc
      common /cdisc/ x0,y0,r0
      common /comex/ idisc

      double precision r2,xp,yp,zp

      cont = get_context()

      IF (fclaw2d_map_is_used(cont)) THEN
         CALL fclaw2d_map_c2m(cont,blockno,x,y,xp,yp,zp)
         r2 = (xp-x0)**2 + (yp-y0)**2
      ELSE
         r2 = (x-x0)**2 + (y-y0)**2
      ENDIF

      fdisc = r2 - r0**2
c
      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
! { dg-do run }
! { dg-options "-fbounds-check" }
! { dg-shouldfail "Value 0 out of range in ORDER argument to RESHAPE intrinsic" }
program main
  implicit none
  integer(kind=4), dimension(6) :: source1 = (/ 1, 2, 3, 4, 5, 6 /)
  integer, dimension(2) :: shape1 = (/ 2, 3/)
  integer(kind=4), dimension(2) :: pad1 = (/ 0, 0/)
  character(len=200) :: l1, l2
  integer :: i1, i2

  l1 = "0 2"
  read(unit=l1,fmt=*) i1, i2
  write (unit=l2,fmt=*) reshape(source1, shape1, pad1, (/i1, i2/)) ! Invalid
end program main
! { dg-output "Fortran runtime error: Value 0 out of range in ORDER argument to RESHAPE intrinsic" }

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
! llarray.f90
module M
   implicit none
!   integer, parameter :: dp = kind([double precision ::])
   type PTRllnode
      type(llnode), pointer :: p => NULL()
   end type PTRllnode
   type llnode
      integer col
!      complex(dp) payload
      type(PTRllnode) next
   end type llnode
end module M

program Q
   use M
   implicit none
   type(PTRllnode), pointer :: p
   type(PTRllnode), allocatable, target :: heads(:)
   integer nrow
   integer row
   integer col
   integer iunit, istat
!   complex(dp) payload

   open(newunit=iunit,file='csriinput.txt',status='old')
   read(iunit,*) nrow
   allocate(heads(nrow))
   do
      read(iunit,*,iostat=istat) row,col !,payload
      if(istat /= 0 .OR. row < lbound(heads,1) .OR. row > ubound(heads,1)) exit
      p => heads(row)
      do while(associated(p%p))
         p => p%p%next
      end do
      allocate(p%p)
      p%p%col = col
!      p%p%payload = payload
   end do

   do row = lbound(heads,1),ubound(heads,1)
      print '(*(g0))', 'row: ',row
      p => heads(row)
      do while(associated(p%p))
         print '(*(g0))', '   col: ',p%p%col
         p => p%p%next
      end do
   end do
end program Q

max_stars_repo_name: MorrisHuang-skipper/Serial-MD

max_issues_repo_name: MorrisHuang-skipper/Serial-MD

max_forks_repo_name: MorrisHuang-skipper/Serial-MD


program example_pointers

implicit none (type, external)

integer :: N=5
integer, pointer :: ptr_N
integer, pointer :: ptr2_N

allocate(ptr_N)
ptr_N=10
ptr2_N=>ptr_N
print*, ptr_N
print*, ptr2_N
print*, N

end program example_pointers


text: The Consortium for Women and Research at UCD is dedicated to the support of research by and on women and on gender in its multiple intersections with race, class, sexual identity, and other categories of identity and analysis.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE CALHTA
C
C CHANGE RECORD
C **  SUBROUTINE CALHTA PERFORMS A HARMONIC ANALYSIS FOR THE M2 TIDE
C **  OVER TWO TIDAL CYCLES
C
      USE GLOBAL
      CHARACTER*80 TITLE1,TITLE2,TITLE3,TITLE4,TITLE11,TITLE12
C
C **  INITIALIZE ON FIRST ENTRY FOR CURRENT ANALYSIS INTERVAL
C
      IF(NHAR.GT.1) GOTO 1000
      DO L=2,LA
        AMCP(L)=0.
        AMSP(L)=0.
        AMCUE(L)=0.
        AMSUE(L)=0.
        AMCVE(L)=0.
        AMSVE(L)=0.
      ENDDO
      DO K=1,KC
        DO L=2,LA
          AMCU(L,K)=0.
          AMSU(L,K)=0.
          AMCV(L,K)=0.
          AMSV(L,K)=0.
        ENDDO
      ENDDO
      OPEN(1,FILE='SURFAMP.OUT',STATUS='UNKNOWN')
      CLOSE(1,STATUS='DELETE')
      OPEN(1,FILE='SURFAMP.OUT',STATUS='UNKNOWN')
      TITLE1='SURFACE DISPLACEMENT AMPLITUDE CONTOURS'
      OPEN(2,FILE='SURFPHA.OUT',STATUS='UNKNOWN')
      CLOSE(2,STATUS='DELETE')
      OPEN(2,FILE='SURFPHA.OUT',STATUS='UNKNOWN')
      TITLE2='SURFACE DISPLACEMENT PHASE CONTOURS'
      OPEN(3,FILE='MAJAXIS.OUT',STATUS='UNKNOWN')
      CLOSE(3,STATUS='DELETE')
      OPEN(3,FILE='MAJAXIS.OUT',STATUS='UNKNOWN')
      TITLE3='MAJOR AXES OF TIDAL ELLIPSES'
      OPEN(4,FILE='MAJAPHA.OUT',STATUS='UNKNOWN')
      CLOSE(4,STATUS='DELETE')
      OPEN(4,FILE='MAJAPHA.OUT',STATUS='UNKNOWN')
      TITLE4='TIDAL ELLIPSES PHASE CONTOURS'
      OPEN(11,FILE='TIDELKC.OUT',STATUS='UNKNOWN')
      CLOSE(11,STATUS='DELETE')
      OPEN(11,FILE='TIDELKC.OUT',STATUS='UNKNOWN')
      TITLE11='SURFACE TIDAL ELLIPSES'
      OPEN(12,FILE='TIDELKB.OUT',STATUS='UNKNOWN')
      CLOSE(12,STATUS='DELETE')
      OPEN(12,FILE='TIDELKB.OUT',STATUS='UNKNOWN')
      TITLE12='BOTTOM TIDAL ELLIPSES'
      LINES=LA-1
      LEVELS=1
      DBS=0.
      WRITE (1,99) TITLE1
      WRITE (1,101)LINES,LEVELS
      WRITE (1,250)DBS
      CLOSE(1)
      WRITE (2,99) TITLE2
      WRITE (2,101)LINES,LEVELS
      WRITE (2,250)DBS
      CLOSE(2)
      WRITE (11,99) TITLE11
      WRITE (11,101)LINES,LEVELS
      WRITE (11,250)DBS
      CLOSE(11)
      DBS=99.
      WRITE (12,99) TITLE12
      WRITE (12,101)LINES,LEVELS
      WRITE (12,250)DBS
      CLOSE(12)
      LEVELS=2
      DBS1=0.
      DBS2=99.
      WRITE (3,99) TITLE3
      WRITE (3,101)LINES,LEVELS
      WRITE (3,250)DBS1,DBS2
      CLOSE(3)
      WRITE (4,99) TITLE4
      WRITE (4,101)LINES,LEVELS
      WRITE (4,250)DBS1,DBS2
      CLOSE(4)
C
C **  ACCUMULATE HARMONIC ANALYSIS
C
 1000 CONTINUE
      DO L=2,LA
        LN=LNC(L)
        AMCP(L)=P(L)*WC(NHAR)+AMCP(L)
        AMSP(L)=P(L)*WS(NHAR)+AMSP(L)
        UTMP1=0.5*(UHDYE(LEAST(L))+UHDYE(L))*HPI(L)/DYP(L)
        VTMP1=0.5*(VHDXE(LN)+VHDXE(L))*HPI(L)/DXP(L)
        UTMP=CUE(L)*UTMP1+CVE(L)*VTMP1
        VTMP=CUN(L)*UTMP1+CVN(L)*VTMP1
        AMCUE(L)=UTMP*WC(NHAR)+AMCUE(L)
        AMSUE(L)=UTMP*WS(NHAR)+AMSUE(L)
        AMCVE(L)=VTMP*WC(NHAR)+AMCVE(L)
        AMSVE(L)=VTMP*WS(NHAR)+AMSVE(L)
      ENDDO
      DO K=1,KC
        DO L=2,LA
          LN=LNC(L)
          UTMP1=0.5*(U(LEAST(L),K)+U(L,K))
          VTMP1=0.5*(V(LN,K)+V(L,K))
          UTMP=CUE(L)*UTMP1+CVE(L)*VTMP1
          VTMP=CUN(L)*UTMP1+CVN(L)*VTMP1
          AMCU(L,K)=UTMP*WC(NHAR)+AMCU(L,K)
          AMSU(L,K)=UTMP*WS(NHAR)+AMSU(L,K)
          AMCV(L,K)=VTMP*WC(NHAR)+AMCV(L,K)
          AMSV(L,K)=VTMP*WS(NHAR)+AMSV(L,K)
        ENDDO
      ENDDO
C
C **  CHECK FOR END OF ANALYSIS INTERVAL
C
      IF(NHAR.LT.NTSPTC2) GOTO 2000
C
C **  COMPLETE HARMONIC ANALYSIS
C
      DO L=2,LA
        AMC=AS*AMCP(L)-ACS*AMSP(L)
        AMS=-ACS*AMCP(L)+AC*AMSP(L)
        AMCP(L)=AMC
        AMSP(L)=AMS
        AMC=AS*AMCUE(L)-ACS*AMSUE(L)
        AMS=-ACS*AMCUE(L)+AC*AMSUE(L)
        AMCUE(L)=AMC
        AMSUE(L)=AMS
        AMC=AS*AMCVE(L)-ACS*AMSVE(L)
        AMS=-ACS*AMCVE(L)+AC*AMSVE(L)
        AMCVE(L)=AMC
        AMSVE(L)=AMS
      ENDDO
      DO K=1,KC
        DO L=2,LA
          AMC=AS*AMCU(L,K)-ACS*AMSU(L,K)
          AMS=-ACS*AMCU(L,K)+AC*AMSU(L,K)
          AMCU(L,K)=AMC
          AMSU(L,K)=AMS
          AMC=AS*AMCV(L,K)-ACS*AMSV(L,K)
          AMS=-ACS*AMCV(L,K)+AC*AMSV(L,K)
          AMCV(L,K)=AMC
          AMSV(L,K)=AMS
        ENDDO
      ENDDO
      NHAR=0
      OPEN(1,FILE='SURFAMP.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (1,100)N
      OPEN(2,FILE='SURFPHA.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (2,100)N
      DO L=2,LA
        SSURFAMP=GI*SQRT(AMCP(L)*AMCP(L)+AMSP(L)*AMSP(L))
        IF(AMCP(L).EQ.0.0.AND.AMSP(L).EQ.0.0)THEN
          PHI=999999.
        ELSE
          PHI=ATAN2(AMSP(L),AMCP(L))
        ENDIF
        SSURFPHS=TIDALP*PHI/(3600.*PI2)
        SSURFPSC=TIDALP*PHI/PI2
        IF(SSURFPSC.LT.0.0)SSURFPSC=SSURFPSC+TIDALP
        WRITE(1,200)IL(L),JL(L),DLON(L),DLAT(L),SSURFAMP,SSURFPSC
        WRITE(2,200)IL(L),JL(L),DLON(L),DLAT(L),SSURFPHS,SSURFPSC
      ENDDO
      CLOSE(1)
      CLOSE(2)
      OPEN(3,FILE='MAJAXIS.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (3,100)N
      OPEN(4,FILE='MAJAPHA.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (4,100)N
      OPEN(11,FILE='TIDELKC.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (11,100)N
      OPEN(12,FILE='TIDELKB.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (12,100)N
      DO L=2,LA
        TERM1=AMCU(L,KC)+AMSV(L,KC)
        TERM2=AMCV(L,KC)-AMSU(L,KC)
        TERM3=AMCU(L,KC)-AMSV(L,KC)
        TERM4=AMCV(L,KC)+AMSU(L,KC)
        RPLUS=0.5*SQRT(TERM1*TERM1+TERM2*TERM2)
        RMINS=0.5*SQRT(TERM3*TERM3+TERM4*TERM4)
        IF(TERM1.EQ.0.0.AND.TERM2.EQ.0.0)THEN
          APLUS=999999.
        ELSE
          APLUS=ATAN2(TERM2,TERM1)
        ENDIF
        IF(TERM3.EQ.0.0.AND.TERM4.EQ.0.0)THEN
          AMINS=999999.
        ELSE
          AMINS=ATAN2(TERM4,TERM3)
        ENDIF
        RRMAJ=RPLUS+RMINS
        RRMIN=ABS(RPLUS-RMINS)
        AACCWX=0.5*(APLUS+AMINS)
        RMAJUKC=RRMAJ*COS(AACCWX)
        RMAJVKC=RRMAJ*SIN(AACCWX)
        IF(RMAJUKC.LT.0.0)THEN
          RMAJUKC=-RMAJUKC
          RMAJVKC=-RMAJVKC
        ENDIF
        PHASEKC=(0.25/PI)*TIDALP*(AMINS-APLUS)/3600.
        WRITE(11,200)IL(L),JL(L),DLON(L),DLAT(L),AACCWX,RRMAJ,RRMIN
        TERM1=AMCU(L,1)+AMSV(L,1)
        TERM2=AMCV(L,1)-AMSU(L,1)
        TERM3=AMCU(L,1)-AMSV(L,1)
        TERM4=AMCV(L,1)+AMSU(L,1)
        RPLUS=0.5*SQRT(TERM1*TERM1+TERM2*TERM2)
        RMINS=0.5*SQRT(TERM3*TERM3+TERM4*TERM4)
        IF(TERM1.EQ.0.0.AND.TERM2.EQ.0.0)THEN
          APLUS=999999.
        ELSE
          APLUS=ATAN2(TERM2,TERM1)
        ENDIF
        IF(TERM3.EQ.0.0.AND.TERM4.EQ.0.0)THEN
          AMINS=999999.
        ELSE
          AMINS=ATAN2(TERM4,TERM3)
        ENDIF
        RRMAJ=RPLUS+RMINS
        RRMIN=ABS(RPLUS-RMINS)
        AACCWX=0.5*(APLUS+AMINS)
        RMAJUKB=RRMAJ*COS(AACCWX)
        RMAJVKB=RRMAJ*SIN(AACCWX)
        IF(RMAJUKB.LT.0.0)THEN
          RMAJUKB=-RMAJUKB
          RMAJVKB=-RMAJVKB
        ENDIF
        PHASEKB=(0.25/PI)*TIDALP*(AMINS-APLUS)/3600.
        WRITE(3,200)IL(L),JL(L),DLON(L),DLAT(L),RMAJUKC,RMAJVKC,
     &      RMAJUKB,RMAJVKB
        WRITE(4,200)IL(L),JL(L),DLON(L),DLAT(L),PHASEKC,PHASEKB
        WRITE(12,200)IL(L),JL(L),DLON(L),DLAT(L),AACCWX,RRMAJ,RRMIN
      ENDDO
      CLOSE(3)
      CLOSE(4)
      CLOSE(11)
      CLOSE(12)
 2000 CONTINUE
      NHAR=NHAR+1
   99 FORMAT(A80)
  100 FORMAT(I10)
  101 FORMAT(2I10)
  200 FORMAT(2I5,1X,6E13.5)
  250 FORMAT(12E12.4)
      RETURN
      END


text:
program foo
  real :: r4_res, r4_x = 2, r4_y = 10
  real(kind = 8) :: r8_res, r8_x = 2.5, r8_y = 3

  integer :: i4_res, i4_x = 3, i4_y = 12

  r4_res = r4_x ** r4_y
  r8_res = r8_x ** r8_y
  print *, r4_res
  print *, r8_res

  r4_res = r4_x ** i4_y
  print *, r4_res
end program foo

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
module model_temperature_gradient
#include <messenger.h>
    use mod_kinds,          only: rk
    use mod_constants,      only: HALF, ONE, TWO
    use mod_fluid,          only: gam, Rgas
    use type_model,         only: model_t
    use type_chidg_worker,  only: chidg_worker_t
    use DNAD_D
    implicit none

    !>  A model for computing shear stress
    !!
    !!  Model Fields:
    !!      : shear_11, shear_12, shear_13
    !!      : shear_21, shear_22, shear_23
    !!      : shear_31, shear_32, shear_33
    !!
    !!  @author Nathan A. Wukie
    !!  @date   12/1/2016
    !!
    !---------------------------------------------------------------------------------------
    type, extends(model_t)  :: temperature_gradient_t

    contains

        procedure   :: init
        procedure   :: compute

    end type temperature_gradient_t
    !***************************************************************************************

contains

    !>  Initialize the model with a name and the model fields it is contributing to.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   12/1/2016
    !!
    !---------------------------------------------------------------------------------------
    subroutine init(self)
        class(temperature_gradient_t), intent(inout)   :: self

        call self%set_name('Temperature Gradient')
        call self%set_dependency('f(Grad(Q))')

        call self%add_model_field('Temperature Gradient - 1')
        call self%add_model_field('Temperature Gradient - 2')
        call self%add_model_field('Temperature Gradient - 3')

    end subroutine init
    !***************************************************************************************

    !>  Routine for computing the Pressure and Temperature.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   12/1/2016
    !!
    !--------------------------------------------------------------------------------------
    subroutine compute(self,worker)
        class(temperature_gradient_t),  intent(in)      :: self
        type(chidg_worker_t),           intent(inout)   :: worker

        type(AD_D), allocatable,    dimension(:) ::         &
            density, mom1, mom2, mom3, energy, invdensity,  &
            grad1_density, grad2_density, grad3_density,    &
            grad1_mom1,    grad2_mom1,    grad3_mom1,       &
            grad1_mom2,    grad2_mom2,    grad3_mom2,       &
            grad1_mom3,    grad2_mom3,    grad3_mom3,       &
            grad1_energy,  grad2_energy,  grad3_energy,     &
            grad1_T,       grad2_T,       grad3_T,          &
            u,             v,             w,       p,       &
            dT_ddensity,   dT_dmom1,   dT_dmom2,   dT_dmom3,   dT_denergy,  &
            dp_ddensity,   dp_dmom1,   dp_dmom2,   dp_dmom3,   dp_denergy,  &
            dke_ddensity,  dke_dmom1,  dke_dmom2,  dke_dmom3, r

        type(AD_D), allocatable, dimension(:,:) :: grad_density, grad_mom1, grad_mom2, grad_mom3, grad_energy

        real(rk)    :: const

        !
        ! Interpolate solution to quadrature nodes
        !
        density = worker%get_field('Density',    'value')
        mom1    = worker%get_field('Momentum-1', 'value')
        mom2    = worker%get_field('Momentum-2', 'value')
        mom3    = worker%get_field('Momentum-3', 'value')
        energy  = worker%get_field('Energy',     'value')

        !
        ! Interpolate gradient to quadrature nodes
        !
        grad1_density = worker%get_field('Density'   , 'grad1')
        grad2_density = worker%get_field('Density'   , 'grad2')
        grad3_density = worker%get_field('Density'   , 'grad3')

        grad1_mom1    = worker%get_field('Momentum-1', 'grad1')
        grad2_mom1    = worker%get_field('Momentum-1', 'grad2')
        grad3_mom1    = worker%get_field('Momentum-1', 'grad3')

        grad1_mom2    = worker%get_field('Momentum-2', 'grad1')
        grad2_mom2    = worker%get_field('Momentum-2', 'grad2')
        grad3_mom2    = worker%get_field('Momentum-2', 'grad3')

        grad1_mom3    = worker%get_field('Momentum-3', 'grad1')
        grad2_mom3    = worker%get_field('Momentum-3', 'grad2')
        grad3_mom3    = worker%get_field('Momentum-3', 'grad3')

        grad1_energy  = worker%get_field('Energy    ', 'grad1')
        grad2_energy  = worker%get_field('Energy    ', 'grad2')
        grad3_energy  = worker%get_field('Energy    ', 'grad3')

        !
        ! Account for cylindrical. Get tangential momentum from angular momentum.
        ! Also convert derivatives from derivatives of angular momentum to tangential.
        !
        ! We want:
        !       (rho * u_theta)  instead of      (r * rho * u_theta)
        !   grad(rho * u_theta)  instead of  grad(r * rho * u_theta)
        !
        !   grad(rho * u_theta) = grad(r * rho * u_theta)/r  -  grad(r)(rho*u_theta)/r
        !
        ! Where grad(r) = [1,0,0]
        !
        if (worker%coordinate_system() == 'Cylindrical') then
            r = worker%coordinate('1')
            mom2       = mom2 / r
            grad1_mom2 = (grad1_mom2/r) - mom2/r
            grad2_mom2 = (grad2_mom2/r)
            grad3_mom2 = (grad3_mom2/r)
        end if

        !
        ! Get model field 'Pressure'
        !
        p = worker%get_field('Pressure', 'value')

        !
        ! Compute velocities
        !
        invdensity = ONE/density
        u          = mom1*invdensity
        v          = mom2*invdensity
        w          = mom3*invdensity

        !
        ! Compute Kinetic Energy Jacobians
        !
        dke_ddensity = -HALF*(u*u + v*v + w*w)
        dke_dmom1    = u
        dke_dmom2    = v
        dke_dmom3    = w

        !
        ! Compute Pressure Jacobians
        !
        dp_ddensity = -(gam-ONE)*dke_ddensity
        dp_dmom1    = -(gam-ONE)*dke_dmom1
        dp_dmom2    = -(gam-ONE)*dke_dmom2
        dp_dmom3    = -(gam-ONE)*dke_dmom3
        dp_denergy  =  dp_dmom3    ! Initialize derivatives
        dp_denergy  =  (gam-ONE)   ! No negative sign

        !
        ! Compute Temperature Jacobians
        !
        const = ONE/Rgas
        dT_ddensity = const*invdensity*dp_ddensity  -  const*invdensity*invdensity*p
        dT_dmom1    = const*invdensity*dp_dmom1
        dT_dmom2    = const*invdensity*dp_dmom2
        dT_dmom3    = const*invdensity*dp_dmom3
        dT_denergy  = const*invdensity*dp_denergy

        !
        ! Compute temperature gradient
        !
        grad1_T = dT_ddensity*grad1_density + dT_dmom1*grad1_mom1 + dT_dmom2*grad1_mom2 + dT_dmom3*grad1_mom3 + dT_denergy*grad1_energy
        grad2_T = dT_ddensity*grad2_density + dT_dmom1*grad2_mom1 + dT_dmom2*grad2_mom2 + dT_dmom3*grad2_mom3 + dT_denergy*grad2_energy
        grad3_T = dT_ddensity*grad3_density + dT_dmom1*grad3_mom1 + dT_dmom2*grad3_mom2 + dT_dmom3*grad3_mom3 + dT_denergy*grad3_energy

        call worker%store_model_field('Temperature Gradient - 1', 'value', grad1_T)
        call worker%store_model_field('Temperature Gradient - 2', 'value', grad2_T)
        call worker%store_model_field('Temperature Gradient - 3', 'value', grad3_T)

    end subroutine compute
    !***************************************************************************************

end module model_temperature_gradient

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
	SUBROUTINE G2T_WDIR1 ( itrnd, kflag, wd_c, iret )
C************************************************************************
C* G2T_WDIR1								*
C*									*
C* This subroutine tries to combine the wind directions from the RANGE	*
C* and EXCEPT parts.  For example, if the wind direction is N and NW	*
C* over RANGE and EXCEPT parts respectively, it may be combined to N 	*
C* to NW under speed constraint.					*
C*									*
C* G2T_WDIR1 ( ITRND, KFLAG, WD_C, IRET )				*
C*									*
C* Input parameters:							*
C*	ITRND		INTEGER		Index for trending block	*
C*									*
C* Output parameters:							*
C*	KFLAG		LOGICAL		Combine wind dir flag		*
C*	WD_C(2)		CHAR*		Combined wind dir		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* T. Lee/SAIC		06/07	Created					*
C************************************************************************
	INCLUDE		'goftxt.cmn'
	CHARACTER*(*)	wd_c(*)
	LOGICAL		eflag, kflag
	CHARACTER	wd_r(2)*2, wd_e(2)*2
C------------------------------------------------------------------------
	iret = 0
	DO ii = 1, 2
	    wd_c ( ii ) = ' '
	END DO
	kflag = .false.
C
C*	If no EXCEPT part, no need to combine wind directions.
C
	eflag = eflag_d ( 2, itrnd )
	IF ( .not. eflag )  RETURN
C
C*	Convert compass direction to wind direction.
C
	DO ii = 1, 2
	    wd_r ( ii ) =  wdir_d ( ii, itrnd )
	    wd_e ( ii ) =  wdir_de ( ii, itrnd )
	END DO
C
	CALL G2T_WDIR ( wd_r, wd_e, kflag, wd_c, ier )
C
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
cc Copyright (C) 2009: Zydrunas Gimbutas and Hong Xiao
cc Contact: Zydrunas Gimbutas <gimbutas@cims.nyu.edu>
cc          Hong Xiao <hxiao@ucdavis.edu>
cc
cc This software is being released under a modified FreeBSD license
cc (see COPYING in home directory).
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c       Orthogonal polynomials on the standard triangle
c
c       (0,2/sqrt(3)), (-1,-1/sqrt(3)), (1,-1/sqrt(3))
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c       This file contains a set of subroutines for the handling
c       of the evaluation of orthogonal polynomials on the standard
c       triangle.  It contains 2 user-callable subroutines.
c
c   ortho2eva - evaluates at the user-supplied point z
c       a collection of polynomials (of x,y) orthogonal on the
c       standard triangle
c
c   ortho2eva3 - evaluates at the user-supplied point z a collection
c       of polynomials (of x,y) orthogonal on the standard triangle,
c       together with their derivatives with respect to x and y.

c
c
      subroutine ortho2eva(mmax,z,pols,w)
c
c     This subroutine evaluates at the user-supplied point z
c     a collection of polynomials (of x,y) orthogonal on the
c     standard triangle.
c
c     The "standard" triangle is the triangle with the
c     vertices
c
c     (0,2/sqrt(3)), (-1,-1/sqrt(3)), (1,-1/sqrt(3)).       (1)
c
c     The polynomials evaluated by this subroutine are all
c     orthogonal polynomials up to order mmax, arranged in the
c     increasing order.
c
c     This subroutine is based on the Koornwinder's representation
c     of the orthogonal polynomials on the right triangle
c
c     (-1,-1), (-1,1), (1,-1)                                       (2)
c
c     given by
c
c     K_mn(x,y) = P_m ((2*x+1+y)/(1-y)) ((1-y)/2)^m P_n^{2m+1,0} (y)
c
c     where P_m are the Legendre polynomials or order m
c     and P_n^{2m+1,0} are the Jacobi polynomials of order n with
c     the parameters alpha=2*m+1 and beta=0.
c
c                   Input parameters:
c
c  mmax - the maximum order to which the polynomials are to be evaluated;
c  z - the location in R^2 where the polynomials are to be evaluated;
c       normally, expected to be inside (including boundary) the
c       standard triangle (1) above.
c  w - the work array. must be at least mmax+1+(mmax+1)**2 long
c
c                   Output parameters:
c
c  pols - the orthogonal polynomials evaluated at the point z
c       ( (mmax+1)*(mmax+2)/2 of them things)
c
c       . . . allocate memory for the work arrays
c
      implicit real *8 (a-h,o-z)
      dimension z(2), pols(1), w(1)
      data init/1/
      save
c
        if( init .eq. 1 ) then
        c0=1/sqrt(3.0d0)*sqrt(sqrt(3.0d0))
        c1=sqrt(2.0d0)*sqrt(sqrt(3.0d0))
        c2=sqrt(2.0d0)*sqrt(sqrt(3.0d0))
        init=0
        endif
c
        if( mmax .eq. 0 ) then
        pols(1)=c0
        return
        endif
c
        if( mmax .eq. 1 ) then
        pols(1)=c0
        pols(2)=z(1)*c1
        pols(3)=z(2)*c2
        return
        endif
c
      iw1=1
      iw2=iw1+mmax+1
      iw3=iw2+(mmax+1)**2
      call ortho2eva0(mmax,z,pols,w(iw1),w(iw2))
      return
      end
c
c
c
c
c
      subroutine ortho2eva0(mmax,z,pols,f1,f2)
c
c     evaluate the orthonormal polynomials on the triangle
c
      implicit real *8 (a-h,o-z)
      dimension z(2), pols(1), f1(1), f2(1)
      data init/1/
      save
c
      if( init.eq.1 ) then
         zero=0
         sqrt2=dsqrt(2.0d0)
         sqrt3=dsqrt(3.0d0)
         init=0
         r11=-1.0d0/3.0d0
         r12=-1.0d0/sqrt3
         r21=1.0d0/9.0d0*(-sqrt3)*sqrt3
         r22=1.0d0/9.0d0*(6.0d0)*sqrt3
      endif
c
      a=z(1)
      b=z(2)
c
c     ... map the standard triangle to the right
c     triangle with the vertices (-1,-1), (1,-1), (-1,1)
c
      x=r11+r12*b+a
      y=r21+r22*b
c
c     ... evaluate the Koornwinder's polynomials
c     the via three terms recursion
c
      par1=(2*x+1+y)/2
      par2=(1-y)/2
      call klegeypols(par1,par2,mmax,f1)
      do 1200 m=0,mmax
         par1=2*m+1
         call kjacopols(y,par1,zero,(mmax-m),f2(1+m*(mmax+1)))
 1200 continue

      kk=0
      do 2200 m=0,mmax
         do 2000 n=0,m
         kk=kk+1
c
c     ... evaluate the polynomial (m-n, n)
c
         pols(kk)=f1(m-n+1)*f2(n+1+(m-n)*(mmax+1))
c
c     ... and normalize it
c
         scale=dsqrt((1.0d0+(m-n)+n)*(1.0d0+(m-n)+(m-n))/sqrt3)
cccc         scale=1
         pols(kk)=pols(kk)*scale
 2000	 continue
 2200 continue
      return
      end
c
c
c
c
c
      subroutine ortho2eva3(mmax,z,pols,dersx,dersy,w)
c
c     This subroutine evaluates at the user-supplied point z
c     a collection of polynomials (of x,y) orthogonal on the
c     standard triangle, together with their derivatives with
c     respect to x and y.
c
c     The "standard" triangle is the triangle with the
c     vertices
c
c     (0,2/sqrt(3)), (-1,-1/sqrt(3), (1,-1/sqrt(3)).       (1)
c
c     The polynomials evaluated by this subroutine are all
c     orthogonal polynomials up to order mmax, arranged in the
c     increasing order.
c
c     This subroutine is based on the Koornwinder's representation
c     of the orthogonal polynomials on the right triangle
c
c     (-1,-1), (-1,1), (1,-1)                                       (2)
c
c     given by
c
c     K_mn(x,y) = P_m ((2*x+1+y)/(1-y)) ((1-y)/2)^m P_n^{2m+1,0} (y)
c
c     where P_m are the Legendre polynomials or order m
c     and P_n^{2m+1,0} are the Jacobi polynomials of order n with
c     the parameters alpha=2*m+1 and beta=0.
c
c                   Input parameters:
c
c  mmax - the maximum order to which the polynomials are to be evaluated;
c  z - the location in R^2 where the polynomials are to be evaluated;
c       normally, expected to be inside (including boundary) the
c       standard triangle (1) above.
c  w - the work array. must be at least mmax+1+(mmax+1)**2 long
c
c                   Output parameters:
c
c  pols - the orthogonal polynomials evaluated at the point z
c       ( (mmax+1)*(mmax+2)/2 of them things)
c  dersx - the derivatives with respect to x of the polynomials
c       returned in array pols
c  dersy - the derivatives with respect to y of the polynomials
c       returned in array pols
c
c       . . . allocate memory for the work arrays
c
      implicit real *8 (a-h,o-z)
      dimension z(2),pols(1),dersx(1),dersy(1),w(1)
      data init/1/
      save
c
        if( init .eq. 1 ) then
        c0=1/sqrt(3.0d0)*sqrt(sqrt(3.0d0))
        c1=sqrt(2.0d0)*sqrt(sqrt(3.0d0))
        c2=sqrt(2.0d0)*sqrt(sqrt(3.0d0))
        init=0
        endif
c
        if( mmax .eq. 0 ) then
        pols(1)=c0
        dersx(1)=0
        dersy(1)=0
        return
        endif
c
        if( mmax .eq. 1 ) then
        pols(1)=c0
        pols(2)=z(1)*c1
        pols(3)=z(2)*c2
        dersx(1)=0
        dersx(2)=c1
        dersx(3)=0
        dersy(1)=0
        dersy(2)=0
        dersy(3)=c2
        return
        endif
c
      iw1=1
      iw2=iw1+mmax+1
      iw3=iw2+(mmax+1)**2
      iw4=iw3+mmax+1
      iw5=iw4+mmax+1
      iw6=iw5+(mmax+1)**2
      iw7=iw6+(mmax+1)**2
      call ortho2eva30(mmax,z,pols,dersx,dersy,
     $     w(iw1),w(iw2),w(iw3),w(iw4),w(iw5),w(iw6))
      return
      end
c
c
c
c
c
      subroutine ortho2eva30(mmax,z,pols,dersx,dersy,
     $     f1,f2,f3,f4,f5,f6)
c
c     evaluate the orthonormal polynomials on the triangle,
c     together with their derivatives with respect to x and y.
c
      implicit real *8 (a-h,o-z)
      dimension z(2), pols(1), dersx(1), dersy(1)
      dimension f1(1), f2(1)
      dimension f3(1), f4(1), f5(1), f6(1)
      data init/1/
      save
c
      if( init.eq.1 ) then
         zero=0
         sqrt2=dsqrt(2.0d0)
         sqrt3=dsqrt(3.0d0)
         init=0
         r11=-1.0d0/3.0d0
         r12=-1.0d0/sqrt3
         r21=1.0d0/9.0d0*(-sqrt3)*sqrt3
         r22=1.0d0/9.0d0*(6.0d0)*sqrt3
      endif
c
      x=z(1)
      y=z(2)
c
c     ... map the standard triangle to the right
c     triangle with the vertices (-1,-1), (1,-1), (-1,1)
c
      a=r11+r12*y+x
      b=r21+r22*y

c
c     ... evaluate the Koornwinder's polynomials
c     the via three terms recursion
c
      par1=(2*a+1+b)/2
      par2=(1-b)/2
      call klegeypols3(par1,par2,mmax,f1,f3,f4)

      do 1200 m=0,mmax
         par1=2*m+1
         call kjacopols2(
     $        b,par1,zero,(mmax-m),f2(1+m*(mmax+1)),f5(1+m*(mmax+1)))
 1200 continue

      kk=0
      do 2200 m=0,mmax
         do 2000 n=0,m
         kk=kk+1
c
c     ... evaluate the polynomial (m-n, n), and their derivatives
c     with respect to x,y
c
         pols(kk)=f1(m-n+1)*f2(n+1+(m-n)*(mmax+1))
c
         dersx(kk)=
     $        f3(m-n+1)*f2(n+1+(m-n)*(mmax+1))
c
         dersy(kk)=
     $        f1(m-n+1)*f5(n+1+(m-n)*(mmax+1))*r22 +
     $        f3(m-n+1)*f2(n+1+(m-n)*(mmax+1))*(r12+r22/2) +
     $        f4(m-n+1)*f2(n+1+(m-n)*(mmax+1))*(-r22/2)

c
c     ... and normalize it
c
         scale=dsqrt((1.0d0+(m-n)+n)*(1.0d0+(m-n)+(m-n))/sqrt3)
cccc         scale=1
         pols(kk)=pols(kk)*scale
         dersx(kk)=dersx(kk)*scale
         dersy(kk)=dersy(kk)*scale
 2000	 continue
 2200 continue

      return
      end
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        this is the end of the debugging code and the beginning
c        of the evaluation of scaled Legendre polynomials
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
        subroutine klegeypols(x,y,n,pols)
        implicit real *8 (a-h,o-z)
        dimension pols(1)
c
c       evaluate a sequence of scaled Legendre polynomials
c       P_n(x/y) y^n, with the parameter y \in [0..1]
c
c       ...
c
        pkm1=1
        pk=x
c
        pk=1
        pkp1=x
c
c
c        if n=0 or n=1 - exit
c
        if(n .ge. 2) goto 1200
        pols(1)=1
        if(n .eq. 0) return
c
        pols(2)=x
        return
 1200 continue
c
        pols(1)=1
        pols(2)=x
c
c       n is greater than 2. conduct recursion
c
        do 2000 k=1,n-1
c
        pkm1=pk
        pk=pkp1
        pkp1=( (2*k+1)*x*pk-k*pkm1*y**2 )/(k+1)
        pols(k+2)=pkp1
 2000 continue
c
        return
        end
c
c
c
c
c
        subroutine klegeypols3(x,y,n,pols,dersx,dersy)
        implicit real *8 (a-h,o-z)
        dimension pols(1),dersx(1),dersy(1)
c
c     evaluate a sequence of scaled Legendre polynomials P_n(x/y) y^n,
c     with the parameter y \in [0..1], together with their derivatives
c     with respect to the parameters x and y.
c
c       ... if n=0 or n=1 - exit
c
        if(n .ge. 2) goto 1200
        pols(1)=1
        dersx(1)=0
        dersy(1)=0
        if(n .eq. 0) return
c
        pols(2)=x
        dersx(2)=1
        dersy(2)=0
        return
 1200 continue
c
c       n is greater than 1. conduct recursion
c
        pols(1)=1
        dersx(1)=0
        dersy(1)=0
        pols(2)=x
        dersx(2)=1
        dersy(2)=0
c
        pkm1=1
        pk=x
        dkm1=0
        dk=1
        ykm1=0
        yk=0
c
        pk=1
        pkp1=x
        dk=0
        dkp1=1
        yk=0
        ykp1=0
c
        do 2000 k=1,n-1
c
        pkm1=pk
        pk=pkp1
        dkm1=dk
        dk=dkp1
        ykm1=yk
        yk=ykp1
c
        pkp1=( (2*k+1)*x*pk-k*pkm1*y*y )/(k+1)
        dkp1=( (2*k+1)*(x*dk+pk)-k*dkm1*y*y )/(k+1)
        ykp1=( (2*k+1)*(x*yk)-k*(pkm1*2*y+ykm1*y*y) )/(k+1)
        pols(k+2)=pkp1
        dersx(k+2)=dkp1
        dersy(k+2)=ykp1
 2000 continue
c
        return
        end
c
c
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c        this is the end of the debugging code and the beginning
c        of the evaluation of Jacobi polynomials
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
c
        subroutine kjacopols(x,a,b,n,pols)
        implicit real *8 (a-h,o-z)
        dimension pols(1)
c
c       evaluates a bunch of Jacobi polynomials
c       at the user-provided point x
c
c       ... if n=0 or n=1 - exit
c
        if(n .ge. 2) goto 1200
        pols(1)=1
        if(n .eq. 0) return
c
        pols(2)=(a/2-b/2)+(1+a/2+b/2)*x
        return
 1200 continue
c
        pols(1)=1
        pols(2)=(a/2-b/2)+(1+a/2+b/2)*x
c
c       n is greater than 2. conduct recursion
c
        pkm1=1
        pk=(a/2-b/2)+(1+a/2+b/2)*x
        pk=1
        pkp1=(a/2-b/2)+(1+a/2+b/2)*x
c
        do 2000 k=2,n
c
        pkm1=pk
        pk=pkp1
        alpha=(2*k+a+b-1)*(a**2-b**2+(2*k+a+b-2)*(2*k+a+b)*x)
        beta=2*(k+a-1)*(k+b-1)*(2*k+a+b)
        pkp1=(alpha*pk-beta*pkm1)/(2*k*(k+a+b)*(2*k+a+b-2))
        pols(k+1)=pkp1
 2000 continue
c
        return
        end
c
c
c
c
c
        subroutine kjacopols2(x,a,b,n,pols,ders)
        implicit real *8 (a-h,o-z)
        dimension pols(1), ders(1)
c
c       evaluates a bunch of Jacobi polynomials (together
c       with their derivatives) at the user-provided point x
c
c       ... if n=0 or n=1 - exit
c
        if(n .ge. 2) goto 1200
        pols(1)=1
        ders(1)=0
        if(n .eq. 0) return
c
        pols(2)=(a/2-b/2)+(1+a/2+b/2)*x
        ders(2)=(1+a/2+b/2)
        return
 1200 continue
c
        pols(1)=1
        ders(1)=0
        pols(2)=(a/2-b/2)+(1+a/2+b/2)*x
        ders(2)=(1+a/2+b/2)
c
c       n is greater than 2. conduct recursion
c
        pkm1=1
        dkm1=0
        pk=(a/2-b/2)+(1+a/2+b/2)*x
        dk=(1+a/2+b/2)
c
        pk=1
        dk=0
        pkp1=(a/2-b/2)+(1+a/2+b/2)*x
        dkp1=(1+a/2+b/2)
c
        do 2000 k=2,n
c
        pkm1=pk
        pk=pkp1
        dkm1=dk
        dk=dkp1
        alpha1=(2*k+a+b-1)*(a**2-b**2)
        alpha2=(2*k+a+b-1)*((2*k+a+b-2)*(2*k+a+b))
        beta=2*(k+a-1)*(k+b-1)*(2*k+a+b)
        gamma=(2*k*(k+a+b)*(2*k+a+b-2))
        pkp1=((alpha1+alpha2*x)*pk-beta*pkm1)/gamma
        dkp1=((alpha1+alpha2*x)*dk-beta*dkm1+alpha2*pk)/gamma
        pols(k+1)=pkp1
        ders(k+1)=dkp1
 2000 continue
c
        return
        end

max_stars_repo_name: zgimbutas/triasymq

max_issues_repo_name: zgimbutas/triasymq

max_forks_repo_name: zgimbutas/triasymq


text:
      Subroutine droplet3(pix,ip,matr,dar,thr,nx,ny,lp)
      integer nx, ny, lp,ix, iy, nxx, nyy,st,ip,jj
      integer matro(0:nx-1,0:ny-1),cxj,cyj,iw
      integer j,cx(0:lp-1),cy(0:lp-1),jjj
      real fs,thr
      real matr(0:nx-1,0:ny-1),dar(0:nx-1,0:ny-1)
      real pix(0:lp-1)
Cf2py intent(out) pix
Cf2py intent(out) ip
Cf2py intent(in,c) matr
Cf2py intent(in,c) dar
Cf2py intent(in) thr
Cf2py intent(in) lp
Cf2py depend(nx,ny) matr
Cf2py depend(nx,ny) dar
Cf2py depend(lp) pix

C     this is to find droplets really good program!!!
      j=0
      ip=0
      nxx=nx-1
      nyy=ny-1
      fs=0
      st=-1
      iw=-1
      do 42 iy=0,nyy,1
       do 41 ix=0,nxx,1
        matr(ix,iy)=matr(ix,iy)-dar(ix,iy)
        if (matr(ix,iy) .LT. thr) then
         matr(ix,iy)=0.0
        endif
        matro(ix,iy)=0
41     continue
42    continue
      do 40 iy=0,nyy,1
       do 30 ix=0,nxx,1
        if (matr(ix,iy) .GT. 0.0 .AND. matro(ix,iy) .EQ. 0) then
         matro(ix,iy)=1
         st=0
         cx(j)=ix
         cy(j)=iy
         iw=-1
        endif
        do while (st .GE. 0)
         iw=iw+1
         cxj=cx(iw)
         cyj=cy(iw)
         st=st-1
         if (cxj-1 .GE. 0) then
          if (matr(cxj-1,cyj) .GT. 0.0 .AND. matro(cxj-1,cyj)
     1       .EQ. 0) then
             matro(cxj-1,cyj)=1
             j=j+1
            cx(j)=cxj-1
            cy(j)=cyj
            st=st+1
          endif
         endif
         if (cxj+1 .LE. nxx) then
          if (matr(cxj+1,cyj) .GT. 0.0 .AND. matro(cxj+1,cyj)
     1       .EQ. 0) then
             matro(cxj+1,cyj)=1
             j=j+1
             cx(j)=cxj+1
             cy(j)=cyj
             st=st+1
          endif
         endif
         if (cyj-1 .GE. 0) then
          if (matr(cxj,cyj-1) .GT. 0.0 .AND. matro(cxj,cyj-1)
     1       .EQ. 0) then
             matro(cxj,cyj-1)=1
             j=j+1
             cx(j)=cxj
             cy(j)=cyj-1
             st=st+1
          endif
         endif
         if (cyj+1 .LE. nyy) then
          if (matr(cxj,cyj+1) .GT. 0.0 .AND. matro(cxj,cyj+1)
     1       .EQ. 0) then
             matro(cxj,cyj+1)=1
             j=j+1
             cx(j)=cxj
             cy(j)=cyj+1
             st=st+1
          endif
         endif
        enddo
        if (iw .GE. 0) then
         do 20 jjj=0,iw,1
          fs=fs+matr(cx(jjj),cy(jjj))
          cx(jjj)=0
          cy(jjj)=0
20       continue
         pix(ip)=fs
         ip=ip+1
        endif
       fs=0
       j=0
       jj=0
       st=-1
       iw=-1
30     continue
40    continue
      return
      end Subroutine

max_stars_repo_name: kif/dynamix

max_issues_repo_name: kif/dynamix

max_forks_repo_name: kif/dynamix


text:
!
!***********************************************************************
!                                                                      *
      PROGRAM jj2lsj2K
!                                                                      *
!     This MAIN program controls the transformation of atomic states,  *
!     which are given in a jj-coupled CSF basis, into an LS-coupled    *
!     basis. The program requires a jj-coupled basis in standard order *
!     where, if both subshells | n j = l-1/2> and | n j = l+1/2>       *
!     of a given shell (nl) occurs, they always follow successively    *
!     in this order. The LS-coupled basis, moreover, is given in       *
!     the same sequence of shells.                                     *
!                                                                      *
!     All LS-jj transformation coefficients are precalculated and      *
!     'stored' in the modules rabs_lsj_data_1, rabs_lsj_data_2 and     *
!     rabs_lsj_data_3.                                                 *
!                                                                      *
!     Calls:  FACTT, SETISO, JJ2LSJ, starttime, stoptime.              *
!                                                                      *
!     Written by G. Gaigalas,                                          *
!     NIST                                                  May 2011   *
!     VILNIUS                                               May 2017   *
!                                                                      *
!     Modified by G. Gaigalas and C. Cychen                     2021   *
!                                                                      *
!***********************************************************************
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE jj2lsj_code
      IMPLICIT NONE
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer           :: ncount1
!-----------------------------------------------
      call starttime (ncount1, 'jj2lsj')
      print *, " "
      print *, "jj2lsj: Transformation of ASFs from a jj-coupled CSF basis"
      print *, "        into an LS-coupled CSF basis  (Fortran 95 version)"
      print *, "        (C) Copyright by   G. Gaigalas and Ch. F. Fischer,"
      print *, "        (2021)."
      print *, "        Input files: name.c, name.(c)m"
      print *, "        Ouput files: name.lsj.lbl,"
      print *, "         (optional)  name.lsj.c, name.lsj.j,"
      print *, "                     name.uni.lsj.lbl, name.uni.lsj.sum"
      print *, " "
!
!  Set up the table of logarithms of factorials
      call factt
!
      CALL setiso('isodata')
      CALL jj2lsj
!
      call stoptime (ncount1, 'jj2lsj')
      stop
      end program jj2lsj2K

max_stars_repo_name: compas/grasp

max_issues_repo_name: compas/grasp

max_forks_repo_name: compas/grasp


!        Generated by TAPENADE     (INRIA, Ecuador team)
!  tapenade 3.x
!
!  Differentiation of lamm in forward (tangent) mode:
!   variations   of useful results: lamm
!   with respect to varying inputs: tlocal lammref
! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.
!> @brief calculate temperature dependent thermal conductivity
!> @param[in] lammref thermal conductivity from input file
!> @param[in] tlocal temperature
!> @param[in] tref reference temperature
!> @param[in] ismpl local sample index
!> @return temperature dependent thermal conductivity
!> @details
!> calculate temperature dependent thermal conductivity\n
!>
!>        lam_zoth/haenel = (770degC/(350degC+T) + 0.7) W/mK
!>
!> If T > 800degC, use the formula from zoth & haenel, 1988. \n
!>
!>        lamm = lam_zoth/haenel
!>
!> If T < 800degC, use the same formula with a factor `fct`
!>
!>        lamm = fct * lam_zoth/haenel
!>
!> The factor `fct` introduces an additional temperature dependence
!> such that\n
!> 1. `lamm(20degC) = lammref`\n
!> 2. `lamm(800degC) = lam_zoth/haenel(800degC)\n\n
!>
!> Thus, the thermal conductivity in the input file should resemble
!> information about the value of the thermal conductivity of the
!> matrix at temperature 20 degC.\n\n
!>
!> Sources: \n
!> Zoth, G., & Haenel, R, Handbook of terrestrial heat-flow density
!> determination, (1988), Appendix 10.1 Thermal
!> Conductivity. http://dx.doi.org/10.1007/978-94-009-2847-3\n\n
!>
!> Lehmann, H., Wang, K., & Clauser, C., Parameter identification and
!> uncertainty analysis for heat transfer at the ktb drill site using
!> a 2-d inverse method, Tectonophysics, 291(1-4), 179–194
!> (1998). Section 2.2 http://dx.doi.org/10.1016/s0040-1951(98)00039-0
!> \n
DOUBLE PRECISION FUNCTION g_LAMM(lammref, g_lammref, tlocal, &
& g_tlocal, tref, ismpl, lamm)
  IMPLICIT NONE
! Sample index
  INTEGER :: ismpl
! Thermal conductivity of matrix at tref=20degC
  DOUBLE PRECISION, INTENT(IN) :: lammref
  DOUBLE PRECISION, INTENT(IN) :: g_lammref
! Local temperature [degC]
  DOUBLE PRECISION, INTENT(IN) :: tlocal
  DOUBLE PRECISION, INTENT(IN) :: g_tlocal
! Reference temperature [20 degC]
  DOUBLE PRECISION, INTENT(IN) :: tref
! Upper limit temperature, where the approximation becomes
! equal to the general Zoth/Haenel formula
  DOUBLE PRECISION, PARAMETER :: tlimit=800.0d0
! lam_zoth/haenel at tlocal
  DOUBLE PRECISION :: lamm_zh
  DOUBLE PRECISION :: g_lamm_zh
! lam_zoth/haenel at tref
  DOUBLE PRECISION :: lamm_zhref
! Reference Interpolation-factor
  DOUBLE PRECISION :: fctref
  DOUBLE PRECISION :: g_fctref
! Weight: Quotient of temperature differences
  DOUBLE PRECISION :: twgt
  DOUBLE PRECISION :: g_twgt
! Interpolation factor
  DOUBLE PRECISION :: fct
  DOUBLE PRECISION :: g_fct
  DOUBLE PRECISION :: lamm
  IF (tlocal .GT. tlimit) THEN
! lam_zoth/haenel at tlocal
    g_lamm = -(770.0d0*g_tlocal/(tlocal+350.0d0)**2)
    lamm = 770.0d0/(350.0d0+tlocal) + 0.7d0
  ELSE
! lam_zoth/haenel at tlocal
    g_lamm_zh = -(770.0d0*g_tlocal/(tlocal+350.0d0)**2)
    lamm_zh = 770.0d0/(350.0d0+tlocal) + 0.7d0
! lam_zoth/haenel at tref
    lamm_zhref = 770.0d0/(350.0d0+tref) + 0.7d0
! Reference Interpolation-factor: Input lamm at tref divided
! by lam_zoth/haenel at tref
    g_fctref = g_lammref/lamm_zhref
    fctref = lammref/lamm_zhref
! Quotient of temperature differences, local minus reference
! divided by limit minus reference
    g_twgt = g_tlocal/(tlimit-tref)
    twgt = (tlocal-tref)/(tlimit-tref)
! Interpolation factor between fctref at tref and 1 at tlimit
    g_fct = (1-twgt)*g_fctref - (fctref-1.0)*g_twgt
    fct = fctref*(1-twgt) + twgt
! Final lam: input at tref and lam_zoth/haenel at tlimit
    g_lamm = lamm_zh*g_fct + fct*g_lamm_zh
    lamm = fct*lamm_zh
  END IF
  RETURN
END FUNCTION g_LAMM


    program demo_maxval
    implicit none
    integer,save :: ints(3,5)= reshape([&
       1,  2,  3,  4,  5, &
      10, 20, 30, 40, 50, &
      11, 22, 33, 44, 55  &
    ],shape(ints),order=[2,1])

       write(*,*) maxval(ints)
       write(*,*) maxval(ints,dim=1)
       write(*,*) maxval(ints,dim=2)
       ! find biggest number less than 30 with mask
       write(*,*) maxval(ints,mask=ints.lt.30)
    end program demo_maxval


text:
Address(Michelangelo Place) is a residential Culdesacs culdesac in the Wildhorse section of East Davis.

 Intersecting Streets

  Leonardo Street


text:
Aaron is a mathematics graduate student.  He has a http://www.math.ucdavis.edu/~arosenberg website, and mailto:arosenberg@ucdavis.edu an email address.

Aaron adopted a pair of 5monthold kittens from the Orphan Kitten Project in August.  Asher (happy in Hebrew) is nearly twice as big as his sister, Tzviya (gazelle in Hebrew).  Tzviya loves chin rubs, chairs, kneading, watching Aaron cook, and watching anything moving on the computer screen.  Tzviya is slowly learning how to play fetch.  Asher loves belly rubs, playing soccer with any objects that will roll or slide, food, bathtubs, and sinks.  Asher is slowly learning that he cant get a belly rub while Aaron is sleeping...therefore he must wake Aaron up to get his 3am belly rub.  ...and 4am.  ...and 5am.  ...and 6am.  ...and 7am.  ...and 8am.

Tzviya and Asher have learned new tricks!  Tzviya is getting the hang of XfetchX play with ball and dont return it and sit.  Asher is becoming increasingly talented at chewing through cords of all shapes and forms from USB and speaker cables to highvoltage, highcurrent surgeprotector cables.  Apparently, kittens have 50% resistance to lightning damage, in addition to the XsixX four lives the mischievous little brat has remaining.

I also now have three baby corn snakes.  They are albino (aka, pink!) and heterozygous for blood corn, so they should have some blood corn babies when they grow up.  Pitriya doesnt like eating.
Vidra (Badger in Hungarian)
Pitriya (Mushroom in Hebrew)
Zmiya (Snake in Russian)

  Mathematical interests:
   topology
    knot theory
     Khovanov homology
      cobordisms: thingies that thing between things
       further abstract abstracted abstraction, as long as its pretty.
    Legendrian knot theory
  Nonmathematical interests:
   stuff
   shiny stuff
   tasty stuff
   Dune
   Truth and Beauty (not the quark varieties!)
   your mom

  Favorite quote written in a paper or problem set so far: For w in S, note that w is never happy unless w has cookies.
  I am firmly convinced that all mathematical terms should be words like happy, perfectly normal, and quandle.
  Cobordism?!  Nah: thingy that things between things.

20051207 13:33:37 nbsp i know you! Users/MichelleAccurso
  I know me too!  er.  wait. nevermind, no I dont. Users/AaronRosenberg AR
   Thats sad. Ill help you with that sometime... Users/MichelleAccurso

20060804 16:51:01 nbsp i like your fractal art work.  for some reason i get an error when i try to access your senior thesis. Users/MattHh

  Im http://www.math.ucdavis.edu/~arosenberg/papers/hcthesissemifinal.pdf still seeing an error.  Yeah, I want to read your thesis. :)  Users/PhilipNeustrom
  Ill see if I can figure out whats going on sometime in the next week or two.  I swear I fixed it and saw it working on multiple computers a couple months ago... Users/AaronRosenberg AR

20070816 18:35:20 nbsp Its good to know that in a world that is corrupted by journalists, the Oxford comma still lives on! Users/WilliamLewis

20071112 16:19:50 nbsp Just curious what methods do you use to generate your fractal art. Is it paint, computers,...? Users/MattHh

   It says at the bottom of http://www.math.ucdavis.edu/~arosenberg/fractal/fractal.html this page the computer programs hes used. Users/PhilipNeustrom

20071130 09:09:26 nbsp Your kitties are too cute!  Any advice on how you taught them to fetch? Users/FeralFeline

20071209 16:10:19 nbsp hello. do you know about the Wiki BBQ Jan 2008? Everyones invited. Were trying to spread the word, so if you want to let other people on the wiki know, go ahead and copy and past this. gracias!  Users/JessicaRockwell

20071210 07:56:32 nbsp Aaron  That sounds like a cool book on cats!  Boo and Taffeta are great cats and it is a shame that they havent gotten adopted yet.  Sometimes being at Petsmart in a cage with all the unfamiliar sounds and actions can be scary for the cats and their true greatness does not shine through. Users/FeralFeline