text:
C
C $Id: srfgk.f,v 1.5 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 SRFGK (X,Y,Z,M,MX,NX,NY,S,STEREO)
C
      DIMENSION       X(NX)      ,Y(NY)      ,Z(MX,NY)   ,M(2,NX,NY) ,
     1                S(6)
      DIMENSION       MXS(2)     ,MXF(2)     ,MXJ(2)     ,MYS(2)     ,
     1                MYF(2)     ,MYJ(2)
      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
      COMMON /PWRZ1S/ XXMIN      ,XXMAX      ,YYMIN      ,YYMAX      ,
     1                ZZMIN      ,ZZMAX      ,DELCRT     ,EYEX       ,
     2                EYEY       ,EYEZ
      COMMON /SRFIP1/ IFR        ,ISTP       ,IROTS      ,IDRX       ,
     1                IDRY       ,IDRZ       ,IUPPER     ,ISKIRT     ,
     2                NCLA       ,THETA      ,HSKIRT     ,CHI        ,
     3                CLO        ,CINC       ,ISPVAL
      DATA JF,IF,LY,LX,ICNST/1,1,2,2,0/
      CALL Q8QST4 ('GRAPHX','SRFACE','SRFGK','VERSION 01')
      BIGEST = R1MACH(2)
      MMXX = MX
      NNXX = NX
      NNYY = NY
      STER = STEREO
      NXP1 = NNXX+1
      NYP1 = NNYY+1
      NLA = NCLA
      NSPVAL = ISPVAL
      NDRZ = IDRZ
      IF (IDRZ .NE. 0)
     1    CALL CLSET (Z,MMXX,NNXX,NNYY,CHI,CLO,CINC,NLA,40,CL,NCL,
     2                ICNST,IOFFP,SPVAL,BIGEST)
      IF (IDRZ .NE. 0) NDRZ = 1-ICNST
      STHETA = SIN(STER*THETA)
      CTHETA = COS(STER*THETA)
      RX = S(1)-S(4)
      RY = S(2)-S(5)
      RZ = S(3)-S(6)
      D1 = SQRT(RX*RX+RY*RY+RZ*RZ)
      D2 = SQRT(RX*RX+RY*RY)
      DX = 0.
      DY = 0.
      IF (STEREO .EQ. 0.) GO TO  20
      D1 = D1*STEREO*THETA
      IF (D2 .GT. 0.) GO TO  10
      DX = D1
      GO TO  20
   10 AGL = ATAN2(RX,-RY)
      DX = D1*COS(AGL)
      DY = D1*SIN(AGL)
   20 IROT = IROTS
      NPIC = 1
      IF (STER .NE. 0.) NPIC = 2
      FACT = 1.
      IF (NRSWT .NE. 0) FACT = RZERO/D1
      IF (ISTP.EQ.0 .AND. STER.NE.0.) IROT = 1
      DO 570 IPIC=1,NPIC
         NUPPER = IUPPER
         IF (IFR .LT. 0) CALL FRAME
C
C SET UP MAPPING FROM FLOATING POINT 3-SPACE TO CRT SPACE.
C
         SIGN1 = IPIC*2-3
         EYEX = S(1)+SIGN1*DX
         POIX = S(4)+SIGN1*DX
         EYEY = S(2)+SIGN1*DY
         POIY = S(5)+SIGN1*DY
         EYEZ = S(3)
         POIZ = S(6)
         LL = 0
         XEYE = EYEX
         YEYE = EYEY
         ZEYE = EYEZ
         CALL TRN32S (POIX,POIY,POIZ,XEYE,YEYE,ZEYE,0)
         LL = IPIC+2*ISTP+3
         IF (STER .EQ. 0.) LL = 1
         IF (NRSWT .NE. 0) GO TO 100
         XXMIN = X(1)
         XXMAX = X(NNXX)
         YYMIN = Y(1)
         YYMAX = Y(NNYY)
         UMIN = BIGEST
         VMIN = BIGEST
         ZZMIN = BIGEST
         UMAX = -UMIN
         VMAX = -VMIN
         ZZMAX = -ZZMIN
         DO  40 J=1,NNYY
            DO  30 I=1,NNXX
               ZZ = Z(I,J)
               IF (IOFFP.EQ.1 .AND. ZZ.EQ.SPVAL) GO TO  30
               ZZMAX = MAX(ZZMAX,ZZ)
               ZZMIN = MIN(ZZMIN,ZZ)
               CALL TRN32S (X(I),Y(J),Z(I,J),UT,VT,DUMMY,1)
               UMAX = MAX(UMAX,UT)
               UMIN = MIN(UMIN,UT)
               VMAX = MAX(VMAX,VT)
               VMIN = MIN(VMIN,VT)
   30       CONTINUE
   40    CONTINUE
         IF (ISKIRT .NE. 1) GO TO  70
         NXSTP = NNXX-1
         NYSTP = NNYY-1
         DO  60 J=1,NNYY,NYSTP
            DO  50 I=1,NNXX,NXSTP
               CALL TRN32S (X(I),Y(J),HSKIRT,UT,VT,DUMMY,1)
               UMAX = MAX(UMAX,UT)
               UMIN = MIN(UMIN,UT)
               VMAX = MAX(VMAX,VT)
               VMIN = MIN(VMIN,VT)
   50       CONTINUE
   60    CONTINUE
   70    CONTINUE
         WIDTH = UMAX-UMIN
         HIGHT = VMAX-VMIN
         DIF = .5*(WIDTH-HIGHT)
         IF (DIF)  80,100, 90
   80    UMIN = UMIN+DIF
         UMAX = UMAX-DIF
         GO TO 100
   90    VMIN = VMIN-DIF
         VMAX = VMAX+DIF
  100    XEYE = EYEX
         YEYE = EYEY
         ZEYE = EYEZ
         CALL TRN32S (POIX,POIY,POIZ,XEYE,YEYE,ZEYE,0)
         DO 120 J=1,NNYY
            DO 110 I=1,NNXX
               CALL TRN32S (X(I),Y(J),Z(I,J),UT,VT,DUMMY,1)
               M(1,I,J) = UT
               M(2,I,J) = VT
  110       CONTINUE
  120    CONTINUE
C
C INITIALIZE UPPER AND LOWER VISIBILITY ARRAYS
C
         DO 130 K=1,1024
            LIMU(K) = 0
            LIML(K) = 1024
  130    CONTINUE
C
C FIND ORDER TO DRAW LINES
C
         NXPASS = 1
         IF (S(1) .GE. X(NNXX)) GO TO 160
         IF (S(1) .LE. X(1)) GO TO 170
         DO 140 I=2,NNXX
            LX = I
            IF (S(1) .LE. X(I)) GO TO 150
  140    CONTINUE
  150    MXS(1) = LX-1
         MXJ(1) = -1
         MXF(1) = 1
         MXS(2) = LX
         MXJ(2) = 1
         MXF(2) = NNXX
         NXPASS = 2
         GO TO 180
  160    MXS(1) = NNXX
         MXJ(1) = -1
         MXF(1) = 1
         GO TO 180
  170    MXS(1) = 1
         MXJ(1) = 1
         MXF(1) = NNXX
  180    NYPASS = 1
         IF (S(2) .GE. Y(NNYY)) GO TO 210
         IF (S(2) .LE. Y(1)) GO TO 220
         DO 190 J=2,NNYY
            LY = J
            IF (S(2) .LE. Y(J)) GO TO 200
  190    CONTINUE
  200    MYS(1) = LY-1
         MYJ(1) = -1
         MYF(1) = 1
         MYS(2) = LY
         MYJ(2) = 1
         MYF(2) = NNYY
         NYPASS = 2
         GO TO 230
  210    MYS(1) = NNYY
         MYJ(1) = -1
         MYF(1) = 1
         GO TO 230
  220    MYS(1) = 1
         MYJ(1) = 1
         MYF(1) = NNYY
C
C PUT ON SKIRT ON FRONT SIDE IF WANTED
C
  230    IF (NXPASS.EQ.2 .AND. NYPASS.EQ.2) GO TO 490
         IF (ISKIRT .EQ. 0) GO TO 290
         IN = MXS(1)
         IF = MXF(1)
         JN = MYS(1)
         JF = MYF(1)
         IF (NYPASS .NE. 1) GO TO 260
         CALL TRN32S (X(1),Y(JN),HSKIRT,UX1,VX1,DUMMY,1)
         CALL TRN32S (X(NNXX),Y(JN),HSKIRT,UX2,VX2,DUMMY,1)
         QU = (UX2-UX1)/(X(NNXX)-X(1))
         QV = (VX2-VX1)/(X(NNXX)-X(1))
         YNOW = Y(JN)
         DO 240 I=1,NNXX
            CALL TRN32S (X(I),YNOW,HSKIRT,RU,RV,DUMMY,1)
            CALL DRAWS (INT(RU),INT(RV),M(1,I,JN),M(2,I,JN),1,1)
  240    CONTINUE
         CALL DRAWS (INT(UX1),INT(VX1),INT(UX2),INT(VX2),1,1)
         IF (IDRY .NE. 0) GO TO 260
         DO 250 I=2,NNXX
            CALL DRAWS (M(1,I-1,JN),M(2,I-1,JN),M(1,I,JN),M(2,I,JN),1,1)
  250    CONTINUE
  260    IF (NXPASS .NE. 1) GO TO 290
         CALL TRN32S (X(IN),Y(1),HSKIRT,UY1,VY1,DUMMY,1)
         CALL TRN32S (X(IN),Y(NNYY),HSKIRT,UY2,VY2,DUMMY,1)
         QU = (UY2-UY1)/(Y(NNYY)-Y(1))
         QV = (VY2-VY1)/(Y(NNYY)-Y(1))
         XNOW = X(IN)
         DO 270 J=1,NNYY
            CALL TRN32S (XNOW,Y(J),HSKIRT,RU,RV,DUMMY,1)
            CALL DRAWS (INT(RU),INT(RV),M(1,IN,J),M(2,IN,J),1,1)
  270    CONTINUE
         CALL DRAWS (INT(UY1),INT(VY1),INT(UY2),INT(VY2),1,1)
         IF (IDRX .NE. 0) GO TO 290
         DO 280 J=2,NNYY
            CALL DRAWS (M(1,IN,J-1),M(2,IN,J-1),M(1,IN,J),M(2,IN,J),1,1)
  280    CONTINUE
C
C PICK PROPER ALGORITHM
C
  290    LI = MXJ(1)
         MI = MXS(1)-LI
         NI = ABS(MI-MXF(1))
         LJ = MYJ(1)
         MJ = MYS(1)-LJ
         NJ = ABS(MJ-MYF(1))
C
C IF NXPASS IS 1, PUT THE I LOOP OUTERMOST.  OTHERWISE, PUT THE J LOOP
C OUTERMOST.
C
         IF (NXPASS.EQ.2) GO TO 360
         IF (ISKIRT.NE.0 .OR. NYPASS.NE.1) GO TO 310
         I = MXS(1)
         DO 300 J=2,NNYY
            CALL DRAWS (M(1,I,J-1),M(2,I,J-1),M(1,I,J),M(2,I,J),0,1)
  300    CONTINUE
  310    DO 350 II=1,NNXX
            I = MI+II*LI
            IPLI = I+LI
            IF (NYPASS .EQ. 1) GO TO 320
            K = MYS(1)
            L = MYS(2)
            IF (IDRX .NE. 0)
     1          CALL DRAWS (M(1,I,K),M(2,I,K),M(1,I,L),M(2,I,L),1,1)
            IF (NDRZ.NE.0 .AND. II.NE.NI)
     1          CALL CTCELL (Z,MMXX,NNXX,NNYY,M,MIN(I,I+LI),K)
  320       DO 340 JPASS=1,NYPASS
               LJ = MYJ(JPASS)
               MJ = MYS(JPASS)-LJ
               NJ = ABS(MJ-MYF(JPASS))
               DO 330 JJ=1,NJ
                  J = MJ+JJ*LJ
                  JPLJ = J+LJ
                  IF (IDRX.NE.0 .AND. JJ.NE.NJ)
     1                CALL DRAWS (M(1,I,J),M(2,I,J),M(1,I,JPLJ),
     2                            M(2,I,JPLJ),1,1)
                  IF (I.NE.MXF(1) .AND. IDRY.NE.0)
     1                CALL DRAWS (M(1,IPLI,J),M(2,IPLI,J),M(1,I,J),
     2                            M(2,I,J),1,1)
                  IF (NDRZ.NE.0 .AND. JJ.NE.NJ .AND. II.NE.NNXX)
     1                CALL CTCELL (Z,MMXX,NNXX,NNYY,M,MIN(I,I+LI),
     2                             MIN(J,J+LJ))
  330          CONTINUE
  340       CONTINUE
  350    CONTINUE
         GO TO 430
  360    IF (ISKIRT.NE.0 .OR. NXPASS.NE.1) GO TO 380
         J = MYS(1)
         DO 370 I=2,NNXX
            CALL DRAWS (M(1,I-1,J),M(2,I-1,J),M(1,I,J),M(2,I,J),0,1)
  370    CONTINUE
  380    DO 420 JJ=1,NNYY
            J = MJ+JJ*LJ
            JPLJ = J+LJ
            IF (NXPASS .EQ. 1) GO TO 390
            K = MXS(1)
            L = MXS(2)
            IF (IDRY .NE. 0)
     1          CALL DRAWS (M(1,K,J),M(2,K,J),M(1,L,J),M(2,L,J),1,1)
            IF (NDRZ.NE.0 .AND. JJ.NE.NJ)
     1          CALL CTCELL (Z,MMXX,NNXX,NNYY,M,K,MIN(J,J+LJ))
  390       DO 410 IPASS=1,NXPASS
               LI = MXJ(IPASS)
               MI = MXS(IPASS)-LI
               NI = ABS(MI-MXF(IPASS))
               DO 400 II=1,NI
                  I = MI+II*LI
                  IPLI = I+LI
                  IF (IDRY.NE.0 .AND. II.NE.NI)
     1                CALL DRAWS (M(1,I,J),M(2,I,J),M(1,IPLI,J),
     2                            M(2,IPLI,J),1,1)
                  IF (J.NE.MYF(1) .AND. IDRX.NE.0)
     1                CALL DRAWS (M(1,I,JPLJ),M(2,I,JPLJ),M(1,I,J),
     2                            M(2,I,J),1,1)
                  IF (NDRZ.NE.0 .AND. II.NE.NI .AND. JJ.NE.NNYY)
     1                CALL CTCELL (Z,MMXX,NNXX,NNYY,M,MIN(I,I+LI),
     2                             MIN(J,J+LJ))
  400          CONTINUE
  410       CONTINUE
  420    CONTINUE
  430    IF (ISKIRT .EQ. 0) GO TO 520
C
C FIX UP IF SKIRT IS USED WITH LINES ONE WAY.
C
         IF (IDRX .NE. 0) GO TO 460
         DO 450 IPASS=1,NXPASS
            IF (NXPASS .EQ. 2) IF = 1+(IPASS-1)*(NNXX-1)
            DO 440 J=2,NNYY
               CALL DRAWS (M(1,IF,J-1),M(2,IF,J-1),M(1,IF,J),M(2,IF,J),
     1                     1,0)
  440       CONTINUE
  450    CONTINUE
  460    IF (IDRY .NE. 0) GO TO 520
         DO 480 JPASS=1,NYPASS
            IF (NYPASS .EQ. 2) JF = 1+(JPASS-1)*(NNYY-1)
            DO 470 I=2,NNXX
               CALL DRAWS (M(1,I-1,JF),M(2,I-1,JF),M(1,I,JF),M(2,I,JF),
     1                     1,0)
  470       CONTINUE
  480    CONTINUE
         GO TO 520
C
C ALL VISIBLE IF VIEWED FROM DIRECTLY ABOVE OR BELOW.
C
  490    IF (NUPPER.GT.0 .AND. S(3).LT.S(6)) GO TO 520
         IF (NUPPER.LT.0 .AND. S(3).GT.S(6)) GO TO 520
         NUPPER = 1
         IF (S(3) .LT. S(6)) NUPPER = -1
         DO 510 I=1,NNXX
            DO 500 J=1,NNYY
               IF (IDRX.NE.0 .AND. J.NE.NNYY)
     1             CALL DRAWS (M(1,I,J),M(2,I,J),M(1,I,J+1),M(2,I,J+1),
     2                         1,0)
               IF (IDRY.NE.0 .AND. I.NE.NNXX)
     1             CALL DRAWS (M(1,I,J),M(2,I,J),M(1,I+1,J),M(2,I+1,J),
     2                         1,0)
               IF (IDRZ.NE.0 .AND. I.NE.NNXX .AND. J.NE.NNYY)
     1             CALL CTCELL (Z,MMXX,NNXX,NNYY,M,I,J)
  500       CONTINUE
  510    CONTINUE
  520    IF (STER .EQ. 0.) GO TO 560
         IF (ISTP) 540,530,550
  530    CALL FRAME
  540    CALL FRAME
         GO TO 570
  550    IF (IPIC .NE. 2) GO TO 570
  560    IF (IFR .GT. 0) CALL FRAME
  570 CONTINUE
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
C example for fsw_classic heating source
C fixed format for fortran77 (programming area 7-72)

C DFLUX for friction stir welding
C Systeme International
C welding speed 300 mm/min
C rotational velocity 600 rpm
C cylinder pin without extra topology

C ---------------------------------------
      SUBROUTINE DFLUX(FLUX,SOL,KSTEP,KINC,TIME,NOEL,NPT,COORDS,
     & JLTYP,TEMP,PRESS,SNAME)
      INCLUDE 'ABA_PARAM.INC'
      DIMENSION FLUX(2), TIME(2), COORDS(3)
      CHARACTER*80 SNAME
      real x,y,z ! coordinates
      real x0,y0,z0 ! initial position
      real x1,y1,z1 ! final position
      real t0,t1 ! dwelling time
      real v0_1,omega0_1 ! welding parameters 0_1
      real r_shoulder,r_pinr,l_pin ! shape parameters
      real P,miu,eff_heat,k ! internal parameters
      real x_c,y_c,z_c ! current position of welding tool
      real r_c ! radius between intergal point and tool
      real force ! frictional force
      real weld_t0_1 ! final time for weld 0_1
      REAL PI; PARAMETER ( PI= 3.14159 )

C ----------------------------------------
C definition of coordinates
      x=COORDS(1) ! coordinate of X
      y=COORDS(2) ! coordinate of Y
      z=COORDS(3) ! coordinate of Z
C initial position of welding tool
      x0=0.
      y0=0.
      z0=0.
C final position of welding tool
      x1=0.025
      y1=0.
      z1=0.
C dwelling time
      dwell_t0=3. ! initial
      dwell_t1=3. ! final
C welding speed and rotational velocity
      v0_1=0.005 ! welding speed
      omega0_1=600.*2.*PI/60. ! rotational velocity
C shape of welding tool
      r_shoulder=0.007 ! radius of shoulder
      r_pinr=0.002 ! radius of pin root
      l_pin=0.003 ! length of pin
C internal parameters
      P=15286242.0 ! axial pressure of welding tool
      miu=0.37 ! coefficient of friction
      eff_heat=0.98 ! coefficient of heat generation
      k=0.33 ! heating ratio of pin to shoulder
C frictional force
      if (SOL.lt.390) then
        force=miu*P
      else if (SOL.lt.450) then
        force=0.577*(98.1-(98.1-47)*(SOL-390)/(450-390))*1e6
      else if (SOL.lt.550) then
        force=0.577*(47-(47-5)*(SOL-450)/(550-450))*1e6
      else
        force=0.577*1e6
      end if

C ----------------------------------------
C initialize the flux value
      FLUX(1)=0.
C real-time position of welding tool
      weld_t0_1=dwell_t0+(sqrt((x1-x0)**2+(y1-y0)**2+
     & (z1-z0)**2)/v0_1)
      if (TIME(2).le.dwell_t0) then
        x_c=x0
        y_c=y0
        z_c=z0
      else if (TIME(2).gt.weld_t0_1) then
        x_c=x1
        y_c=y0
        z_c=z0
      else
        x_c=x0+(v0_1*(TIME(2)-dwell_t0))
        y_c=y0
        z_c=z0
      end if
C distance between tool and integral point
      r_c=sqrt((x-x_c)**2+(y-y_c)**2+(z-z_c)**2)
C fsw_classic heating source
      if ((JLTYP.eq.0).and.(r_c.le.r_shoulder)) then
      	FLUX(1)=force*r_c*omega0_1*eff_heat
      else if ((JLTYP.EQ.1)
     &            .and.(r_c.le.r_pinr)
     &            .and.(z.ge.(z0-l_pin))) then
      	FLUX(1)=(2.*k*force*omega0_1*eff_heat*
     &            (r_shoulder**3-r_pinr**3))/
     &            (3.*l_pin*r_pinr**2)
      else
      	FLUX(1)=0.
      end if

      return
      END

ext: for


text:
wiki:WikiPedia:Lake_Tahoe Lake Tahoe can be reached by getting on the Interstate 80 80 and either staying on it or switching to the Highway 50 50. If you stay on 80 you can take 267 (North Shore Boulevard) or 89 (West River Road) to the north shore. 50 takes you to wiki:southlaketahoe South Lake Tahoe.

UC Davis maintains a research station that houses the Tahoe Research Group to study the lake.  In July 2010, UC Davis opened its Eriksson Education Center in Tahoe City.  The location has been used as a research center since 1975, but the main building and surrounding area have undergone substantial renovation since then.

Tahoe Lodging
Plenty of places to stay around Tahoe

     Reno(lots of casinos)
     South Lake Tahoe(some casinos)
     Truckee Tahoe City
     Vaction rental around the lake

Reno is the cheapest stay in Tahoe region.  South Lake Tahoe is the 2nd cheapest stay.

If you want to stay in a cozy place, find a vacation rental at http://www.craigslist.com Craigslist, http://www.vrbo.com VRBO.com , or http://www.tahoevacationrentalsite.com TahoeVacationRentalSite.com

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Subroutine BuildNumber(nbuild)

Logical res
Character(255) FilObj,FilInc
!Character(36) str

Include 'BuildNumber.fi'

! BuildNumber saves the build number of the current program version
! using the BuilNumber.fi file
!
! Usage in the program (no additional code is needed):
!
!		.....
!		Call Build Number(nBuild)
!		Write(*,'('' Build number: '',i5')nBuild
!		....
!
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Before use, adjust the directory names where two your signal files are located!
FilObj='d:\Projects\ReadExcitedStates\ReadExcitedStates\Debug\BuildNumber.obj'
FilInc='d:\Projects\ReadExcitedStates\ReadExcitedStates\BuildNumber.fi'		!<--    Source\BuildNumber.fi ???
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!

INQUIRE(file=FilObj,EXIST=res)
If (.not.res) Return

INQUIRE(file=FilInc,EXIST=res)
If (.not.res) Return

Open(98,file=FilInc)
Read(98,'(10x,i6)')nBuild0
If (nBuild0/=nbuild) Return
nbuild=nbuild+1
Rewind(98)
Write(98,'(''nBuild   ='',i6)')nbuild
Close(98)

Open(98,file=FilObj)
Close(98,status='delete')

!str='del '//FilObj
!res=systemqq(str)

End

max_stars_repo_name: skignatov/ReadExcitedStates

max_issues_repo_name: skignatov/ReadExcitedStates

max_forks_repo_name: skignatov/ReadExcitedStates


text:
*> \brief \b CEBCHVXX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*      SUBROUTINE CEBCHVXX( THRESH, PATH )
*
*     .. Scalar Arguments ..
*      REAL               THRESH
*      CHARACTER*3        PATH
*       ..
*
*  Purpose
*  ======
*
*> \details \b Purpose:
*> \verbatim
*>
*>  CEBCHVXX will run CGESVXX on a series of Hilbert matrices and then
*>  compare the error bounds returned by CGESVXX to see if the returned
*>  answer indeed falls within those bounds.
*>
*>  Eight test ratios will be computed.  The tests will pass if they are .LT.
*>  THRESH.  There are two cases that are determined by 1 / (SQRT( N ) * EPS).
*>  If that value is .LE. to the component wise reciprocal condition number,
*>  it uses the guaranteed case, other wise it uses the unguaranteed case.
*>
*>  Test ratios:
*>     Let Xc be X_computed and Xt be X_truth.
*>     The norm used is the infinity norm.
*>
*>     Let A be the guaranteed case and B be the unguaranteed case.
*>
*>       1. Normwise guaranteed forward error bound.
*>       A: norm ( abs( Xc - Xt ) / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ) and
*>          ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N),10) * EPS.
*>          If these conditions are met, the test ratio is set to be
*>          ERRBND( *, nwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS.
*>       B: For this case, CGESVXX should just return 1.  If it is less than
*>          one, treat it the same as in 1A.  Otherwise it fails. (Set test
*>          ratio to ERRBND( *, nwise_i, bnd_i ) * THRESH?)
*>
*>       2. Componentwise guaranteed forward error bound.
*>       A: norm ( abs( Xc(j) - Xt(j) ) ) / norm (Xt(j)) .LE. ERRBND( *, cwise_i, bnd_i )
*>          for all j .AND. ERRBND( *, cwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS.
*>          If these conditions are met, the test ratio is set to be
*>          ERRBND( *, cwise_i, bnd_i ) / MAX(SQRT(N), 10).  Otherwise it is 1/EPS.
*>       B: Same as normwise test ratio.
*>
*>       3. Backwards error.
*>       A: The test ratio is set to BERR/EPS.
*>       B: Same test ratio.
*>
*>       4. Reciprocal condition number.
*>       A: A condition number is computed with Xt and compared with the one
*>          returned from CGESVXX.  Let RCONDc be the RCOND returned by CGESVXX
*>          and RCONDt be the RCOND from the truth value.  Test ratio is set to
*>          MAX(RCONDc/RCONDt, RCONDt/RCONDc).
*>       B: Test ratio is set to 1 / (EPS * RCONDc).
*>
*>       5. Reciprocal normwise condition number.
*>       A: The test ratio is set to
*>          MAX(ERRBND( *, nwise_i, cond_i ) / NCOND, NCOND / ERRBND( *, nwise_i, cond_i )).
*>       B: Test ratio is set to 1 / (EPS * ERRBND( *, nwise_i, cond_i )).
*>
*>       6. Reciprocal componentwise condition number.
*>       A: Test ratio is set to
*>          MAX(ERRBND( *, cwise_i, cond_i ) / CCOND, CCOND / ERRBND( *, cwise_i, cond_i )).
*>       B: Test ratio is set to 1 / (EPS * ERRBND( *, cwise_i, cond_i )).
*>
*>     .. Parameters ..
*>     NMAX is determined by the largest number in the inverse of the hilbert
*>     matrix.  Precision is exhausted when the largest entry in it is greater
*>     than 2 to the power of the number of bits in the fraction of the data
*>     type used plus one, which is 24 for single precision.
*>     NMAX should be 6 for single and 11 for double.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup complex_lin
*
*  =====================================================================
      SUBROUTINE CEBCHVXX( THRESH, PATH )
      IMPLICIT NONE
*     .. Scalar Arguments ..
      REAL               THRESH
      CHARACTER*3        PATH

      INTEGER            NMAX, NPARAMS, NERRBND, NTESTS, KL, KU
      PARAMETER          (NMAX = 6, NPARAMS = 2, NERRBND = 3,
     $                    NTESTS = 6)

*     .. Local Scalars ..
      INTEGER            N, NRHS, INFO, I ,J, k, NFAIL, LDA,
     $                   N_AUX_TESTS, LDAB, LDAFB
      CHARACTER          FACT, TRANS, UPLO, EQUED
      CHARACTER*2        C2
      CHARACTER(3)       NGUAR, CGUAR
      LOGICAL            printed_guide
      REAL               NCOND, CCOND, M, NORMDIF, NORMT, RCOND,
     $                   RNORM, RINORM, SUMR, SUMRI, EPS,
     $                   BERR(NMAX), RPVGRW, ORCOND,
     $                   CWISE_ERR, NWISE_ERR, CWISE_BND, NWISE_BND,
     $                   CWISE_RCOND, NWISE_RCOND,
     $                   CONDTHRESH, ERRTHRESH
      COMPLEX            ZDUM

*     .. Local Arrays ..
      REAL               TSTRAT(NTESTS), RINV(NMAX), PARAMS(NPARAMS),
     $                   S(NMAX), R(NMAX),C(NMAX),RWORK(3*NMAX),
     $                   DIFF(NMAX, NMAX),
     $                   ERRBND_N(NMAX*3), ERRBND_C(NMAX*3)
      INTEGER            IPIV(NMAX)
      COMPLEX            A(NMAX,NMAX),INVHILB(NMAX,NMAX),X(NMAX,NMAX),
     $                   WORK(NMAX*3*5), AF(NMAX, NMAX),B(NMAX, NMAX),
     $                   ACOPY(NMAX, NMAX),
     $                   AB( (NMAX-1)+(NMAX-1)+1, NMAX ),
     $                   ABCOPY( (NMAX-1)+(NMAX-1)+1, NMAX ),
     $                   AFB( 2*(NMAX-1)+(NMAX-1)+1, NMAX )

*     .. External Functions ..
      REAL               SLAMCH

*     .. External Subroutines ..
      EXTERNAL           CLAHILB, CGESVXX, CSYSVXX, CPOSVXX,
     $                   CGBSVXX, CLACPY, LSAMEN
      LOGICAL            LSAMEN

*     .. Intrinsic Functions ..
      INTRINSIC          SQRT, MAX, ABS, REAL, AIMAG

*     .. Statement Functions ..
      REAL               CABS1
*     ..
*     .. Statement Function Definitions ..
      CABS1( ZDUM ) = ABS( REAL( ZDUM ) ) + ABS( AIMAG( ZDUM ) )

*     .. Parameters ..
      INTEGER            NWISE_I, CWISE_I
      PARAMETER          (NWISE_I = 1, CWISE_I = 1)
      INTEGER            BND_I, COND_I
      PARAMETER          (BND_I = 2, COND_I = 3)

*  Create the loop to test out the Hilbert matrices

      FACT = 'E'
      UPLO = 'U'
      TRANS = 'N'
      EQUED = 'N'
      EPS = SLAMCH('Epsilon')
      NFAIL = 0
      N_AUX_TESTS = 0
      LDA = NMAX
      LDAB = (NMAX-1)+(NMAX-1)+1
      LDAFB = 2*(NMAX-1)+(NMAX-1)+1
      C2 = PATH( 2: 3 )

*     Main loop to test the different Hilbert Matrices.

      printed_guide = .false.

      DO N = 1 , NMAX
         PARAMS(1) = -1
         PARAMS(2) = -1

         KL = N-1
         KU = N-1
         NRHS = n
         M = MAX(SQRT(REAL(N)), 10.0)

*        Generate the Hilbert matrix, its inverse, and the
*        right hand side, all scaled by the LCM(1,..,2N-1).
         CALL CLAHILB(N, N, A, LDA, INVHILB, LDA, B,
     $        LDA, WORK, INFO, PATH)

*        Copy A into ACOPY.
         CALL CLACPY('ALL', N, N, A, NMAX, ACOPY, NMAX)

*        Store A in band format for GB tests
         DO J = 1, N
            DO I = 1, KL+KU+1
               AB( I, J ) = (0.0E+0,0.0E+0)
            END DO
         END DO
         DO J = 1, N
            DO I = MAX( 1, J-KU ), MIN( N, J+KL )
               AB( KU+1+I-J, J ) = A( I, J )
            END DO
         END DO

*        Copy AB into ABCOPY.
         DO J = 1, N
            DO I = 1, KL+KU+1
               ABCOPY( I, J ) = (0.0E+0,0.0E+0)
            END DO
         END DO
         CALL CLACPY('ALL', KL+KU+1, N, AB, LDAB, ABCOPY, LDAB)

*        Call C**SVXX with default PARAMS and N_ERR_BND = 3.
         IF ( LSAMEN( 2, C2, 'SY' ) ) THEN
            CALL CSYSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
     $           IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
     $           RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
     $           PARAMS, WORK, RWORK, INFO)
         ELSE IF ( LSAMEN( 2, C2, 'PO' ) ) THEN
            CALL CPOSVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
     $           EQUED, S, B, LDA, X, LDA, ORCOND,
     $           RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
     $           PARAMS, WORK, RWORK, INFO)
         ELSE IF ( LSAMEN( 2, C2, 'HE' ) ) THEN
            CALL CHESVXX(FACT, UPLO, N, NRHS, ACOPY, LDA, AF, LDA,
     $           IPIV, EQUED, S, B, LDA, X, LDA, ORCOND,
     $           RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
     $           PARAMS, WORK, RWORK, INFO)
         ELSE IF ( LSAMEN( 2, C2, 'GB' ) ) THEN
            CALL CGBSVXX(FACT, TRANS, N, KL, KU, NRHS, ABCOPY,
     $           LDAB, AFB, LDAFB, IPIV, EQUED, R, C, B,
     $           LDA, X, LDA, ORCOND, RPVGRW, BERR, NERRBND,
     $           ERRBND_N, ERRBND_C, NPARAMS, PARAMS, WORK, RWORK,
     $           INFO)
         ELSE
            CALL CGESVXX(FACT, TRANS, N, NRHS, ACOPY, LDA, AF, LDA,
     $           IPIV, EQUED, R, C, B, LDA, X, LDA, ORCOND,
     $           RPVGRW, BERR, NERRBND, ERRBND_N, ERRBND_C, NPARAMS,
     $           PARAMS, WORK, RWORK, INFO)
         END IF

         N_AUX_TESTS = N_AUX_TESTS + 1
         IF (ORCOND .LT. EPS) THEN
!        Either factorization failed or the matrix is flagged, and 1 <=
!        INFO <= N+1. We don't decide based on rcond anymore.
!            IF (INFO .EQ. 0 .OR. INFO .GT. N+1) THEN
!               NFAIL = NFAIL + 1
!               WRITE (*, FMT=8000) N, INFO, ORCOND, RCOND
!            END IF
         ELSE
!        Either everything succeeded (INFO == 0) or some solution failed
!        to converge (INFO > N+1).
            IF (INFO .GT. 0 .AND. INFO .LE. N+1) THEN
               NFAIL = NFAIL + 1
               WRITE (*, FMT=8000) C2, N, INFO, ORCOND, RCOND
            END IF
         END IF

*        Calculating the difference between C**SVXX's X and the true X.
         DO I = 1,N
            DO J =1,NRHS
               DIFF(I,J) = X(I,J) - INVHILB(I,J)
            END DO
         END DO

*        Calculating the RCOND
         RNORM = 0
         RINORM = 0
         IF ( LSAMEN( 2, C2, 'PO' ) .OR. LSAMEN( 2, C2, 'SY' ) .OR.
     $        LSAMEN( 2, C2, 'HE' ) ) THEN
            DO I = 1, N
               SUMR = 0
               SUMRI = 0
               DO J = 1, N
                  SUMR = SUMR + S(I) * CABS1(A(I,J)) * S(J)
                  SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (S(J) * S(I))
               END DO
               RNORM = MAX(RNORM,SUMR)
               RINORM = MAX(RINORM,SUMRI)
            END DO
         ELSE IF ( LSAMEN( 2, C2, 'GE' ) .OR. LSAMEN( 2, C2, 'GB' ) )
     $           THEN
            DO I = 1, N
               SUMR = 0
               SUMRI = 0
               DO J = 1, N
                  SUMR = SUMR + R(I) * CABS1(A(I,J)) * C(J)
                  SUMRI = SUMRI + CABS1(INVHILB(I, J)) / (R(J) * C(I))
               END DO
               RNORM = MAX(RNORM,SUMR)
               RINORM = MAX(RINORM,SUMRI)
            END DO
         END IF

         RNORM = RNORM / CABS1(A(1, 1))
         RCOND = 1.0/(RNORM * RINORM)

*        Calculating the R for normwise rcond.
         DO I = 1, N
            RINV(I) = 0.0
         END DO
         DO J = 1, N
            DO I = 1, N
               RINV(I) = RINV(I) + CABS1(A(I,J))
            END DO
         END DO

*        Calculating the Normwise rcond.
         RINORM = 0.0
         DO I = 1, N
            SUMRI = 0.0
            DO J = 1, N
               SUMRI = SUMRI + CABS1(INVHILB(I,J) * RINV(J))
            END DO
            RINORM = MAX(RINORM, SUMRI)
         END DO

!        invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
!        by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
         NCOND = CABS1(A(1,1)) / RINORM

         CONDTHRESH = M * EPS
         ERRTHRESH = M * EPS

         DO K = 1, NRHS
            NORMT = 0.0
            NORMDIF = 0.0
            CWISE_ERR = 0.0
            DO I = 1, N
               NORMT = MAX(CABS1(INVHILB(I, K)), NORMT)
               NORMDIF = MAX(CABS1(X(I,K) - INVHILB(I,K)), NORMDIF)
               IF (INVHILB(I,K) .NE. 0.0) THEN
                  CWISE_ERR = MAX(CABS1(X(I,K) - INVHILB(I,K))
     $                            /CABS1(INVHILB(I,K)), CWISE_ERR)
               ELSE IF (X(I, K) .NE. 0.0) THEN
                  CWISE_ERR = SLAMCH('OVERFLOW')
               END IF
            END DO
            IF (NORMT .NE. 0.0) THEN
               NWISE_ERR = NORMDIF / NORMT
            ELSE IF (NORMDIF .NE. 0.0) THEN
               NWISE_ERR = SLAMCH('OVERFLOW')
            ELSE
               NWISE_ERR = 0.0
            ENDIF

            DO I = 1, N
               RINV(I) = 0.0
            END DO
            DO J = 1, N
               DO I = 1, N
                  RINV(I) = RINV(I) + CABS1(A(I, J) * INVHILB(J, K))
               END DO
            END DO
            RINORM = 0.0
            DO I = 1, N
               SUMRI = 0.0
               DO J = 1, N
                  SUMRI = SUMRI
     $                 + CABS1(INVHILB(I, J) * RINV(J) / INVHILB(I, K))
               END DO
               RINORM = MAX(RINORM, SUMRI)
            END DO
!        invhilb is the inverse *unscaled* Hilbert matrix, so scale its norm
!        by 1/A(1,1) to make the scaling match A (the scaled Hilbert matrix)
            CCOND = CABS1(A(1,1))/RINORM

!        Forward error bound tests
            NWISE_BND = ERRBND_N(K + (BND_I-1)*NRHS)
            CWISE_BND = ERRBND_C(K + (BND_I-1)*NRHS)
            NWISE_RCOND = ERRBND_N(K + (COND_I-1)*NRHS)
            CWISE_RCOND = ERRBND_C(K + (COND_I-1)*NRHS)
!            write (*,*) 'nwise : ', n, k, ncond, nwise_rcond,
!     $           condthresh, ncond.ge.condthresh
!            write (*,*) 'nwise2: ', k, nwise_bnd, nwise_err, errthresh
            IF (NCOND .GE. CONDTHRESH) THEN
               NGUAR = 'YES'
               IF (NWISE_BND .GT. ERRTHRESH) THEN
                  TSTRAT(1) = 1/(2.0*EPS)
               ELSE
                  IF (NWISE_BND .NE. 0.0) THEN
                     TSTRAT(1) = NWISE_ERR / NWISE_BND
                  ELSE IF (NWISE_ERR .NE. 0.0) THEN
                     TSTRAT(1) = 1/(16.0*EPS)
                  ELSE
                     TSTRAT(1) = 0.0
                  END IF
                  IF (TSTRAT(1) .GT. 1.0) THEN
                     TSTRAT(1) = 1/(4.0*EPS)
                  END IF
               END IF
            ELSE
               NGUAR = 'NO'
               IF (NWISE_BND .LT. 1.0) THEN
                  TSTRAT(1) = 1/(8.0*EPS)
               ELSE
                  TSTRAT(1) = 1.0
               END IF
            END IF
!            write (*,*) 'cwise : ', n, k, ccond, cwise_rcond,
!     $           condthresh, ccond.ge.condthresh
!            write (*,*) 'cwise2: ', k, cwise_bnd, cwise_err, errthresh
            IF (CCOND .GE. CONDTHRESH) THEN
               CGUAR = 'YES'
               IF (CWISE_BND .GT. ERRTHRESH) THEN
                  TSTRAT(2) = 1/(2.0*EPS)
               ELSE
                  IF (CWISE_BND .NE. 0.0) THEN
                     TSTRAT(2) = CWISE_ERR / CWISE_BND
                  ELSE IF (CWISE_ERR .NE. 0.0) THEN
                     TSTRAT(2) = 1/(16.0*EPS)
                  ELSE
                     TSTRAT(2) = 0.0
                  END IF
                  IF (TSTRAT(2) .GT. 1.0) TSTRAT(2) = 1/(4.0*EPS)
               END IF
            ELSE
               CGUAR = 'NO'
               IF (CWISE_BND .LT. 1.0) THEN
                  TSTRAT(2) = 1/(8.0*EPS)
               ELSE
                  TSTRAT(2) = 1.0
               END IF
            END IF

!     Backwards error test
            TSTRAT(3) = BERR(K)/EPS

!     Condition number tests
            TSTRAT(4) = RCOND / ORCOND
            IF (RCOND .GE. CONDTHRESH .AND. TSTRAT(4) .LT. 1.0)
     $         TSTRAT(4) = 1.0 / TSTRAT(4)

            TSTRAT(5) = NCOND / NWISE_RCOND
            IF (NCOND .GE. CONDTHRESH .AND. TSTRAT(5) .LT. 1.0)
     $         TSTRAT(5) = 1.0 / TSTRAT(5)

            TSTRAT(6) = CCOND / NWISE_RCOND
            IF (CCOND .GE. CONDTHRESH .AND. TSTRAT(6) .LT. 1.0)
     $         TSTRAT(6) = 1.0 / TSTRAT(6)

            DO I = 1, NTESTS
               IF (TSTRAT(I) .GT. THRESH) THEN
                  IF (.NOT.PRINTED_GUIDE) THEN
                     WRITE(*,*)
                     WRITE( *, 9996) 1
                     WRITE( *, 9995) 2
                     WRITE( *, 9994) 3
                     WRITE( *, 9993) 4
                     WRITE( *, 9992) 5
                     WRITE( *, 9991) 6
                     WRITE( *, 9990) 7
                     WRITE( *, 9989) 8
                     WRITE(*,*)
                     PRINTED_GUIDE = .TRUE.
                  END IF
                  WRITE( *, 9999) C2, N, K, NGUAR, CGUAR, I, TSTRAT(I)
                  NFAIL = NFAIL + 1
               END IF
            END DO
      END DO

c$$$         WRITE(*,*)
c$$$         WRITE(*,*) 'Normwise Error Bounds'
c$$$         WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,nwise_i,bnd_i)
c$$$         WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,nwise_i,cond_i)
c$$$         WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,nwise_i,rawbnd_i)
c$$$         WRITE(*,*)
c$$$         WRITE(*,*) 'Componentwise Error Bounds'
c$$$         WRITE(*,*) 'Guaranteed error bound: ',ERRBND(NRHS,cwise_i,bnd_i)
c$$$         WRITE(*,*) 'Reciprocal condition number: ',ERRBND(NRHS,cwise_i,cond_i)
c$$$         WRITE(*,*) 'Raw error estimate: ',ERRBND(NRHS,cwise_i,rawbnd_i)
c$$$         print *, 'Info: ', info
c$$$         WRITE(*,*)
*         WRITE(*,*) 'TSTRAT: ',TSTRAT

      END DO

      WRITE(*,*)
      IF( NFAIL .GT. 0 ) THEN
         WRITE(*,9998) C2, NFAIL, NTESTS*N+N_AUX_TESTS
      ELSE
         WRITE(*,9997) C2
      END IF
 9999 FORMAT( ' C', A2, 'SVXX: N =', I2, ', RHS = ', I2,
     $     ', NWISE GUAR. = ', A, ', CWISE GUAR. = ', A,
     $     ' test(',I1,') =', G12.5 )
 9998 FORMAT( ' C', A2, 'SVXX: ', I6, ' out of ', I6,
     $     ' tests failed to pass the threshold' )
 9997 FORMAT( ' C', A2, 'SVXX passed the tests of error bounds' )
*     Test ratios.
 9996 FORMAT( 3X, I2, ': Normwise guaranteed forward error', / 5X,
     $     'Guaranteed case: if norm ( abs( Xc - Xt )',
     $     ' / norm ( Xt ) .LE. ERRBND( *, nwise_i, bnd_i ), then',
     $     / 5X,
     $     'ERRBND( *, nwise_i, bnd_i ) .LE. MAX(SQRT(N), 10) * EPS')
 9995 FORMAT( 3X, I2, ': Componentwise guaranteed forward error' )
 9994 FORMAT( 3X, I2, ': Backwards error' )
 9993 FORMAT( 3X, I2, ': Reciprocal condition number' )
 9992 FORMAT( 3X, I2, ': Reciprocal normwise condition number' )
 9991 FORMAT( 3X, I2, ': Raw normwise error estimate' )
 9990 FORMAT( 3X, I2, ': Reciprocal componentwise condition number' )
 9989 FORMAT( 3X, I2, ': Raw componentwise error estimate' )

 8000 FORMAT( ' C', A2, 'SVXX: N =', I2, ', INFO = ', I3,
     $     ', ORCOND = ', G12.5, ', real RCOND = ', G12.5 )

      END

max_stars_repo_name: sydha/Lapack

max_issues_repo_name: sydha/Lapack

max_forks_repo_name: sydha/Lapack


text:
#include "../defines.inc"
module nucleosome
!   --------------------------------------------------------------
!
!    This module is designed to handle the various geometrical
!    considerations of nucleosomes.
!
!   --------------------------------------------------------------

   use precision, only: dp, pi
   implicit none
   real(dp), parameter :: basePairsPerTurn = 10.5_dp
   real(dp), dimension(10, WLC_P__MAX_BACEPAIRS_PER_BEAD) :: multiParams
   real(dp), dimension(3, 3, 147) :: nucleosomeROT
   real(dp), dimension(3, 147) :: nucleosomeTran

   private :: nucleosomeROT, nucleosomeTran
contains

   subroutine nucleosome_prop(Uin, Vin, Rin, linkBP, wrapBP, Uout, Vout, Rout)
!Calculates the final position and orientation
!Rout,Uout,Vout
!based on the the incomming position and orientation
!Rin, Uin, Vin
!for a nucleosome with wrapBP bace paris of DNA wrapped around it.
!The intrisic rotation of a linkBP long linker that follows is also included.
      use vector_utils, only: cross
      implicit none
      real(dp), intent(in), dimension(3) :: Uin ! Entry direction along DNA
      real(dp), intent(in), dimension(3) :: Vin ! Entry tangent vector
      real(dp), intent(in), dimension(3) :: Rin ! Position of entry
      real(dp), intent(in) :: linkBP
      integer, intent(in) :: wrapBP
      real(dp), intent(out), dimension(3) :: Uout ! Exit position
      real(dp), intent(out), dimension(3) :: Vout ! Exit tangent vector
      real(dp), intent(out), dimension(3) :: Rout ! Exit position

      real(dp), dimension(3, 3) :: linkRot
      real(dp), dimension(3, 3) :: mtrx
      real(dp), parameter :: angle = 2*pi/basePairsPerTurn ! intrinsic rotation/bp

      mtrx(:, 1) = Vin
      mtrx(:, 2) = cross(Uin, Vin)
      mtrx(:, 3) = Uin

      Rout = Rin + MATMUL(mtrx, nucleosomeTran(:, wrapBP))

      linkRot(1, :) = [cos(angle*linkBP), -sin(angle*linkBP), 0.0_dp]
      linkRot(2, :) = [sin(angle*linkBP), cos(angle*linkBP), 0.0_dp]
      linkRot(3, :) = [0.0_dp, 0.0_dp, 1.0_dp]

      mtrx = MATMUL(MATMUL(mtrx, nucleosomeROT(:, :, wrapBP)), linkRot)

      Uout = mtrx(:, 3)/norm2(mtrx(:, 3))
      Vout = mtrx(:, 1)/norm2(mtrx(:, 1))
      return

   end subroutine nucleosome_prop

   function nucleosome_energy(RP1, R, UP1, U, VP1, V, linkBP, wrapBP)
!Caltulate the bending energy of a nucleosome
!at position R and orientation U, V
!with wrapBP of DNA wrapped around it followed by a liner linkBP
!bace pairs long that ends at RP1, UP1, and VP1.
!Assumes a SSWLCWT linker.
      use MC_wlc, only: E_SSWLCWT
      real(dp), intent(in), dimension(3) :: R ! R of bead i
      real(dp), intent(in), dimension(3) :: RP1 ! R of bead i+1
      real(dp), intent(in), dimension(3) :: U ! U of bead i
      real(dp), intent(in), dimension(3) :: UP1 ! U of bead i+1
      real(dp), intent(in), dimension(3) :: V ! V of bead i
      real(dp), intent(in), dimension(3) :: VP1 ! V of bead i+1
      real(dp), intent(in) :: linkBP
      integer, intent(in) :: wrapBP
      real(dp) nucleosome_energy(4)
      real(dp) EB, EPAR, EPERP, GAM, ETA, XIR, XIU, sigma, etwist, simtype

      real(dp) Rtemp(3)
      real(dp) Utemp(3)
      real(dp) Vtemp(3)

      call nucleosome_prop(U, V, R, linkBP, wrapBP, Utemp, Vtemp, Rtemp)

      !  need to interpolate between basepairs if fractional
      call get_params(linkBP, EB, EPAR, EPERP, GAM, ETA, XIR, XIU, sigma, etwist, simtype)

      nucleosome_energy = E_SSWLCWT(RP1, Rtemp, UP1, Utemp, VP1, Vtemp, &
                                    EB, &!multiParams(1,linkBP),&   ! EB
                                    EPAR, &!multiParams(2,linkBP),&   ! EPAR
                                    EPERP, &!multiParams(3,linkBP),&   ! EPERP
                                    ETA, &!multiParams(5,linkBP),&   ! ETA
                                    GAM, &!multiParams(4,linkBP),&   ! GAM
                                    ETWIST)!multiParams(9,linkBP))   ! etwist
   end function nucleosome_energy

   function internucleosome_energy(RI, RJ, UI, UJ, VI, VJ)
! internucleosome attraction energy. this is based off the work out of
! de pablos group where they look at the pairwise potential between 2
! nucleosome dependent on their orientation/geometry and distance
! from each other. parameters come from their counterion condensation
! assumptions. i try to simplify their findings and create 3 main
! classes of interactions: 1) face-face where the nucleosomes are
! aligned such that their core histones can interact 2) face-side
! where the nucleosomes are aligned perpendicular such that the core
! histones of one nucleosome interacts with the DNA wrapping the exterior
! of the other nucleosome and 3) side-side where the nucleosome DNA wrappings
! are in line with each other.
      use MC_wlc, only: E_SSWLCWT
      use vector_utils, only: cross
      use energies, only: energyOf, internucleosome_
      real(dp), intent(in), dimension(3) :: RI ! R of nuc i
      real(dp), intent(in), dimension(3) :: RJ! R of nuc j
      real(dp), intent(in), dimension(3) :: UI ! U of nuc i
      real(dp), intent(in), dimension(3) :: UJ ! U of nuc j
      real(dp), intent(in), dimension(3) :: VI ! V of nuc i
      real(dp), intent(in), dimension(3) :: VJ ! V of nuc j
      real(dp), parameter :: tau_faceface = 1.38
      real(dp), parameter :: e_faceface = 3.712*WLC_P__INTERNUCLEOSOME_ENERGY
      real(dp), parameter :: tau_faceside = 0.82
      real(dp), parameter :: e_faceside = 1.476*WLC_P__INTERNUCLEOSOME_ENERGY
      real(dp), parameter :: tau_sideside = 2.0
      real(dp), parameter :: e_sideside = 1.64*WLC_P__INTERNUCLEOSOME_ENERGY
      real(dp), dimension(3), parameter :: center = [4.8455, -2.4445, 0.6694]
      real(dp), dimension(3, 3) :: mtrxI, mtrxJ
      real(dp), dimension(3) :: polyI, faceI, faceItop, faceIbot
      real(dp), dimension(3) :: polyJ, faceJ, faceJtop, faceJbot
      real(dp), dimension(4, 3) :: faceDistList
      real(dp), dimension(4) :: distList
      real(dp), dimension(3) :: distS, distC, dist
      real(dp) cospsi, costhetaI, costhetaJ, cosphiI, cosphiJ
      integer i, indList(1)
      real(dp) internucleosome_energy

      ! initialiaze
      internucleosome_energy = 0

      ! construct matrices
      mtrxI(:, 1) = VI
      mtrxI(:, 2) = cross(UI, VI)
      mtrxI(:, 3) = UI
      mtrxJ(:, 1) = VJ
      mtrxJ(:, 2) = cross(UJ, VJ)
      mtrxJ(:, 3) = UJ

      ! center of nucs
      polyI = RI + MATMUL(mtrxI, center)
      polyJ = RJ +  MATMUL(mtrxJ, center)

      ! construct face I normal vector (pointing up through face)
      faceItop = RI + MATMUL(mtrxI, center + [0.0_dp, WLC_P__NUCLEOSOME_HEIGHT/2,0.0_dp])
      faceIbot = RI + MATMUL(mtrxI, center + [0.0_dp, -WLC_P__NUCLEOSOME_HEIGHT/2,0.0_dp])
      faceI = faceItop - faceIbot
      ! construct face J normal vector (pointing up through face)
      faceJtop = RJ + MATMUL(mtrxJ, center + [0.0_dp, WLC_P__NUCLEOSOME_HEIGHT/2,0.0_dp])
      faceJbot = RJ + MATMUL(mtrxJ, center + [0.0_dp, -WLC_P__NUCLEOSOME_HEIGHT/2,0.0_dp])
      faceJ = faceJtop - faceJbot
      cospsi = dot_product(faceI/norm2(faceI), faceJ/norm2(faceJ))

      ! list of combinatorial face attractions
      faceDistList(1, :) = faceItop - faceJbot
      faceDistList(2, :) = faceItop - faceJtop
      faceDistList(3, :) = faceIbot - faceJbot
      faceDistList(4, :) = faceIbot - faceJtop

      ! find the closest faces to define the face oritentation vector
      do i = 1, 4
         distList(i) = norm2(faceDistList(i, :))
      enddo
      indList = minloc(distList)
      distS = faceDistList(indList(1), :)
      costhetaI = dot_product(distS/norm2(distS), faceI/norm2(faceI))
      costhetaJ = dot_product(distS/norm2(distS), faceJ/norm2(faceJ))

      ! face-face (histone-histone attraction)
      if (norm2(distS) <= tau_faceface) then
         internucleosome_energy = internucleosome_energy &
                                  - e_faceface*(costhetaI**2)*(costhetaJ**2)/tau_faceface
      else
         internucleosome_energy = internucleosome_energy &
                                  - e_faceface*(costhetaI**2)*(costhetaJ**2)/norm2(distS)
      endif
      distC = polyI - polyJ
      cosphiI = dot_product(distC/norm2(distC), faceI/norm2(faceI))
      cosphiJ = dot_product(distC/norm2(distC), faceJ/norm2(faceJ))

      ! face-side (histone-DNA wrapping attraction)
      dist = norm2(distC) - WLC_P__NUCLEOSOME_HEIGHT/2 - WLC_P__NUCLEOSOME_RADIUS
      if (norm2(dist) <= tau_faceside) then
         internucleosome_energy = internucleosome_energy &
                                  - e_faceside*(1 - cospsi**2)*(cosphiI**2 + cosphiJ**2)/tau_faceside
      else
         internucleosome_energy = internucleosome_energy &
                                  - e_faceside*(1 - cospsi**2)*(cosphiI**2 + cosphiJ**2)/norm2(dist)
      endif

      ! side-side (DNA wrapping-DNA wrapping attraction)
      dist = norm2(distC)-2*WLC_P__NUCLEOSOME_RADIUS
      if (norm2(dist) <= tau_sideside) then
         internucleosome_energy = internucleosome_energy &
                                  - e_sideside*(1 - cosphiI**2)*(1 - cosphiJ**2)/tau_sideside
      else
         internucleosome_energy = internucleosome_energy &
                                  - e_sideside*(1 - cosphiI**2)*(1 - cosphiJ**2)/norm2(dist)
      endif

   end function internucleosome_energy

! ---------------------------------------------------------------------
!  Return parameters for a linker i bace pairs long
! ---------------------------------------------------------------------
   subroutine get_params(i, EB, EPAR, EPERP, GAM, ETA, XIR, XIU, sigma, etwist, simtype)
!Return parameters for a linker that is i bace pairs long
      implicit none
      real(dp), intent(in) :: i ! Number of bace pairs in linder
      real(dp), intent(out) :: EB ! Bending modulus
      real(dp), intent(out) :: EPAR ! Stretch modulus
      real(dp), intent(out) :: EPERP ! Shear modulus
      real(dp), intent(out) :: GAM ! ground state segment length
      real(dp), intent(out) :: ETA ! bend-shear coupling
      real(dp), intent(out) :: XIR, XIU, sigma
      real(dp), intent(out) :: etwist ! Twist coefficient
      real(dp), intent(out) :: simtype ! WLC, SSWLC, or Gaussian Chain
      integer indUp, indDown
      real(dp) ratio, offratio

      ! interpolate between points
      indDown = floor(i)
      indUp = ceiling(i)
      ratio = i/indUp
      offratio = 1 - ratio
      EB = ratio*multiParams(1, indUp) + offratio*multiParams(1, indDown)
      EPAR = ratio*multiParams(2, indUp) + offratio*multiParams(2, indDown)
      EPERP = ratio*multiParams(3, indUp) + offratio*multiParams(3, indDown)
      GAM = ratio*multiParams(4, indUp) + offratio*multiParams(4, indDown)
      ETA = ratio*multiParams(5, indUp) + offratio*multiParams(5, indDown)
      XIR = ratio*multiParams(6, indUp) + offratio*multiParams(6, indDown)
      XIU = ratio*multiParams(7, indUp) + offratio*multiParams(7, indDown)
      sigma = ratio*multiParams(8, indUp) + offratio*multiParams(8, indDown)
      etwist = ratio*multiParams(9, indUp) + offratio*multiParams(9, indDown)
      simtype = int(ratio*multiParams(10, indUp) + offratio*multiParams(10, indDown) + 0.1_dp) ! the +0.1 is to prevent round down due to lack of precision
   end subroutine get_params

! ----------------------------------------------------------------------
!
!          Look up elastic constants for linkers,
!          Look up translation and rotation for nucleosomes.
!
! ----------------------------------------------------------------------
   subroutine setup_nucleosome_constants()
      use precision, only: nan
      use MC_wlc, only: calc_elastic_constants
      implicit none
      integer i, j, simtype
      real(dp) DEL
      real(dp) EB, EPAR, EPERP, GAM, ETA, XIR, XIU, sigma, etwist, dt

      multiParams = nan
      do i = 2, WLC_P__MAX_BACEPAIRS_PER_BEAD
         DEL = WLC_P__LENGTH_PER_BP*i/WLC_P__LP
         call calc_elastic_constants(DEL, WLC_P__LP, WLC_P__LT, EB, EPAR, GAM, XIR, EPERP, ETA, XIU, DT, &
                                     SIGMA, ETWIST, simtype)
         if (simtype .ne. 2) then
            print *, "Error!  nucleosomes only setup for SSWLCWT"
            print *, "Simtype ", simtype, " encountered"
            stop
         endif
         multiParams(1, i) = EB
         multiParams(2, i) = EPAR
         multiParams(3, i) = EPERP
         multiParams(4, i) = GAM
         multiParams(5, i) = ETA
         multiParams(6, i) = XIR
         multiParams(7, i) = XIU
         multiParams(8, i) = sigma
         multiParams(9, i) = etwist
         multiParams(10, i) = real(simtype, dp)

      enddo

      open (UNIT=5, FILE="input/nucleosomeR", STATUS="OLD")
      do i = 1, 147
         do j = 1, 3
            read (5, *) nucleosomeROT(j, :, 148 - i)
         enddo
      enddo
      close (5)

      if (WLC_P__INCLUDE_NUC_TRANS) then
         open (UNIT=5, FILE="input/nucleosomeT", STATUS="OLD")
         do i = 1, 147
            read (5, *) nucleosomeTran(:, 148 - i)
         enddo
         close (5)
      else
         nucleosomeTran = 0.0_dp
      endif
   end subroutine setup_nucleosome_constants

   subroutine load_nucleosome_positions(wlc_nucleosomeWrap, wlc_basepairs)
! if WLC_P__INCLUDE_DISCRETIZE_LINKER is turned on, then this determines the
! discretization scheme of the beads throughout the chain. it sets a real
! (can be non-integer) value for the discretization as a function of contour length
! and WLC_P__NB. if WLC_P__INCLUDE_DISCRETIZE_LINKER is off, then each bead is
! if WLC_P__INCLUDE_DISCRETIZE_LINKER is turn on, then this determines the
! discretization scheme of the beads throughout the chain. it defaults to trying
! to set integer values for each bead discretization, but not has the capability to
! use real numbers. if WLC_P__INCLUDE_DISCRETIZE_LINKER is off, then each bead is
! default assumed to be a nucleosome, so there are no additional beads used to model
! the fluctuations of the linker geometry with more detail. nucleosomes are intialized
! on the chain according to WLC_P__LINKER_TYPE, where the default is 'phased', i.e.
! the nuclesomes are separated by a constant linker length throughout the chain
      use precision, only: nan
      ! sterics testing !
      use GJKAlgorithm, only: GJK, sameShapeTest, noIntersectX, intersectX, tangentX, runtimeTest5, runtimeTest6, &
                              noIntersectY, intersectY, tangentY, runtimeTest1, runtimeTest3, &
                              noIntersectZ, intersectZ, tangentZ, runtimeTest2, runtimeTest4
      use polydispersity, only: first_bead_of_chain, last_bead_of_chain
      implicit none
      integer, intent(out) :: wlc_nucleosomeWrap(WLC_P__NT)
      real(dp), intent(out) :: wlc_basepairs(WLC_P__NT)
      real(dp) discretization, num_link_beads
      real(dp), dimension(2, 33) :: LL_dist
      real(dp) cumlinker, linker
      real(dp) urand(3)  ! random vector
      integer iter, i, j, k

      if (WLC_P__INCLUDE_NUC_TRANS) then
         if (WLC_P__INCLUDE_DISCRETIZE_LINKER) then
            ! figure out main discretization scheme
            discretization = WLC_P__LINKER_DISCRETIZATION
            num_link_beads = WLC_P__NBPL
            if (WLC_P__LINKER_TYPE == 'phased' ) then
               do i = 1, WLC_P__NP
                  ! initialize
                  iter = first_bead_of_chain(i)
                  ! set first linker later
                  iter = iter + num_link_beads
                  ! set middle linkers
                  do j = 2, WLC_P__NUM_NUCLEOSOMES ! hanging linker off nucleosomes
                     wlc_nucleosomeWrap(iter) = 147
                     wlc_basepairs(iter) = discretization
                     iter = iter + 1
                     if (iter + num_link_beads - 2 <= last_bead_of_chain(i) - num_link_beads) then
                        wlc_basepairs(iter:iter + num_link_beads - 2) = discretization
                        wlc_nucleosomeWrap(iter:iter + num_link_beads - 2) = 1
                        iter = iter + num_link_beads - 1
                     endif
                  enddo
                  ! first linker
                  wlc_nucleosomeWrap(first_bead_of_chain(i):first_bead_of_chain(i) + num_link_beads - 1) = 1
                  wlc_basepairs(first_bead_of_chain(i):first_bead_of_chain(i) + num_link_beads - 1) = discretization
                  ! last linker
                  wlc_basepairs(iter) = discretization
                  wlc_nucleosomeWrap(iter) = 147
                  iter = iter + 1
                  wlc_nucleosomeWrap(iter:iter + num_link_beads - 1) = 1
                  wlc_basepairs(iter:iter + num_link_beads - 1) = discretization
                  ! set last wlc_basepairs to 0 as reminder that this is not an actual extension
                  wlc_basepairs(last_bead_of_chain(i)) = 0
               enddo
            else if (WLC_P__LINKER_TYPE == 'voong') then
               ! read in the LL distribution
               open (UNIT = 5, FILE = "input/Voong_LL_23to55.txt", STATUS = "OLD")
               do i = 1, 33
                  read(5, *) LL_dist(:, i)
               enddo
               close (5)
               do k = 1, WLC_P__NP
                  linker = 0
                  do while (linker <= 23 .or. linker >= 55)
                     ! initialize
                     cumlinker = 0
                     iter = first_bead_of_chain(k)
                     ! set first linker later
                     iter = iter + num_link_beads
                     ! set middle linkers
                     do i = 2, WLC_P__NUM_NUCLEOSOMES ! hanging linker off nucleosomes
                        call random_number(urand)
                        do j = 1, 33
                           if (LL_dist(2, j) >= urand(1)) then
                              exit
                           endif
                        enddo
                        linker =  LL_dist(1, j)
                        discretization = linker/num_link_beads
                        print*, linker, discretization, num_link_beads
                        wlc_nucleosomeWrap(iter) = 147
                        wlc_basepairs(iter) = discretization
                        iter = iter + 1
                        if (iter + num_link_beads - 2 <= last_bead_of_chain(k) - num_link_beads) then
                           wlc_basepairs(iter:iter + num_link_beads - 2) = discretization
                           wlc_nucleosomeWrap(iter:iter + num_link_beads - 2) = 1
                           iter = iter + num_link_beads - 1
                        endif
                        cumlinker = cumlinker + linker
                     enddo
                     ! set last (and first) linker to conserve contour length
                     linker = (WLC_P__LL*(1 + WLC_P__NUM_NUCLEOSOMES) - cumlinker)/2
                     print*, 'attempted to sample from linker distribution for chain'
                  enddo
                  discretization = linker/num_link_beads
                  print*, linker, discretization, num_link_beads
                  ! first linker
                  wlc_nucleosomeWrap(first_bead_of_chain(k):first_bead_of_chain(k) + num_link_beads - 1) = 1
                  wlc_basepairs(first_bead_of_chain(k):first_bead_of_chain(k) + num_link_beads - 1) = discretization
                  ! last linker
                  wlc_basepairs(iter) = discretization
                  wlc_nucleosomeWrap(iter) = 147
                  iter = iter + 1
                  wlc_nucleosomeWrap(iter:iter + num_link_beads - 1) = 1
                  wlc_basepairs(iter:iter + num_link_beads - 1) = discretization
                  ! set last wlc_basepairs to 0 as reminder that this is not an actual extension
                  wlc_basepairs(last_bead_of_chain(k)) = 0
               enddo
            else
               print *, 'not recognized linker length type'
               stop
            endif
            ! testing sterics here !
            if (WLC_P__GJK_STERICS) then
               do iter = 1, 10000 ! check to make sure GJK checkng is not stochastic
                  call sameShapeTest()
                  call noIntersectX()
                  call intersectX()
                  call tangentX()
                  call noIntersectY()
                  call intersectY()
                  call tangentY()
                  call noIntersectZ()
                  call intersectZ()
                  call tangentZ()
                  call runtimeTest1()
                  call runtimeTest2()
                  call runtimeTest3()
                  call runtimeTest4()
                  call runtimeTest5()
                  call runtimeTest6()
               enddo
               print *, "SUCCESS: successful completion of all GJK collision unit tests"
            endif
         endif
      else
         wlc_nucleosomeWrap = 147
         wlc_basepairs = WLC_P__LL
      endif
      print*, wlc_basepairs
      print*, wlc_nucleosomeWrap

   end subroutine

end module nucleosome

max_stars_repo_name: SpakowitzLab/BasicWLC

max_issues_repo_name: SpakowitzLab/BasicWLC

max_forks_repo_name: SpakowitzLab/BasicWLC


text: ${SV_BFMS}/generic_rom/generic_rom_bfm.sv


text:
C***********************************************************************
C***********************************************************************
C
C   Version:        0.3
C   Last modified:  January 12, 1995
C   Authors:        Esmond G. Ng and Barry W. Peyton
C
C   Mathematical Sciences Section, Oak Ridge National Laboratory
C
C***********************************************************************
C***********************************************************************
C******     BTREE2 ..... BINARY TREE REPRESENTATION OF ETREE     *******
C***********************************************************************
C***********************************************************************
C
C   PURPOSE:
C       TO DETERMINE A BINARY TREE REPRESENTATION OF THE ELIMINATION
C       TREE, FOR WHICH EVERY "LAST CHILD" HAS THE MAXIMUM POSSIBLE
C       COLUMN NONZERO COUNT IN THE FACTOR.  THE RETURNED REPRESENTATION
C       WILL BE GIVEN BY THE FIRST-SON AND BROTHER VECTORS.  THE ROOT OF
C       THE BINARY TREE IS ALWAYS NEQNS.
C
C   INPUT PARAMETERS:
C       NEQNS           -   NUMBER OF EQUATIONS.
C       PARENT          -   THE PARENT VECTOR OF THE ELIMINATION TREE.
C                           IT IS ASSUMED THAT PARENT(I) > I EXCEPT OF
C                           THE ROOTS.
C       COLCNT          -   COLUMN NONZERO COUNTS OF THE FACTOR.
C
C   OUTPUT PARAMETERS:
C       FSON            -   THE FIRST SON VECTOR.
C       BROTHR          -   THE BROTHER VECTOR.
C
C   WORKING PARAMETERS:
C       LSON            -   LAST SON VECTOR.
C
C***********************************************************************
C
      SUBROUTINE  BTREE2 (  NEQNS , PARENT, COLCNT, FSON  , BROTHR,
     &                      LSON    )
C
C***********************************************************************
C
        INTEGER             BROTHR(*)     , COLCNT(*)     ,
     &                      FSON(*)       , LSON(*)       ,
     &                      PARENT(*)
C
        INTEGER             NEQNS
C
C***********************************************************************
C
        INTEGER*4           LROOT , NODE  , NDLSON, NDPAR
C
C***********************************************************************
C
        IF  ( NEQNS .LE. 0 )  RETURN
C
        DO  100  NODE = 1, NEQNS
            FSON(NODE) = 0
            BROTHR(NODE) = 0
            LSON(NODE) = 0
  100   CONTINUE
        LROOT = NEQNS
C       ------------------------------------------------------------
C       FOR EACH NODE := NEQNS-1 STEP -1 DOWNTO 1, DO THE FOLLOWING.
C       ------------------------------------------------------------
        IF  ( NEQNS .LE. 1 )  RETURN
        DO  300  NODE = NEQNS-1, 1, -1
            NDPAR = PARENT(NODE)
            IF  ( NDPAR .LE. 0  .OR.  NDPAR .EQ. NODE )  THEN
C               -------------------------------------------------
C               NODE HAS NO PARENT.  GIVEN STRUCTURE IS A FOREST.
C               SET NODE TO BE ONE OF THE ROOTS OF THE TREES.
C               -------------------------------------------------
                BROTHR(LROOT) = NODE
                LROOT = NODE
            ELSE
C               -------------------------------------------
C               OTHERWISE, BECOMES FIRST SON OF ITS PARENT.
C               -------------------------------------------
                NDLSON = LSON(NDPAR)
                IF  ( NDLSON .NE. 0 )  THEN
                    IF  ( COLCNT(NODE) .GE. COLCNT(NDLSON) )  THEN
                        BROTHR(NODE) = FSON(NDPAR)
                        FSON(NDPAR) = NODE
                    ELSE
                        BROTHR(NDLSON) = NODE
                        LSON(NDPAR) = NODE
                    ENDIF
                ELSE
                    FSON(NDPAR) = NODE
                    LSON(NDPAR) = NODE
                ENDIF
            ENDIF
  300   CONTINUE
        BROTHR(LROOT) = 0
C
        RETURN
      END

max_stars_repo_name: Pangqiyuangh/SeIInv

max_issues_repo_name: Pangqiyuangh/SeIInv

max_forks_repo_name: Pangqiyuangh/SeIInv


text:
module git2_remote
    use, intrinsic :: iso_c_binding
    use :: git2_proxy
    implicit none
    private

    ! git_fetch_prune_t
    enum, bind(c)
        enumerator :: GIT_FETCH_PRUNE_UNSPECIFIED
        enumerator :: GIT_FETCH_PRUNE
        enumerator :: GIT_FETCH_NO_PRUNE
    end enum

    ! git_remote_autotag_option_t
    enum, bind(c)
        enumerator :: GIT_REMOTE_DOWNLOAD_TAGS_UNSPECIFIED = 0
        enumerator :: GIT_REMOTE_DOWNLOAD_TAGS_AUTO
        enumerator :: GIT_REMOTE_DOWNLOAD_TAGS_NONE
        enumerator :: GIT_REMOTE_DOWNLOAD_TAGS_ALL
    end enum

    integer(kind=c_int), parameter, public :: GIT_FETCH_OPTIONS_VERSION    = 1
    integer(kind=c_int), parameter, public :: GIT_REMOTE_CALLBACKS_VERSION = 1
    integer(kind=c_int), parameter, public :: GIT_REMOTE_CALLBACKS_INIT(1) = [ GIT_REMOTE_CALLBACKS_VERSION ]

    integer(kind=c_int), parameter, public :: GIT_FETCH_OPTIONS_INIT(6) = [ GIT_FETCH_OPTIONS_VERSION, &
                                                                            GIT_REMOTE_CALLBACKS_INIT, &
                                                                            GIT_FETCH_PRUNE_UNSPECIFIED, &
                                                                            1, &
                                                                            GIT_REMOTE_DOWNLOAD_TAGS_UNSPECIFIED, &
                                                                            GIT_PROXY_OPTIONS_INIT ]
end module git2_remote

max_stars_repo_name: interkosmos/fortran-git

max_issues_repo_name: interkosmos/fortran-git

max_forks_repo_name: interkosmos/fortran-git


text:
subroutine diffuse(x1, y1, z1, ib, jb, kb, dt,imt,jmt,km,dxv,dyu,dzt,h,ah,av,do3d,doturb)

!============================================================================
! Add a small displacement to a particle s.t. it is still in the model area.
!
!  Input:
!    dt             : time step of trajectory iteration in seconds
!    imt,jmt,km     : grid index sizing constants in (x,y,z), are for
!                     horizontal and vertical rho grid [scalar]
!    dxv            : Horizontal grid cell walls areas in x direction [imt,jmt-1]
!    dyu            : Horizontal grid cell walls areas in y direction [imt-1,jmt]
!    dzt            : Height of k-cells in 3 dim in meters on rho vertical grid. [imt,jmt,km]
!    h              : Depths [imt,jmt]
!    ah             : horizontal diffusion in m^2/s. Input parameter.
!                   : For -turb and -diffusion flags
!    av             : vertical diffusion in m^2/s. Input parameter. For -diffusion flag.
!    do3d           : Flag to set whether to use 3d velocities or not
!    doturb         : Flag to set whether or not to use turb/diff and which kind if so
!
!  Input/Output:
!    x1, y1, z1     : Current positions of the particle. Are updated by the subroutine.
!    ib, jb, kb     : Indices for current box for the particle. Are updated by the subroutine.
!
!  Other parameters used in function:
!    tmpX, tmpY, tmpZ   : Temporal position
!    xd, yd, zd         : Displacement
!    itno               : Number of iterations
!    tmpi, tmpj, tmpk   : Temporal box indices
!    tryAgain           : Tells whether to continue displace
!============================================================================

implicit none

integer,        intent(in)                                  :: do3d, doturb
real(kind=8),   intent(in)                                  :: ah,av,dt
real*8,         intent(in),     dimension(imt-1,jmt)        :: dyu
real*8,         intent(in),     dimension(imt,jmt-1)        :: dxv
real*8,         intent(in),     dimension(imt,jmt)          :: h
real(kind=8),   intent(in),     dimension(imt,jmt,km,2)     :: dzt
integer,        intent(in)                                  :: imt,jmt,km
real(kind=8),   intent(in out)                              :: x1, y1, z1
integer,        intent(in out)                              :: ib,jb,kb
integer                                                     :: itno,tmpi, tmpj, tmpk
real(kind=8)                                                :: xd, yd, zd, tmpX, tmpY, tmpZ
logical                                                     :: tryAgain

! print *,''
! print *,'start of diffusion'
! print *,'x1=',x1,' y1=',y1,' ib=',ib,' jb=',jb

tryAgain = .FALSE.

! Is particle within model area?
! KMT: I do not understand how the vertical checks are correct here, so I am changing them

!   if(ib>=1 .AND. ib<=imt .AND. jb>=1 .AND. jb<=jmt .AND. km+1-kmt(ib,jb)<=kb .AND. kb>=1 ) then
if(ib>1 .AND. ib<imt .AND. jb>1 .AND. jb<jmt .AND. km>=kb .AND. kb>=1 ) then
    tryAgain = .TRUE.
else
!     print *,'outside model domain in diffusion',ib,jb,km,kb
! !     print *,'outside model domain in diffusion',ib,jb,km+1-kmt(ib,jb),kb
end if

itno=0
do while(tryAgain)
    itno=itno+1
    ! find random displacement
    call displacement(xd, yd, zd, ib, jb, dt,imt,jmt,h,ah,av,dxv,dyu,do3d,doturb)
!     CALL displacement(xd, yd, zd, ib, jb, kb, dt,imt,jmt,kmt,h,ah,av)
!     print *,'xd=',xd,' yd=',yd,' zd=',zd
!     print *,'dzt(ib,jb,:,1)=',dzt(ib,jb,:,1)
    ! Convert displacement from meters to model coordinates
    ! KMT note: following example in master code for orc project
    xd = xd/dxv(ib,jb)
    yd = yd/dyu(ib,jb)
    zd = zd/dzt(ib,jb,kb,1) !KMT: this should be better than dz()
!     print *,'dzt(ib,jb,kb,1)=',dzt(ib,jb,kb,1)

    ! TO DO: improve so that dzt is calculated including free surface and for correct time
!     zd = zd/dz(kb) !should be replaced for bottom box and include ssh (original note)

!     print *,'xd=',xd,' yd=',yd,' zd=',zd

    ! Update position temporarily
    tmpX = x1 + xd
    tmpY = y1 + yd
    tmpZ = z1 + zd
    ! Update box number temporarily
    tmpi = int(tmpX) + 1
    tmpj = int(tmpY) + 1
    tmpk = int(tmpZ) + 1
!     print *,'temp dzt(ib,jb,kb,1)=',dzt(tmpi,tmpj,tmpk,1)

!     ! Check if particle is on an open boundary. Using rco example.
!     if(tmpi==1 .AND. tmpj>=1 .AND. tmpj<=jmt .AND. km+1-kmt(tmpi,tmpj)<=tmpk .AND. tmpk>=1 ) then
!         tryAgain = .FALSE.
!     end if

    ! check that column is deep enough
    ! KMT: this is the same check as above that looks wrong and I changed
!     if( 1<=tmpi .AND. tmpi<=imt .AND. 1<=tmpj .AND. tmpj<=jmt .AND. km+1-kmt(tmpi,tmpj)<=tmpk .AND. tmpk>=1 ) then
    if(tmpi>1 .AND. tmpi<imt .AND. tmpj>1 .AND. tmpj<jmt .AND. km>=tmpk .AND. tmpk>=1 ) then
!         print *,'km=',km,' km+1-kmt(tmpi,tmpj)=',km+1-kmt(tmpi,tmpj),' tmpk=',tmpk
        tryAgain = .FALSE.
!         print *,'found new position for drifter that is within model domain, so move on.'
        ! if false then a new position for the particle has been found and we exit the loop
    end if

    ! KMT: Need to check to see if new position is on masked land.
    ! Can do this by checking dzt at the new position, which is zero when on land.
    ! If we are on land, keep tryAgain = .TRUE. so that this displacement choice is
    ! not used. This can alter the value of tryAgain just given above.
    if(dzt(tmpi,tmpj,tmpk,1)==0.) then
        tryAgain = .TRUE.
    endif

    ! If tryAgain is still true, the particle is outside model area. The
    ! displacement is not saved, but we make a new try to displace.

    ! "Infinite loop?"
    if(itno>=100000 .AND. tryAgain) then
        tryAgain = .FALSE.
        write(*,*)"Particle stuck in infinite diffusion loop. No diffusion added.",ib,jb,kb
        stop 34956
        tmpX=x1 ; tmpY=y1 ; tmpZ=z1
        tmpi=ib ; tmpj=jb ; tmpk=kb
    end if

  enddo

! Update return position
! if(abs(x1-tmpX).gt.5. .and. abs(x1-tmpX).lt.1000.) print *,x1,tmpX,xd,ib,jb,dxv(ib,jb)
x1 = tmpX
y1 = tmpY
ib = tmpi
jb = tmpj

! print *,'After diffusion'
! print *,'x1=',x1,' y1=',y1,' ib=',ib,' jb=',jb

if(do3d==1) then
    z1 = tmpZ
    kb = tmpk
endif

if(x1.ne.dble(idint(x1))) ib=idint(x1)+1 ! make sure the index corresponds to the right grid cell

!print *,'slut',itno,kmt(ib,jb),ib,jb,kb,x1,y1,z1

end subroutine diffuse

!===============================================================================
! Calculate a random displacement
! (sqrt(-4ah*dt*log(1-q1))*cos(2PIq2),
!  sqrt(-4ah*dt*log(1-q1))*sin(2PIq2),
!  sqrt(-4av*dt*log(1-q3))*cos(2PIq4))
! where av and ah are set in run.in, dt is the model time step and q1,q2,q3,q4 are random
! numbers between 0 and 1.
!
! Arguments :
! xd, yd, zd : Variables in which the displacement will be stored
! dt: Model time step
!===============================================================================
subroutine displacement(xd, yd, zd, ib, jb, dt,imt,jmt,h,ah,av,dxv,dyu,do3d,doturb)

implicit none

integer,        intent(in)                          :: ib,jb,imt,jmt   ! box indices
integer,        intent(in)                          :: do3d, doturb
real(kind=8),   intent(in)                          :: dt,ah,av
real*8,         intent(in), dimension(imt,jmt)      :: h
real*8,         intent(in), dimension(imt-1,jmt)    :: dyu
real*8,         intent(in), dimension(imt,jmt-1)    :: dxv
real(kind=8),   intent(out)                         :: xd, yd, zd
real(kind=8),   parameter                           :: pi = 3.14159265358979323846
real(kind=8)                                        :: q1, q2, q3, q4, r, elip, xx, yy
real(kind=8)                                        :: grdx, grdy, grad, theta
integer                                             :: ip,im,jp,jm

! random generated numbers between 0 and 1
q1 = rand()
q2 = rand()
q3 = rand()
q4 = rand()

! Horizontal displacements in meters
R = sqrt(-4*ah*dt*log(1-q1))
xd = R * cos(2*PI*q2)
yd = R * sin(2*PI*q2)

! Vertical displacement in meters
if(do3d==1) then
    R = sqrt(-4*av*dt*log(1-q3))
    zd = R*cos(2*PI*q4)
else
    zd = 0.
endif

! print *,' '
! print *,'circle'
! print *,'xd=',xd,' yd=',yd

if(doturb==3) then
    !_______________________________________________________________________________
    ! The diffusion is here set on an ellipse instead of a circle
    ! so that the diffusion is higher along the isobaths and weaker
    ! in the perpendicular direction (up/downhill)

    ! for index in i direction for one cell positive from the current cell
    ip=ib+1
    ! KMT: Is this a periodic BC? Replace with just the edge value so it won't
    ! elliptical but just at the edge.
    ! I think the check should happen for ip==imt since in the x direction for u
    ! the grid only goes up to imt-1 for the arakawa c grid.
    if(ip.eq.imt) ip = imt-1 ! use ib instead of ib+1, and here ib=imt-1
!     if(ip.eq.imt+1) ip=1

    ! for index in i direction for one cell negative from the current cell
    im=ib-1
    ! KMT: Is this a periodic BC? Replace with just the edge value so it won't
    ! elliptical but just at the edge.
    if(im.eq.0) im = 1 ! use ib instead of ib-1, and here ib=1
!     if(im.eq.0) im=imt

     ! for index in j direction for one cell positive from the current cell
    jp=jb+1
    ! KMT: In v they don't use period BC.
    ! I think the check should happen for jp==jmt since in the y direction for v
    ! the grid only goes up to jmt-1 for the arakawa c grid.
    if(jp.eq.jmt) jp = jmt-1 ! use jb instead of jb+1, and here jb=jmt-1
!     if(jp.eq.jmt+1) jp=jmt

    ! for index in j direction for one cell negative from the current cell
    jm=jb-1
    ! KMT: In v they don't use period BC and this check is fine as is.
    if(jm.eq.0) jm = 1

    ! just depth gradient
    ! grdx=float(kmt(ip,jb)-kmt(im,jb)) ! zonal      depth gradient (zw should be used)
    ! grdy=float(kmt(ib,jp)-kmt(ib,jm)) ! meridional depth gradient
    ! KMT: switching the bathymetry/depths array and grid distances in here
    ! since user manual (pg 21) shows grdx=\Delta h/\Delta x, for example
!     grdx=(h(ip,jb)-h(im,jb))/(0.5*(dxv(im,jb)+dxv(ib,jb))+0.5*(dxv(ib,jb)+dxv(ip,jb))) ! zonal depth gradient
    ! Same as easier-to-read previous line but more efficient:
    grdx=(h(ip,jb)-h(im,jb))/(0.5*dxv(im,jb) + 0.5*dxv(ip,jb) + dxv(ib,jb)) ! zonal depth gradient

!     grdy=(h(ib,jp)-h(ib,jm))/(0.5*(dyu(ib,jm)+dyu(ib,jb)) + 0.5*(dyu(ib,jb)+dyu(ib,jp))) ! meridional depth gradient
    ! Same as easier-to-read previous line but more efficient:
    grdy=(h(ib,jp)-h(ib,jm))/(0.5*dyu(ib,jm) + dyu(ib,jb) + 0.5*dyu(ib,jp)) ! meridional depth gradient

    grad=dsqrt(grdx**2+grdy**2)       ! total depth gradient

    ! angle between the eastward direction and
    ! the constant depth direction (=0 if only meridional slope)
    if(grad.eq.0.d0) then
        ! theta=0.d0 + pi/2.d0
        theta=0.d0
    else
        ! theta=dasin(grdx/grad) +pi/2.d0 ! varför 90 grader mer? ny2
        theta=dasin(grdx/grad) ! ny1
    endif

    ! elliptic horizontal distribution of the diffusion
    !elip=dabs(1.d2*grad)+1.d0 ! för 1/h
    !elip=dabs(1.d-2*grad)+1.d0 ! för h
    elip=dabs(grad)+1.d0 ! för h

    !elip=amin1(10.d0,grad) ! gives slightly higher relative dispersion
    elip=amin1(5.d0,elip) !

    !print *,'elip=',elip,grad,theta

    !The circular disk is stretched into an elliptic disk
    ! with unchanges surface area by
    xx=xd*elip
    yy=yd/elip

    ! coordinate transformation to put the diffusion on an
    ! ellipse with the maxium diffusion along the isobaths
    xd= xx*dcos(theta)-yy*dsin(theta)
    yd=-xx*dsin(theta)+yy*dcos(theta)

!     print *,'ellipse'
!     print *,'xd=',xd,' yd=',yd
    ! print *,'theta=',theta,' xx=',xx,' yy=',yy,' grdx=',grdx,' grdy=',grdy

    !print *,'theta=',theta*180./pi,elip,grdx/grad

    !if(jb.gt.400) stop 3096
    !_______________________________________________________________________________
endif

end subroutine displacement

max_stars_repo_name: dcherian/tracpy

max_issues_repo_name: dcherian/tracpy

max_forks_repo_name: dcherian/tracpy


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program 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.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      SUBROUTINE GETDIM ( FNAME  , ITYPE  , DIM    , IPRDEP , ITMDEP ,
     *                             LOCDEP , NDIM   , IERROR , OPTION )
!
!
!     Deltares        MARINE & COASTAL MANAGEMENT
!
!     CREATED            : May '96  by L. Postma
!
!     MODIFIED           :
!
!     FUNCTION           : ODS GETDIM routine for DELWAQ HIS-files
!
!     SUBROUTINES CALLED :
!
!     LOGICAL UNITS      :
!
!     PARAMETERS    :
!
!     NAME    KIND     LENGTH     FUNCT.  DESCRIPTION
!     ---------------------------------------------------------
!     FNAME   CHAR*256   3        IN/LOC  Complete file name
!     ITYPE   INTEGER    1        INPUT   File type
!     DIM     CHAR*3     1        INPUT   Wanted dimension
!     IPRDEP  INTEGER    1        INPUT   Par code for dimensions
!     ITMDEP  INTEGER    1        INPUT   Time code for dimensions
!     LOCDEP  INTEGER    1        INPUT   Loc code for dimensions
!     NDIM    INTEGER    5        OUTPUT  Wanted dimensions
!     IERROR  INTEGER    1        OUTPUT  Error code
!     OPTION  CHAR*256   1        IN/OUT  For future use
!
!     NOTE1: FNAME(3) is used as local character space
!     NOTE2: NDIM is NOT according to ODS specs, it returns:
!            NDIM(1) = nr of substances in the file
!            NDIM(2) = nr of locations  in the file
!            NDIM(3) = nr of time steps in the file
!
!
      CHARACTER*256 FNAME(3) , OPTION
      CHARACTER*3   DIM
      DIMENSION     NDIM(5)
      character*256         :: ext     ! file extension
      integer               :: extpos  ! position of extension
      integer               :: extlen  ! length of file extension
      logical               :: mapfil  ! true if map file extension

!
!         Open the DELWAQ .HIS file
!
      CALL DHOPNF ( 10 , FNAME(1) , 24 , 2 , IERROR )
      IF ( IERROR .NE. 0 ) RETURN

      ! map or his

      call dhfext(fname(1), ext, extpos, extlen)
      call dhucas(ext, ext, extlen)
      if ( ext .eq. 'MAP' ) then
         mapfil = .true.
      else
         mapfil = .false.
      endif
!
!         Read primary system characteristics
!
      READ ( 10 , ERR=100 )   FNAME(3)(1:160)
      READ ( 10 , ERR=110 )   NOTOT, NODUMP
      READ ( 10 , ERR=120 ) ( FNAME(3)(181:200) , K = 1,NOTOT )
      if ( .not. mapfil ) then
         READ ( 10 , ERR=130 ) ( IDUMMY, FNAME(3)(221:240) , K = 1,NODUMP )
      endif
!
!         Read the values at all times
!
      NTT   = NODUMP*NOTOT
      NOTIM = 0
cjvb  NTT = 0
   10 READ ( 10 , ERR=140 , END=20 )   IDUMMY, ( ADUMMY , K=1,NTT )
      NOTIM = NOTIM + 1
      GOTO 10
!
!         Supply the desired statistics
!
   20 NDIM(1) = NOTOT
      NDIM(2) = NODUMP
      NDIM(3) = NOTIM
      GOTO 200
!
!         Supply the desired statistics
!
  100 IERROR = 10
      GOTO 200
  110 IERROR = 11
      GOTO 200
  120 IERROR = 12
      GOTO 200
  130 IERROR = 13
      GOTO 200
  140 IERROR = 14
!
  200 CLOSE ( 10 )
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
! { dg-do compile }
! PR 34549 - a real value was accepted for shift.
program main
  implicit none
  real, dimension(2,2) :: r
  data r /1.0, 2.0, 3.0, 4.0/
  print *,cshift(r,shift=2.3,dim=1) ! { dg-error "must be INTEGER" }
end program main

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
        integer function sigprod(n)
        integer n, i
        double precision s, w
        s = 1.0
        w = 1.0
        do 10, i=1, n
           w = w / i
           s = s + w
   10   continue
        sigprod = s
        end

        double precision function ff2(n)
        integer n, i, fact
        double precision sum1
        sum1 = 0.0
        fact = 1.0
        do 60, i=1, n
           fact = fact * i
           sum1 = sum1 + 1/fact
   60   continue
        ff2 = sum1
        end

        integer function sigsum(n)
        integer n, i
        double precision s, w
        s = 1.0
        w = 1.0
        do 20, i=1, n
           w = w * i
           s = s + w
   20   continue
        sigsum = s
        end

        integer function sigsum1(n)
        integer n, i, s
        s=n
        do 30, i=n-1, 1, -1
           s = i * (1 + s)
   30   continue
        sigsum = s+1
        end

        integer function chebyshev(n,x)
        integer n, x, i, u0, u1, v
        u0 = 1
        u1 = x
        do 50, i=2, n
           v = u1
           u1 = -u0 + 2*x*u1
           u0 = v
   50   continue
        chebyshev = u1
        end

        double precision function bessel(z, nu, m)
        integer i, z, nu, m
        double precision res, u
        res = 0
        u = 1
        do 40, i=0, m-1
           res = res + u
           u = -u*z**2 / (4*(i+nu+1)*(i+1))
   40   continue
        bessel = res
        end

        double precision function mystery(z, nu, m)
        integer z, nu, m, i
        double precision res, t
        res = 0
        t = 1
        do 70, i=0, m
           res = res + t
           t = -t*z**2 / (4*(i+nu+1)*(i+1))
   70   continue
        mystery = res
        end

max_stars_repo_name: dacmot/RevEngTools

max_issues_repo_name: dacmot/RevEngTools

max_forks_repo_name: dacmot/RevEngTools


text:
      double precision hmass,hwidth
      common/hmass/hmass,hwidth

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


MODULE mymodule

 IMPLICIT NONE
 CONTAINS

  FUNCTION diff(param1, param2)
    INTEGER, INTENT(IN) :: param1
    INTEGER, INTENT(IN) :: param2
    INTEGER :: diff
    diff = param1 - param2
  END FUNCTION diff

END MODULE mymodule

PROGRAM driver
  USE mymodule
  IMPLICIT NONE
  INTEGER :: a, b

  b = 5678
  a = diff(param2 = b, param1 = 1234)

  PRINT *, 'A is ', a
END PROGRAM driver


text:
      subroutine far(y,c,s,r,in,r2,sig,n,p,po,pv,ier,io)
c********************************************
c  Multi level iterative AR fit which loops back to discard lags
c   whose coefficients have prob value of abs(t) statistics < pv
c  Input: y - n observations, p - order in, pv - two sided p-value
c  Output: c - po AR coefficients, s - standard errors of c coeff's
c    r(1),...,r(n) - residuals
c    r2 - R square,  sig - standard deviation of residuals
c    in - array of lag values for significant lags, po - order out
c  Version 3-12-2005
c   ier = 1: A is not positive definite, ier = 2: p < 1, ier = 3: p > n/2
c============================================
      allocatable::a,b,d,w
      integer p,pi,po
      integer in(p)
      real y(n),r(n),s(p),pv
      real*8 c(p),a(:),b(:),w(:),d(:),sy,ssy
      intent(in) y,n,p,pv
      intent(out) c,s,r,in,po
c****************************
c  Trap p
      if(p<1) then
       write(io,*) 'p in AR(p) ',p,' < 1'
       ier=2
       return
        elseif(p>n/2) then
         write(io,*) 'p in AR(p) ',p,' > n/2 ',n/2,' n = ',n
         ier=3
         return
      endif
c-------
      allocate(a(0:p*p-1),b(0:p),d(0:p),w(0:p),stat=ibad)
      if(ibad/=0) then
       write(io,*) 'Unable to allocate work space for arrays in far'
       return
      endif
c  Covariances in a
      call cr(y,a,b,w,sy,ssy,n,p,io)
      do k=1,p
       c(k)=0.d0
       in(k)=k
      enddo
      pi=p
      po=p-1
      dowhile(po>0)
c  Iteration on ols
       call ar(y,c,s,r,a,b,d,w,n,p,sy,ssy,r2,sig,pv,ier,in,pi,po,io)
       if(po<pi) then
        pi=po
         else
          exit
       endif
      enddo
c-------
      deallocate(a,b,d,w)
      return
      end
c
c============================================
c
      subroutine ar(y,c,s,r,a,b,d,w,n,p,sy,ssy,r2,sig,pv,ier,in,pi,po,io
     *)
c********************************************
c  Computes estimates of parameters of an AR(p) model and the
c   residuals of the OLS fit.
c  Input: y - data array of n observations, order p
c  Output: r - array of e(k)=y(k)-c1*y(k-1)-c2*y(k-1) ...
c   c - p AR coefficients, s - standard errors of coefficients
c   sig- mean and standard deviation of residuals, r2 - R squared
c   po - Number of p values of t statistics < pv
c============================================
      integer p,pi,po,t,n
      integer in(p)
      real y(n),r(n),s(p),sig,pv
      real*8 a(p*p),b(p),d(p),w(p),c(p),sy,ssy,sse,sm,hd,dy,dnp
      intent(in) y,n,p,pi,pv,sy,ssy
      intent(out) c,s,r,sig,r2,po
      intent(in out) in
c********************************************
c  Cholesky decomposition of a
      call dcholdc(a,pi,d,ier,io)
      if(ier/=0) then
       write(io,*) 'DCHOLDC failed, A is not positive definite'
       write(io,*) 'pi = ',pi
       return
      endif
c-------
      call dcholsl(a,pi,d,b,c)
c  Compute residuals, mean, and standard deviation
      sse=0.d0
      hd=0.
      do t=in(pi)+1,n
       sm=0.d0
       do k=1,pi
        sm=sm+c(k)*(dble(y(t-in(k)))-sy)
       enddo
       dy=dble(y(t))-sy-sm
       r(t)=sngl(dy)
       hd=hd+dy
       sse=sse+dy*dy
      enddo
c  Zero out first pi residuals
      do t=1,in(pi)
       r(t)=0.
      enddo
      dnp=dfloat(n-pi)
      hd=hd/dnp
      sse=sse-dnp*hd*hd
c  R squared
      r2=1.-sngl(sse/ssy)
      r2=r2-float((pi-1)/(n-pi))*(1-r2)
      sig=sngl(dsqrt(sse/dfloat(n-pi-1)))
c-------
c  Inverse of A
      call dinverse(a,pi,d)
c  Standard errors of c's
      prb=1.-pv
      po=0
      do k=1,pi
       s(k)=sig*sngl(dsqrt(a((k-1)*pi+k)))
       ch=sngl(c(k))/s(k)
       call cdchi(ch*ch,1.,pc,io)
       if(pc>prb) then
        po=po+1
        in(po)=in(k)
       endif
      enddo
c-------
      if(po<pi .and. po>0) then
c  Reduce down the covariance matrix
       call reduce(a,b,w,ssy,p,po,in,io)
      endif
c-------
      return
      end
c
c===========================================
c
      subroutine cr(y,a,b,w,sy,ssy,n,p,io)
c********************************************
c  Covariances of y in a & variances in b,w
c  sy - Mean of y, ssy - Sum of squares of y
c============================================
      integer p,t
      real y(n)
      real*8 a(p*p),b(p),w(p),sy,ssy,ac
      intent(in) y,n,p
      intent(out) a,b,w,sy,ssy
c********************************************
c  Compute mean
      sy=0.d0
      do t=1,n
       sy=sy+dble(y(t))
      enddo
      sy=sy/dfloat(n)
c-------
c  Compute sums
      ssy=(dble(y(n))-sy)*(dble(y(n))-sy)
      do k=1,p
       ac=0.d0
c  Compute covariances
       do t=1,n-k
        ac=ac+(dble(y(t))-sy)*(dble(y(t+k))-sy)
        if(k==1) then
         ssy=ssy+(dble(y(t))-sy)*(dble(y(t))-sy)
        endif
       enddo
       b(k)=ac
       w(k)=ac
      enddo
c-------
c  Put in upper triangle of matrix a
      a(p*p)=ssy
      do j=1,p-1
       do k=j+1,p
        a((k-1)*p+j)=b(k-j)
       enddo
       a((j-1)*p+j)=ssy
      enddo
c-------
      return
      end
c
c===========================================
c
      subroutine reduce(a,b,w,ssy,p,po,in,io)
c********************************************
c  Rearrange a & b for reduced least squares fit
c============================================
      integer p,po
      integer in(p)
      real*8 a(p*p),b(p),w(p),ssy
      intent(in) w,ssy,in,p,po
      intent(out) a,b
c********************************************
      do k1=1,po
       b(k1)=w(in(k1))
       a((k1-1)*po+k1)=ssy
       if(po>0) then
        do k2=k1+1,po
         a((k2-1)*po+k1)=w(in(k2)-in(k1))
        enddo
       endif
      enddo
c-------
      return
      end

max_stars_repo_name: emammendes/bispectrum

max_issues_repo_name: emammendes/bispectrum

max_forks_repo_name: emammendes/bispectrum


text:
asset: %buffer
    %buffer svar buffer.data
    %buffer svar buffer.size

: recreate-buffer  ( buffer - )
    >r  r@ buffer.size @ allocate throw  r> buffer.data ! ;

: unload-buffer  ( buffer - )
    buffer.data @ free throw ;

: init-buffer  ( size buffer - )
    >r  dup r@ buffer.size !  allocate throw  r@ buffer.data !
    ['] recreate-buffer ['] unload-buffer r> register ;

: buffer  ( size - )
    %buffer *struct  init-buffer ;

: buffer:   ( size - <name> )
    create  buffer  does>  buffer.data @ ;

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine


text:
subroutine dump_pickup()
    use global, only: casename,tt,DumpClock,tsg,xyz
    implicit none
    character(len=64) :: fn
    write(fn,"(I10.10)") int(tt/DumpClock)+1
    print*, fn
    open(0,file=trim(casename)//'.pickup.'//trim(fn)//'.data',form='unformatted',&
        access='stream',convert='BIG_ENDIAN',status='replace')
    write(0) xyz,tsg,tt
    close(0)
end subroutine dump_pickup

max_stars_repo_name: jinbow/Octupus

max_issues_repo_name: jinbow/Octupus

max_forks_repo_name: jinbow/Octupus


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program 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.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      INTEGER FUNCTION SETDLM(DLMCHR, NDLM, IGNCHR, NIGN)
!
      IMPLICIT INTEGER (A-Z)
      CHARACTER*1 DLMCHR(1), IGNCHR(1)
      COMMON /ZDLMTZ/DLM(256), IGN(256), IDLM, IIGN, PUSH(256), PTR
      DATA CHR/0/
!
      IF (NDLM.LE.0 .OR. NDLM.GT.256) GOTO 1001
      IF (NIGN.LE.0 .OR. NIGN.GT.256) GOTO 1002
      DO 10 I=1,256
      DLM(I) = 0
10    IGN(I) = 0
      IDLM = 0
      IIGN = 0
      DO 20 I=1,NDLM
! *** Convert CHARACTER to INTEGER.
      CHR = ICHAR(DLMCHR(I))
! *** Conversion completed.
      IF (DLM(CHR+1).EQ.0) IDLM = IDLM + 1
20    DLM(CHR+1) = 1
      DO 30 I=1,NIGN
! *** Convert CHARACTER to INTEGER.
      CHR = ICHAR(IGNCHR(I))
! *** Conversion completed.
      IF (IGN(CHR+1).EQ.0) IIGN = IIGN + 1
30    IGN(CHR+1) = 1
      SETDLM = 0
      RETURN
1001  SETDLM = 1
      RETURN
1002  SETDLM = 2
      RETURN
      END