text:
	SUBROUTINE RS_GGRP ( ob, here, endob, numb, shift, iret )
C************************************************************************
C* RS_GGRP								*
C*									*
C* This subroutine finds a group of numbers in OB, updates the index	*
C* pointer to OB, and returns the difference between the new position	*
C* and the old position.  A number is defined as any digit between 0 	*
C* and 9, or a slash (/).  The OB must be a continuous string of digits	*
C* and spaces.  The numbers need not be surrounded by spaces.		*
C*									*
C* RS_GGRP  ( OB, HERE, ENDOB, NUMB, SHIFT, IRET )			*
C*									*
C* Input parameters:							*
C*	OB		CHAR*		Encoded observation		*
C*	ENDOB		INTEGER		Length of OB			*
C*	NUMB		INTEGER		Number of consecutive integers	*
C*					 to search for			*
C*									*
C* Input and output parameters:						*
C*	HERE		INTEGER		Current position in OB		*
C*									*
C* Output parameters:							*
C*	SHIFT		INTEGER		Total change of HERE		*
C*	IRET		INTEGER		Return code			*
C*					 -1 = Failed to find group	*
C*					  0 = Successful		*
C*					  3 = Found '333' group instead	*
C*					  4 = Found '444' group instead *
C*					  5 = Found '555' group instead *
C**									*
C* Log:									*
C* J. Nielsen/MIT	10/86						*
C* J. Nielsen/MIT	 2/89						*
C* F. Cuq/UCLA		 8/90	Integer*2 -> Integer*4			*
C* J. Nielsen/TAMU	 2/92	Gempacized				*
C************************************************************************
	CHARACTER*(*)	ob
	INTEGER		endob, numb
C*
	INTEGER		here
C*
	INTEGER		shift, iret
C*
	INTEGER		here1
	LOGICAL		RS_NMER, done
C------------------------------------------------------------------------
	shift = 0
	iret = 0
	here1 = here
	done = .false.
C
C*	Loop to look for NUMB consectutive digits
C
	DO WHILE  ( .not. done )
C
C*	  Search for first non-blank character
C
	  DO WHILE  ( ob ( here:here ) .eq. ' ' )
C
C*	    Make sure there's enough space left for digits
C
	    IF  ( ( here + numb - 1 ) .gt. endob )  THEN
		iret = -1
		here = here1
		shift = 0
		RETURN
	    ENDIF
C
C*	    Try next character
C
	    here = here + 1
	    shift = shift + 1
	  END DO
C
C*	  Look for enough numbers and/or slashes
C
	  i = 0
	  DO WHILE  ( ( i .lt. numb ) .and.
     +		      ( ( RS_NMER ( ob ( here+i : here+i ), 1 ) ) .or.
     +		        ( ob ( here+i : here+i ) .eq. '/' ) ) )
	    i = i + 1
	  END DO
C
C*	  Find enough?
C
	  IF  ( i .eq. numb )   THEN
	    done = .true.
	  ELSE
C
C*	    Check for special 3-digit code group
C
	    IF  ( i .eq. 3 )  THEN
	      IF  ( ob ( here:here+2 ) .eq. '333' )  THEN
		iret = 3
		RETURN
	      ELSE IF  ( ob ( here:here+2 ) .eq. '444' )  THEN
		iret = 4
		RETURN
	      ELSE IF  ( ob ( here:here+2 ) .eq. '555' )  THEN
		iret = 5
		RETURN
	      ENDIF
	    ENDIF
C
C*	    Look for next space, to start over
C
	    here = here + i
	    shift = shift + i
	    DO WHILE  ( ob ( here:here ) .ne. ' ' )
C
C*	      Make sure there's enough space left for digits
C
	      IF  ( ( here + numb ) .gt. endob )  THEN
		iret = -1
		here = here1
		shift = 0
		RETURN
	      ENDIF
C
C*	      Try next character
C
	      here = here + 1
	      shift = shift + 1
	    END DO
	  ENDIF
C
C*	Go back to the beginning, if necessary
C
	END DO
C
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine clawpack46_setaux(maxmx,maxmy,mbc,mx,my,
     &      xlower,ylower,dx,dy,maux,aux)
      implicit none

      integer maxmx, maxmy,mbc,mx,my,maux
      double precision xlower,ylower,dx,dy
      double precision aux(1-mbc:maxmx+mbc,1-mbc:maxmy+mbc, maux)

      integer i,j
      double precision yj

      do i = 1-mbc,mx+mbc
         do j = 1-mbc,my+mbc
            yj = ylower + (j-0.5d0)*dy
            aux(i,j,1) = yj
         enddo
      enddo

      return
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


text:
      SUBROUTINE HOCAPD(
     &         CCENT1, CCENT2, POINT, DIST2, CUVEC, CLEN)
      IMPLICIT NONE
      SAVE
C ********************************************************************
C *                                                                  *
C * This software is an unpublished work containing confidential and *
C * proprietary information of Birkbeck College. Use, disclosure,    *
C * reproduction and transfer of this work without the express       *
C * written consent of Birkbeck College are prohibited. This notice  *
C * must be attached to all copies or extracts of the software.      *
C *                                                                  *
C * (c) 1995 Oliver Smart & Birkbeck College, All rights reserved    *
C * (c) 1996 Oliver Smart & Birkbeck College, All rights reserved    *
C *                                                                  *
C ********************************************************************
C
C Modification history:
C
C Date  Author          Modification
C 10/95 O.S. Smart      Original version
C 07/96 O.S.S.		trap zero distance between centres
C
C This s/r is for use with a capsule. Centres CCENT1 and CCENT2
C it returns the distance squared from POINT to the closest
C point on the capsule centre line.
C
C To save time CUVEC the unit vector in the direction of the
C channel and CLEN can be supplied - if they are found to be
C zero these will be found in this routine
      DOUBLE PRECISION		CCENT1(3), CCENT2(3), POINT(3),
     &				DIST2, CUVEC(3), CLEN

C internal vbles ***************

C local copies of CUVEC and CLEN
      DOUBLE PRECISION		LCUVEC(3), LCLEN

C the vector from ccent1 to point
      DOUBLE PRECISION		PVEC(3), PLEN

C end of decs ******************

C are CUVEC and CLEN supplied
      IF (CLEN.NE. 0.) THEN
        LCLEN = CLEN
        LCUVEC(1) = CUVEC(1)
        LCUVEC(2) = CUVEC(2)
        LCUVEC(3) = CUVEC(3)
      ELSE
C work these out
        LCUVEC(1) = CCENT2(1) - CCENT1(1)
        LCUVEC(2) = CCENT2(2) - CCENT1(2)
        LCUVEC(3) = CCENT2(3) - CCENT1(3)
        LCLEN = SQRT(LCUVEC(1)**2+LCUVEC(2)**2+LCUVEC(3)**2)
        IF (LCLEN.GT.0.) THEN
          LCUVEC(1) =  LCUVEC(1)/LCLEN
          LCUVEC(2) =  LCUVEC(2)/LCLEN
          LCUVEC(3) =  LCUVEC(3)/LCLEN
        ENDIF
      ENDIF

C find vector from CCENT1 to POINT
      PVEC(1) = POINT(1) - CCENT1(1)
      PVEC(2) = POINT(2) - CCENT1(2)
      PVEC(3) = POINT(3) - CCENT1(3)

C now work out the dot product of this vector to ccent
      CALL dDOT( LCUVEC, PVEC, PLEN)

C If plen is smaller than zero then point is closest to ccent1
C also trap zero dist
      IF (PLEN.LE.0.) THEN
        DIST2 = PVEC(1)**2 + PVEC(2)**2 + PVEC(3)**2
      ELSEIF (PLEN.GE.LCLEN) THEN
C this means that the point is closest to ccent2
        DIST2 =  (POINT(1)-CCENT2(1))**2 +
     &           (POINT(2)-CCENT2(2))**2 +
     &           (POINT(3)-CCENT2(3))**2
      ELSE
C this means point is closest to a position on the centre line
C which is PLEN fro CCENT1.  The distance squared is:
        PVEC(1) = PVEC(1) - PLEN*LCUVEC(1)
        PVEC(2) = PVEC(2) - PLEN*LCUVEC(2)
        PVEC(3) = PVEC(3) - PLEN*LCUVEC(3)
        DIST2 = PVEC(1)**2 + PVEC(2)**2 + PVEC(3)**2
      ENDIF

      RETURN
      END


text:
Jeff MacDonald, a.k.a. Suzie, The Suze, Suzapalooze, etc., is a 4thyear Music major at UC Davis.  He lived in the Lefant House from 20042005 with Users/GregorySmith Greg Smith and Users/BenMcCord Brocktune.  He currently lives somewhere in South Davis.

Jeff played bass in Keep Your Distance and Outgunned.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      DOUBLE PRECISION FUNCTION PHO_BESSI0(X)
C**********************************************************************
C
C      Bessel Function I0
C
C**********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION ax , X , y
      SAVE

      ax = ABS(X)
      IF ( ax.LT.3.75D0 ) THEN
         y = (X/3.75D0)**2
         PHO_BESSI0 = 1.0D0 +
     &                y*(3.5156229D0+y*(3.0899424D0+y*(1.2067492D0+
     &                y*(0.2659732D0+y*(0.360768D-1+y*0.45813D-2)))))
      ELSE
         y = 3.75D0/ax
         PHO_BESSI0 = (EXP(ax)/SQRT(ax))
     &                *(0.39894228D0+y*(0.1328592D-1+y*
     &                (0.225319D-2+y*(-0.157565D-2+
     &                y*(0.916281D-2+y*(-0.2057706D-1+
     &                y*(0.2635537D-1+y*(-0.1647633D-1+y*0.392377D-2))))
     &                ))))
      END IF

      END FUNCTION

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
      SUBROUTINE uadded(BUF,LBUF)
!                                PROCESS UADD
!                                ------- ----
!
!  Document Date: 15 May 1996
!
!  PROCESS UADD adds user given points to every specified seismic trace.
!  The user may specify the start time of the addition.
!
!  A common use of process uadd is to subtract the average trace from all
!  the data traces in GPR (ground penetrating radar).
!
!  Spatial interpolation is not available.
!
!  Each parameter list must be terminated with the word END.  The entire
!  set of UADD parameters must be terminated by the word END.
!
!  THE PARAMETER DICTIONARY
!  --- --------- ----------
!
!  PTS   - A list of points to be added to each seismic trace.  The
!          values must be separated by a blank, tab or new line.
!          The points may be given on as many lines as necessary.
!          When the values are given in an exponential format,
!          positive exponents must have a + rather than a blank.
!          Required. e.g. pts 1 2 3 4 5 6 7 8 9 10 11
!
!  FNO   - The first shot (or rp) to be added.  Shot (rp) numbers
!          must increase monotonically.
!          Preset = 0
!
!  LNO   - The last shot (rp) number to be added.  LNO must be
!          larger than FNO in each list and must increase list to list.
!          Preset = 999999999
!
!  STIME - The start time, in seconds, of the addition.  i.e. the time of
!          the first point given in PTS.
!          Preset = the delay of the trace . e.g. stime 1
!
!  END    - Terminates each parameter list.
!
!
!  Written and copyrighted (c) by:
!  Paul Henkart, Scripps Institution of Oceanography, May 1996
!  ALL RIGHTS RESERVED.
!
!  THE PARAMETER LIST ON DISC IS:
!  WORD 1)  FNO
!       2)  LNO
!       3)  NPTS
!       4)  STIME
!       5)  LPRINT
!       6)  NPTS POINTS    ****  THIS IS A VALIABLE FROM LIST TO LIST!!!
!
!
!  ARGUMENTS:
!  BUF    - A SCRATCH ARRAY AT LEAST 500 32 BIT WORDS LONG.
!  LBUF   - THE SAME ARRAY BUT THE 32 BIT INTEGER EQUIVALENT.  NEEDED
!           BECAUSE PRIME FORTRAN DOESN'T ALLOW EQUIVALENCING OF ARGUMENTS.
!
      PARAMETER (NPARS=6)                                               ! THE NUMBER OF USER PARAMETERS
      PARAMETER (MULTIV=6)                                              ! THE PARAMETER NUMBER OF THE FIRST MULTI-VALUED PARAMETER
      COMMON /UADD/ MUNIT,NLISTS
      DIMENSION BUF(111),LBUF(111)
      CHARACTER*6 NAMES(NPARS)
      CHARACTER*1 TYPES(NPARS)
      DIMENSION LENGTH(NPARS)
      CHARACTER*80 TOKEN
      DIMENSION VALS(NPARS),LVALS(NPARS)
      COMMON /EDITS/ IERROR,IWARN,IRUN,NOW,ICOMPT
      COMMON /sioln1/ cbuf
      COMMON /sioln2/ idum1, idum2, iprint
      CHARACTER*100 cbuf
      INTEGER FNO
!
!
      EQUIVALENCE (FNO,LVALS(1)),
     2            (LNO,LVALS(2)),
     3            (LPRINT,LVALS(3)),
     4            (NPTS,LVALS(4)),
     5            (STIME,VALS(5)),
     6            (PTS,VALS(6))
      DATA NAMES/'FNO   ','LNO   ','LPRINT','NPTS  ','STIME ',
     *           'PTS   '/
      DATA TYPES/'L','L','L','L','F','F'/
      DATA LENGTH/3,3,6,4,5,3/
!****
!****      SET THE PRESETS
!****
      FNO=0
      LNO=999999999
      NPTS=0
      STIME = -1.
      LPRINT=0
      IADDWB=0
      LLNO = -1
      NLISTS=0
      NS=0
      isavep = iprint
!****
!****   GET A PARAMETER FILE
!****
      CALL GETFIL(1,MUNIT,TOKEN,ISTAT)
!****
!****   THE CURRENT COMMAND LINE IN THE SYSTEM BUFFER MAY HAVE THE PARAMETERS.
!****   GET A PARAMETER LIST FROM THE USER.
!****
      NTOKES=1
  100 CONTINUE
      CALL GETOKE(TOKEN,NCHARS)                                         ! GET A TOKEN FROM THE USER PARAMETER LINE
      CALL UPCASE(TOKEN,NCHARS)                                         ! CONVERT THE TOKEN TO UPPERCASE
      IF(NCHARS.GT.0) GO TO 150
      IF(NOW.EQ.1) PRINT 140
  140 FORMAT(' <  ENTER PARAMETERS  >')
      CALL RDLINE                                                       ! GET ANOTHER USER PARAMETER LINE
      NTOKES=0
      GO TO 100
  150 CONTINUE
      NTOKES=NTOKES+1
      DO 190 I=1,NPARS                                                  ! SEE IF IT IS A PARAMETER NAME
      LEN=LENGTH(I)                                                     ! GET THE LEGAL PARAMETER NAME LENGTH
      IPARAM=I                                                          ! SAVE THE INDEX
      IF(TOKEN(1:NCHARS).EQ.NAMES(I)(1:LEN).AND.NCHARS.EQ.LEN) GO TO 200
  190 CONTINUE                                                          ! STILL LOOKING FOR THE NAME
      IF(TOKEN(1:NCHARS).EQ.'END'.AND.NCHARS.EQ.3) GO TO 1000           ! END OF PARAM LIST?
      IF(NS.NE.0) GO TO 230
      PRINT 191, TOKEN(1:NCHARS)
  191 FORMAT(' ***  ERROR  *** UADD DOES NOT HAVE A PARAMETER ',
     *  'NAMED ',A10)
      IERROR=IERROR+1
      GO TO 100
!****
!****    FOUND THE PARAMETER NAME, NOW FIND THE VALUE
!****
  200 NS=0
      iprint = isavep
      NPARAM=IPARAM
  210 CONTINUE                                                          !  NOW FIND THE VALUE
      CALL GETOKE(TOKEN,NCHARS)
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 230                                         ! END OF LINE?
      IF(NOW.EQ.1) PRINT 140                                            ! THIS ALLOWS A PARAMETER TO BE ON A DIFFERENT LINE FROM THE NAME
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 210
  230 CONTINUE
      IF(TYPES(NPARAM).NE.'A') GO TO 240
      IF(NAMES(NPARAM).EQ.'ADDWB'.AND.TOKEN(1:NCHARS).EQ.'YES')
     *    IADDWB=1
      GO TO 100
  240 CONTINUE
      CALL DCODE(TOKEN,NCHARS,AREAL,ISTAT)                              ! TRY AND DECODE IT
      IF(ISTAT.EQ.2) GO TO 420                                          ! =2 MEANS IT IS A NUMERIC
      IERROR=IERROR+1                                                   ! DCODE PRINTED AN ERROR
      GO TO 100
  420 IF(TYPES(NPARAM).EQ.'L') GO TO 500
      IF(NPARAM.LT.MULTIV) GO TO 490                                    !  IS IT A MULTIVALUED PARAMETER
      iprint = 0
      NS=NS+1                                                           !  THE TOKEN WAS A MULTI-VALUED PARAMETER
      npts = ns
      ITEMP=MULTIV
      BUF(NS+ITEMP-1)=AREAL
      GO TO 100
  490 VALS(NPARAM)=AREAL                                                !  FLOATING POINT VALUES
      GO TO 100
  500 CONTINUE                                                          !  32 BIT INTEGER VALUES
      LVALS(NPARAM)=AREAL
      GO TO 100
!****
!****   FINISHED A LIST, NOW DO THE ERROR CHECKING
!****
 1000 CONTINUE
      iprint = isavep
      IF(FNO.GT.LLNO) GO TO 1020
      PRINT 1010
 1010 FORMAT(' ***  ERROR  ***  SHOT AND RP NUMBERS MUST INCREASE.')
      IERROR=IERROR+1
 1020 IF(LNO.GE.FNO) GO TO 1030                                         ! DO THEY INCREASE IN THIS LIST
      PRINT 1010
      IERROR=IERROR+1
 1030 LLNO=LNO
!****
!****      WRITE THE PARAMETER LIST TO DISC
!****
      LBUF(1)=FNO
      LBUF(2)=LNO
      LBUF(3)=NPTS
      BUF(4) = stime
      LBUF(5)=LPRINT
      NWRDS=NPTS+5
      IF(IAND(LPRINT,1).EQ.1)  PRINT *,(LBUF(I),I=1,3),buf(4),
     *   lbuf(5),(BUF(J),J=6,NWRDS)
      CALL WRDISC(MUNIT,BUF,NWRDS)
      NLISTS=NLISTS+1
      LLNO=LNO
      NS=0
      NPTS=0
 2020 CALL GETOKE(TOKEN,NCHARS)                                         ! GET THE NEXT TOKEN
      CALL UPCASE(TOKEN,NCHARS)
      NTOKES=NTOKES+1
      IF(NCHARS.GT.0) GO TO 2030                                        !WAS IT THE END OF A LINE?
      IF(NOW.EQ.1) PRINT 140
      CALL RDLINE                                                       ! GET ANOTHER LINE
      NTOKES=0
      GO TO 2020
 2030 IF(TOKEN(1:NCHARS).NE.'END'.OR.NCHARS.NE.3) GO TO 150
      RETURN                                                            !  FINISHED ALL OF THE PARAMETERS!!!
      END


! Code by: Jostein Brandshoi

! ============================================================================================
! Compile/run on Windows:
! > gfortran -o nonlinear_swe.exe nonlinear_swe.f90
! > nonlinear_swe.exe <output_filename>
!
! Compile/run on Linux:
! > gfortran -o nonlinear_swe nonlinear_swe.f90
! > ./nonlinear_swe <output_filename>
!
! Command line arguments:
!     output_filename : Give a filename (*.dat) to store numerical solutions
!
! Purpose of program:
!   Solve nonlinear rotating shallow water equations
!
! Equations:
!   dU/dt + u * du/dx - fV = -gH_0 * dh/dx
!   dV/dt + u * dv/dx fU = 0
!   dh/dt + h * du/dx + u * dh/dx = 0
!
! Scheme (method of characteristics):
!   u_star = (u_P + 2*c_P + u_Q - 2*c_Q + ((f*dt)/2)*(2*v_R -&
!       f*dt*u_R + v_Q + v_P)) / (2 + ((f*f*dt*dt)/2))
!   v_star = v_R - ((f*dt)/2)*(u_star + u_R)
!   c_star = 0.5*(u_star - u_Q + 2*c_Q - ((f*dt)/2)*(v_star + v_Q))
!
!   u(j, n_next) = (1 - alpha_j) * u_star + alpha_j * 0.0
!   v(j, n_next) = (1 - alpha_j) * v_star + alpha_j * 0.0
!   c(j, n_next) = (1 - alpha_j) * c_star + alpha_j * dsqrt(g * H)
! ============================================================================================

PROGRAM advection_solver
    IMPLICIT NONE
    ! ------------------------------- VARIABLE DECLARATIONS -------------------------------
    double precision, parameter :: g = 10.0              ! Accelereation of gravity [m/s^2]
    double precision, parameter :: f = 1E-4              ! Coriolis parameter [s^(-1)]
    double precision, parameter :: H = 100.0             ! Equillibrium height [m]
    double precision, parameter :: dH = H / 2.0          ! Coefficient in initial condition [m]
    double precision, parameter :: L = 2000E+3           ! Length of solution domain [m]
    double precision, parameter :: kappa = 10.0 / L      ! Coefficient in initial condition [m^(-1)]

    double precision, parameter :: dx = L / 100.0        ! Distance between spatial gridpoints [m]
    double precision, parameter :: dt = 60.0             ! Time between temporal gridpoints [s]

    integer :: j, n                                      ! Counter variables (j: space, n: time)
    integer :: n_temp, n_curr = 1, n_next = 2                          ! Indices for curr and next time
    integer, parameter :: FRS_length = 10                              ! Number of points in FRS-zone
    integer, parameter :: j_max = idnint(L / dx) + 1 + 2 * FRS_length  ! Total number of space points
    integer, parameter :: n_max = 50000                  ! Some arbitrary upper bound for simulation time

    double precision :: C_Pjn, C_Qjn, C_Rjn              ! To hold local "Courant" numbers
    double precision :: u_P, v_P, c_P                    ! To hold variable values at point P
    double precision :: u_Q, v_Q, c_Q                    ! To hold variable values at point Q
    double precision :: u_R, v_R                         ! To hold variable values at point R
    double precision :: alpha_j                          ! "Relaxation function"
    double precision :: u_star, v_star, c_star           ! To hold predictor values

    double precision, dimension(j_max, 2) :: u = 0.0     ! To store u at times n and n + 1 (also sets IC)
    double precision, dimension(j_max, 2) :: v = 0.0     ! To store v at times n and n + 1 (also sets IC)
    double precision, dimension(j_max, 2) :: c = 0.0     ! To store c at times n and n + 1

    character(len = 32) :: output_filename               ! Variable to hold output filename for results
    integer, parameter :: unit_number = 10               ! Unit number for identifying output file
    ! -------------------------------------------------------------------------------------

    ! ---------- GETTING ARGS FROM CMD-LINE AND OPEN .DAT FILE FOR RESULT OUTPUT ----------
    CALL GET_COMMAND_ARGUMENT(1, output_filename)                ! Get cmd-line arg for output file
    open(unit = unit_number, file = output_filename, form = "formatted")  ! Open file that will store results
    ! -------------------------------------------------------------------------------------

    print*, j_max, n_max

    ! ------------------------------ INITIAL CONDITION FOR c ------------------------------
    do j = 1, j_max
        c(j, n_curr) = dsqrt(g * (H + dH*dtanh(kappa*(-FRS_length*dx - L/2 + (j - 1)*dx))))
    end do
    ! -------------------------------------------------------------------------------------

    ! ----------------------------- METHOD OF CHARACTERISTICS -----------------------------
    do n = 0, n_max - 1         ! Loop, in time to compute all three variables u, v, c
        do j = 1, j_max
            C_Pjn = (dt / dx) * dabs(u(j, n_curr) + c(j, n_curr))         ! Coefficient for point P
            C_Qjn = (dt / dx) * dabs((u(j, n_curr) - c(j, n_curr)))   ! Coefficient for point Q
            C_Rjn = (dt / dx) * dabs(u(j, n_curr))                    ! Coefficient for point R

            if (C_Pjn > 1 .or. C_Qjn > 1 .or. C_Rjn > 1) then    ! In case we cannot interpolate closest
                print*, "x_P, x_Q or x_R outside! Exiting..."    ! Give user some info what's happening
                call exit(1)                                     ! Terminate execution of program
            end if

            u_P = (1 - C_Pjn) * u(j, n_curr) + C_Pjn * u(j - 1, n_curr)   ! Value of u at point P
            v_P = (1 - C_Pjn) * v(j, n_curr) + C_Pjn * v(j - 1, n_curr)   ! Value of v at point P
            c_P = (1 - C_Pjn) * c(j, n_curr) + C_Pjn * c(j - 1, n_curr)   ! Value of c at point P

            u_Q = (1 - C_Qjn) * u(j, n_curr) + C_Qjn * u(j + 1, n_curr)   ! Value of u at point Q
            v_Q = (1 - C_Qjn) * v(j, n_curr) + C_Qjn * v(j + 1, n_curr)   ! Value of v at point Q
            c_Q = (1 - C_Qjn) * c(j, n_curr) + C_Qjn * c(j + 1, n_curr)   ! Value of c at point Q

            if (u(j, n_curr) >= 0) then     ! If characteristic has positive slope
                u_R = (1 - C_Rjn) * u(j, n_curr) + C_Rjn * u(j - 1, n_curr)   ! Value of u at point R
                v_R = (1 - C_Rjn) * v(j, n_curr) + C_Rjn * v(j - 1, n_curr)   ! Value of v at point R
            else                            ! If characteristic has negative slope
                u_R = (1 - C_Rjn) * u(j, n_curr) + C_Rjn * u(j + 1, n_curr)   ! Value of u at point R
                v_R = (1 - C_Rjn) * v(j, n_curr) + C_Rjn * v(j + 1, n_curr)   ! Value of v at point R
            end if

            if (j >= 2 .and. j <= j_max - 1) then   ! Compute predictor for all points except outer bounds.
                u_star = (u_P + 2*c_P + u_Q - 2*c_Q + ((f*dt)/2)*(2*v_R -&
                    f*dt*u_R + v_Q + v_P)) / (2 + ((f*f*dt*dt)/2))                  ! Predictor for u
                v_star = v_R - ((f*dt)/2)*(u_star + u_R)                            ! Predictor for u
                c_star = 0.5*(u_star - u_Q + 2*c_Q - ((f*dt)/2)*(v_star + v_Q))     ! Predictor for u
            end if

            if (j < FRS_length + 1) then            ! For left FRS-zone
                alpha_j = 1 - dtanh((dble(j) - 1.0) / 2.0)      ! From 1 at outer bound. towards zero approchnig OB
            else if (j > j_max - FRS_length) then   ! For right FRS-zone
                alpha_j = 1 - dtanh((dble(j_max - j)) / 2.0)    ! From 1 at outer bound. towards zero approchnig OB
            else                                    ! For "real" points
                alpha_j = 0.0   ! Zero for all "real" points including open boundaries
            end if

            u(j, n_next) = (1 - alpha_j) * u_star + alpha_j * 0.0            ! Corrector step for u
            v(j, n_next) = (1 - alpha_j) * v_star + alpha_j * 0.0            ! Corrector step for v

            if (j < FRS_length + 1) then
                c(j, n_next) = (1 - alpha_j) * c_star + alpha_j * dsqrt(g * (H-dH))  ! Corrector c in left FRS-zone
            else if (j > j_max - FRS_length) then
                c(j, n_next) = (1 - alpha_j) * c_star + alpha_j * dsqrt(g * (H+dH))  ! Corrector c in right FRS-zone
            else
                c(j, n_next) = c_star   ! Corrector c in left "real" domain
            end if
        end do

        if (modulo(n, 10) == 0) then    ! Write line-wise to file. Write only every 10'th time level
            write(unit = 10, fmt = "(303f20.14)") u(FRS_length + 1:j_max -&
                FRS_length, n_curr),&
                v(FRS_length + 1:j_max - FRS_length, n_curr),&
                c(FRS_length + 1:j_max - FRS_length, n_curr)
        end if

        n_temp = n_next     ! Set help index to next time step column column
        n_next = n_curr     ! Set next time step index to current ime step column
        n_curr = n_temp     ! Set current time step index to next time step column
    end do
    ! -------------------------------------------------------------------------------------

    close(unit = unit_number)    ! Close file after writing

END PROGRAM advection_solver


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

MODULE files

!-------------------------------------------------------------------------------
! Name:     Files
! Purpose:  Contains FORTRAN units and file names.
! Revised:  10 Sep 2001  Original version.  (T. Otte)
!           03 Oct 2001  Added unit for GRIDDESC file.  (T. Otte)
!           07 Jan 2002  Added explicit file names for input meteorology
!                        file, namelist, grid description, and meteorology
!                        header.  (T. Otte)
!           09 Jun 2003  Removed GRIDBDY2D, GRIDBDY3D, and METBDY2D from
!                        output.  (T. Otte)
!           01 Jul 2004  Restored GRIDBDY2D.  Fixed lengths of output file
!                        names to be strictly CHARACTER*16.  Added flag to
!                        create static output (grid) files.  (T. Otte)
!           29 Nov 2004  Added optional input MM5 "TERRAIN" file to get
!                        fractional land use fields.  (T. Otte)
!           26 May 2005  Removed unused (input) variables GRIDCRO2D, GRIDCRO3D,
!                        LANDCRO2D, METCRO2D, METCRO3D, and METDOT3D.  Changed
!                        names for output files by removing "_G1".  Added new
!                        parameter to define maximum number of input meteorology
!                        files, and increased default maximum to 100.  (T. Otte)
!           09 Apr 2007  Added IMPLICIT NONE.  (T. Otte)
!           02 May 2008  Increased MAX_MM to 300 to accommodate 5-minute output
!                        for one day.  (T. Otte)
!           30 Aug 2011  Changed name of routine from FILE to FILES to avoid
!                        conflict with F90 protected intrinsic names.  Changed
!                        F77 character declarations to F90 standard.  (T. Otte)
!           07 Sep 2011  Updated disclaimer.  (T. Otte)
!           31 May 2012  Increased MAX_MM from 300 to 367 to allow for a full
!                        year to be run (including leap year and previous day
!                        for precipitation calculation).  (T. Otte)
!           10 Feb 2018  Added new optional output files for fractional land
!                        use, layered soil data, and mosaic data.  (T. Spero)
!           18 Dec 2018  Removed support for MM5v3 input files.  Added file
!                        names for netCDF output.  Removed runtime option
!                        to not output time-independent files.  (T. Spero)
!-------------------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER                           :: iutmm
  INTEGER,            PARAMETER     :: iutgd      =  4
  INTEGER,            PARAMETER     :: iutmmi     = 10
  INTEGER,            PARAMETER     :: iutnml     =  8

  INTEGER,            PARAMETER     :: max_mm     = 367

  CHARACTER(LEN=1024)                :: file_gd
  CHARACTER(LEN=1024)                :: file_mm    ( max_mm )
  CHARACTER(LEN=1024), PARAMETER     :: file_nml   = 'namelist.mcip'
  CHARACTER(LEN=1024)                :: file_geo

  CHARACTER(LEN=16),  PARAMETER     :: gridbdy2d  = 'GRID_BDY_2D     '
  CHARACTER(LEN=16),  PARAMETER     :: gridcro2d  = 'GRID_CRO_2D     '
  CHARACTER(LEN=16),  PARAMETER     :: griddot2d  = 'GRID_DOT_2D     '
  CHARACTER(LEN=16),  PARAMETER     :: lufraccro  = 'LUFRAC_CRO      '
  CHARACTER(LEN=16),  PARAMETER     :: metbdy3d   = 'MET_BDY_3D      '
  CHARACTER(LEN=16),  PARAMETER     :: metcro2d   = 'MET_CRO_2D      '
  CHARACTER(LEN=16),  PARAMETER     :: metcro3d   = 'MET_CRO_3D      '
  CHARACTER(LEN=16),  PARAMETER     :: metdot3d   = 'MET_DOT_3D      '
  CHARACTER(LEN=16),  PARAMETER     :: mosaiccro  = 'MOSAIC_CRO      '
  CHARACTER(LEN=16),  PARAMETER     :: soicro     = 'SOI_CRO         '

  CHARACTER(LEN=16),  PARAMETER     :: mcipbdyncf = 'mcip_bdy.nc     '
  CHARACTER(LEN=16),  PARAMETER     :: mcipncf    = 'mcip.nc         '

  INTEGER                           :: cdfid_b
  INTEGER                           :: cdfid_m

END MODULE files

max_stars_repo_name: ClimateIT/CMAQ

max_issues_repo_name: ClimateIT/CMAQ

max_forks_repo_name: ClimateIT/CMAQ


! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! 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 and documentation.

! 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.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE RENORM_ON_MASS ( NVC, EPS1 )

! Renormalizes eigenvectors to unit generalized mass

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ERR, F04, F06
      USE SCONTR, ONLY                :  NDOFL, BLNK_SUB_NAM, WARN_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE PARAMS, ONLY                :  EPSIL, SUPINFO, SUPWARN
      USE SUBR_BEGEND_LEVELS, ONLY    :  RENORM_ON_MASS_BEGEND
      USE CONSTANTS_1, ONLY           :  ZERO, ONE
      USE EIGEN_MATRICES_1 , ONLY     :  GEN_MASS, EIGEN_VEC
      USE MODEL_STUF, ONLY            :  EIG_NORM, MAXMIJ, MIJ_COL, MIJ_ROW
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG

      USE RENORM_ON_MASS_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'RENORM_ON_MASS'

      INTEGER(LONG), INTENT(IN)       :: NVC               ! Number of eigenvectors to be renormalized.
      INTEGER(LONG)                   :: I,J               ! DO loop index
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = RENORM_ON_MASS_BEGEND

      REAL(DOUBLE) , INTENT(IN)       :: EPS1              ! Small number to compare variables against zero
      REAL(DOUBLE)                    :: DEN               ! Normalizing factor in gen mass matrix normalization

      INTRINSIC DSQRT,DABS

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      IF (EIG_NORM /= 'MASS    ') THEN
         WRITE(ERR,1001) EIG_NORM
         IF (SUPINFO == 'N') THEN
            WRITE(F06,1001) EIG_NORM
         ENDIF
      ENDIF

      DO I=1,NVC
         IF (DABS(GEN_MASS(I)) < EPS1) THEN
            WARN_ERR = WARN_ERR + 1
            WRITE(ERR,4301) I, GEN_MASS(I)
            IF (SUPWARN == 'N') THEN
               WRITE(F06,4301) I, GEN_MASS(I)
            ENDIF
            RETURN
         ENDIF
      ENDDO

! Adjust MAXMIJ, the largest off-diag gen mass term. It was originally calculated in subr CALC_GEN_MASS and will change if the
! gen masses have changed as a result of this renormalization

      MAXMIJ = MAXMIJ/(GEN_MASS(MIJ_ROW)*GEN_MASS(MIJ_COL))! NOTE: all gen mass terms checked above for > 0.

! Normalize the eigenvectors so that they produce unit generalized mass and reset the gen masses to unity

      DO J=1,NVC
         DEN = DSQRT(GEN_MASS(J))
         DO I=1,NDOFL
            EIGEN_VEC(I,J) = EIGEN_VEC(I,J)/DEN
         ENDDO
         GEN_MASS(J) = ONE                                 ! Now reset generalized masses to unity
      ENDDO

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
 1001 FORMAT(' *INFORMATION: EIGENVECTORS WILL BE RENORMALIZED BASED ON GEN MASS IN LINK4. THEY WILL BE RENORMALIZED TO ',         &
                             A,' LATER IN LINK5',/)

 4301 FORMAT(' *WARNING    : THE GENERALIZED MASS MATRIX HAS A DIAGONAL TERM THAT IS TOO SMALL TO ALLOW RENORMALIZATION OF THE',   &
                           ' EIGENVECTORS.'                                                                                &
                    ,/,14X,' THE SMALL TERM IS FOR EIGENVECTOR ',I8,' AND ITS VALUE IS ',1ES9.2                                    &
                    ,/,14X,' EIGENVECTORS WILL NOT BE RENORMALIZED TO UNIT MASS IN SUBR RENORM_ON_MASS')

99001 FORMAT(1X,10(1ES13.6))

! **********************************************************************************************************************************

      END SUBROUTINE RENORM_ON_MASS


!
!  Copyright (C) 2013, Northwestern University
!  See COPYRIGHT notice in top-level directory.
!

      module runtime_m
      implicit none

      integer i_time        !time step counter
      integer method        ! 0: blocking APIs, 1: nonblocking
      logical restart       ! new run or restart switch:
                            ! .FALSE. for new run, .TRUE. for restart
      integer i_time_end    !ending time step

      logical io_one_species_at_a_time ! whether read/write one species at a time

      double precision time          !current time
      double precision tstep         !time step (non-dimensional)
      double precision time_save     !time at which to write savefiles (seconds)
      double precision time_save_inc !increment by which to write savefiles (seconds)

      double precision time_ref      !reference time (s)

      character(len=20) run_title    !unique title of run

      end module runtime_m


text:
-makelib ies_lib/xpm -sv \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_cdc/hdl/xpm_cdc.sv" \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib ies_lib/xpm \
  "C:/Xilinx/Vivado/2020.2/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies_lib/microblaze_v11_0_4 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/9285/hdl/microblaze_v11_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_microblaze_0_0/sim/mb_design_microblaze_0_0.vhd" \
-endlib
-makelib ies_lib/lmb_v10_v3_0_11 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/c2ed/hdl/lmb_v10_v3_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_dlmb_v10_0/sim/mb_design_dlmb_v10_0.vhd" \
  "../../../bd/mb_design/ip/mb_design_ilmb_v10_0/sim/mb_design_ilmb_v10_0.vhd" \
-endlib
-makelib ies_lib/lmb_bram_if_cntlr_v4_0_19 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/0b98/hdl/lmb_bram_if_cntlr_v4_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_dlmb_bram_if_cntlr_0/sim/mb_design_dlmb_bram_if_cntlr_0.vhd" \
  "../../../bd/mb_design/ip/mb_design_ilmb_bram_if_cntlr_0/sim/mb_design_ilmb_bram_if_cntlr_0.vhd" \
-endlib
-makelib ies_lib/blk_mem_gen_v8_4_4 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/2985/simulation/blk_mem_gen_v8_4.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_lmb_bram_0/sim/mb_design_lmb_bram_0.v" \
-endlib
-makelib ies_lib/generic_baseblocks_v2_1_0 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/b752/hdl/generic_baseblocks_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_infrastructure_v1_1_0 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/ec67/hdl/axi_infrastructure_v1_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_register_slice_v2_1_22 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/af2c/hdl/axi_register_slice_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_5 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/276e/simulation/fifo_generator_vlog_beh.v" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_5 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/276e/hdl/fifo_generator_v13_2_rfs.vhd" \
-endlib
-makelib ies_lib/fifo_generator_v13_2_5 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/276e/hdl/fifo_generator_v13_2_rfs.v" \
-endlib
-makelib ies_lib/axi_data_fifo_v2_1_21 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/54c0/hdl/axi_data_fifo_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/axi_crossbar_v2_1_23 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/bc0a/hdl/axi_crossbar_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_xbar_0/sim/mb_design_xbar_0.v" \
-endlib
-makelib ies_lib/axi_lite_ipif_v3_0_4 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/66ea/hdl/axi_lite_ipif_v3_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/axi_intc_v4_1_15 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/47b8/hdl/axi_intc_v4_1_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_microblaze_0_axi_intc_0/sim/mb_design_microblaze_0_axi_intc_0.vhd" \
-endlib
-makelib ies_lib/xlconcat_v2_1_4 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/4b67/hdl/xlconcat_v2_1_vl_rfs.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_microblaze_0_xlconcat_0/sim/mb_design_microblaze_0_xlconcat_0.v" \
-endlib
-makelib ies_lib/mdm_v3_2_19 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/8715/hdl/mdm_v3_2_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_mdm_1_0/sim/mb_design_mdm_1_0.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_clk_wiz_1_0/mb_design_clk_wiz_1_0_clk_wiz.v" \
  "../../../bd/mb_design/ip/mb_design_clk_wiz_1_0/mb_design_clk_wiz_1_0.v" \
-endlib
-makelib ies_lib/lib_cdc_v1_0_2 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/ef1e/hdl/lib_cdc_v1_0_rfs.vhd" \
-endlib
-makelib ies_lib/proc_sys_reset_v5_0_13 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/8842/hdl/proc_sys_reset_v5_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_rst_clk_wiz_1_100M_0/sim/mb_design_rst_clk_wiz_1_100M_0.vhd" \
-endlib
-makelib ies_lib/interrupt_control_v3_1_4 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/a040/hdl/interrupt_control_v3_1_vh_rfs.vhd" \
-endlib
-makelib ies_lib/axi_gpio_v2_0_24 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/4318/hdl/axi_gpio_v2_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_axi_gpio_0_0/sim/mb_design_axi_gpio_0_0.vhd" \
  "../../../bd/mb_design/ip/mb_design_axi_gpio_1_0/sim/mb_design_axi_gpio_1_0.vhd" \
  "../../../bd/mb_design/ip/mb_design_axi_gpio_2_0/sim/mb_design_axi_gpio_2_0.vhd" \
-endlib
-makelib ies_lib/lib_pkg_v1_0_2 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/0513/hdl/lib_pkg_v1_0_rfs.vhd" \
-endlib
-makelib ies_lib/lib_srl_fifo_v1_0_2 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/51ce/hdl/lib_srl_fifo_v1_0_rfs.vhd" \
-endlib
-makelib ies_lib/axi_uartlite_v2_0_26 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/5edb/hdl/axi_uartlite_v2_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_axi_uartlite_0_0/sim/mb_design_axi_uartlite_0_0.vhd" \
-endlib
-makelib ies_lib/axi_timer_v2_0_24 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/70d6/hdl/axi_timer_v2_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_axi_timer_0_0/sim/mb_design_axi_timer_0_0.vhd" \
-endlib
-makelib ies_lib/fit_timer_v2_0_10 \
  "../../../../CustomHardwareDisplays.gen/sources_1/bd/mb_design/ipshared/e81b/hdl/fit_timer_v2_0_vh_rfs.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../bd/mb_design/ip/mb_design_fit_timer_0_0/sim/mb_design_fit_timer_0_0.vhd" \
  "../../../bd/mb_design/ipshared/adf4/hdl/CustomDisplays_v1_0_S00_AXI.vhd" \
  "../../../bd/mb_design/ipshared/adf4/src/Nexys4DispDriver.vhd" \
  "../../../bd/mb_design/ipshared/adf4/hdl/CustomDisplays_v1_0.vhd" \
  "../../../bd/mb_design/ip/mb_design_CustomDisplays_0_0/sim/mb_design_CustomDisplays_0_0.vhd" \
  "../../../bd/mb_design/sim/mb_design.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  glbl.v
-endlib

max_stars_repo_name: danielvcorreia/cr

max_issues_repo_name: danielvcorreia/cr

max_forks_repo_name: danielvcorreia/cr


text:
*> \brief \b CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*> \htmlonly
*> Download CLARFGP + dependencies
*>
*> [TGZ]
*>
*> [ZIP]
*>
*> [TXT]
*> \endhtmlonly
*
*  Definition:
*  ===========
*
*       SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
*
*       .. Scalar Arguments ..
*       INTEGER            INCX, N
*       COMPLEX            ALPHA, TAU
*       ..
*       .. Array Arguments ..
*       COMPLEX            X( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> CLARFGP generates a complex elementary reflector H of order n, such
*> that
*>
*>       H**H * ( alpha ) = ( beta ),   H**H * H = I.
*>              (   x   )   (   0  )
*>
*> where alpha and beta are scalars, beta is real and non-negative, and
*> x is an (n-1)-element complex vector.  H is represented in the form
*>
*>       H = I - tau * ( 1 ) * ( 1 v**H ) ,
*>                     ( v )
*>
*> where tau is a complex scalar and v is a complex (n-1)-element
*> vector. Note that H is not hermitian.
*>
*> If the elements of x are all zero and alpha is real, then tau = 0
*> and H is taken to be the unit matrix.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the elementary reflector.
*> \endverbatim
*>
*> \param[in,out] ALPHA
*> \verbatim
*>          ALPHA is COMPLEX
*>          On entry, the value alpha.
*>          On exit, it is overwritten with the value beta.
*> \endverbatim
*>
*> \param[in,out] X
*> \verbatim
*>          X is COMPLEX array, dimension
*>                         (1+(N-2)*abs(INCX))
*>          On entry, the vector x.
*>          On exit, it is overwritten with the vector v.
*> \endverbatim
*>
*> \param[in] INCX
*> \verbatim
*>          INCX is INTEGER
*>          The increment between elements of X. INCX > 0.
*> \endverbatim
*>
*> \param[out] TAU
*> \verbatim
*>          TAU is COMPLEX
*>          The value tau.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date December 2016
*
*> \ingroup complexOTHERauxiliary
*
*  =====================================================================
      SUBROUTINE CLARFGP( N, ALPHA, X, INCX, TAU )
*
*  -- LAPACK auxiliary 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            INCX, N
      COMPLEX            ALPHA, TAU
*     ..
*     .. Array Arguments ..
      COMPLEX            X( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               TWO, ONE, ZERO
      PARAMETER          ( TWO = 2.0E+0, ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            J, KNT
      REAL               ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
      COMPLEX            SAVEALPHA
*     ..
*     .. External Functions ..
      REAL               SCNRM2, SLAMCH, SLAPY3, SLAPY2
      COMPLEX            CLADIV
      EXTERNAL           SCNRM2, SLAMCH, SLAPY3, SLAPY2, CLADIV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, AIMAG, CMPLX, REAL, SIGN
*     ..
*     .. External Subroutines ..
      EXTERNAL           CSCAL, CSSCAL
*     ..
*     .. Executable Statements ..
*
      IF( N.LE.0 ) THEN
         TAU = ZERO
         RETURN
      END IF
*
      XNORM = SCNRM2( N-1, X, INCX )
      ALPHR = REAL( ALPHA )
      ALPHI = AIMAG( ALPHA )
*
      IF( XNORM.EQ.ZERO ) THEN
*
*        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
*
         IF( ALPHI.EQ.ZERO ) THEN
            IF( ALPHR.GE.ZERO ) THEN
*              When TAU.eq.ZERO, the vector is special-cased to be
*              all zeros in the application routines.  We do not need
*              to clear it.
               TAU = ZERO
            ELSE
*              However, the application routines rely on explicit
*              zero checks when TAU.ne.ZERO, and we must clear X.
               TAU = TWO
               DO J = 1, N-1
                  X( 1 + (J-1)*INCX ) = ZERO
               END DO
               ALPHA = -ALPHA
            END IF
         ELSE
*           Only "reflecting" the diagonal entry to be real and non-negative.
            XNORM = SLAPY2( ALPHR, ALPHI )
            TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
            DO J = 1, N-1
               X( 1 + (J-1)*INCX ) = ZERO
            END DO
            ALPHA = XNORM
         END IF
      ELSE
*
*        general case
*
         BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
         SMLNUM = SLAMCH( 'S' ) / SLAMCH( 'E' )
         BIGNUM = ONE / SMLNUM
*
         KNT = 0
         IF( ABS( BETA ).LT.SMLNUM ) THEN
*
*           XNORM, BETA may be inaccurate; scale X and recompute them
*
   10       CONTINUE
            KNT = KNT + 1
            CALL CSSCAL( N-1, BIGNUM, X, INCX )
            BETA = BETA*BIGNUM
            ALPHI = ALPHI*BIGNUM
            ALPHR = ALPHR*BIGNUM
            IF( ABS( BETA ).LT.SMLNUM )
     $         GO TO 10
*
*           New BETA is at most 1, at least SMLNUM
*
            XNORM = SCNRM2( N-1, X, INCX )
            ALPHA = CMPLX( ALPHR, ALPHI )
            BETA = SIGN( SLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
         END IF
         SAVEALPHA = ALPHA
         ALPHA = ALPHA + BETA
         IF( BETA.LT.ZERO ) THEN
            BETA = -BETA
            TAU = -ALPHA / BETA
         ELSE
            ALPHR = ALPHI * (ALPHI/REAL( ALPHA ))
            ALPHR = ALPHR + XNORM * (XNORM/REAL( ALPHA ))
            TAU = CMPLX( ALPHR/BETA, -ALPHI/BETA )
            ALPHA = CMPLX( -ALPHR, ALPHI )
         END IF
         ALPHA = CLADIV( CMPLX( ONE ), ALPHA )
*
         IF ( ABS(TAU).LE.SMLNUM ) THEN
*
*           In the case where the computed TAU ends up being a denormalized number,
*           it loses relative accuracy. This is a BIG problem. Solution: flush TAU
*           to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
*
*           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
*           (Thanks Pat. Thanks MathWorks.)
*
            ALPHR = REAL( SAVEALPHA )
            ALPHI = AIMAG( SAVEALPHA )
            IF( ALPHI.EQ.ZERO ) THEN
               IF( ALPHR.GE.ZERO ) THEN
                  TAU = ZERO
               ELSE
                  TAU = TWO
                  DO J = 1, N-1
                     X( 1 + (J-1)*INCX ) = ZERO
                  END DO
                  BETA = -SAVEALPHA
               END IF
            ELSE
               XNORM = SLAPY2( ALPHR, ALPHI )
               TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
               DO J = 1, N-1
                  X( 1 + (J-1)*INCX ) = ZERO
               END DO
               BETA = XNORM
            END IF
*
         ELSE
*
*           This is the general case.
*
            CALL CSCAL( N-1, ALPHA, X, INCX )
*
         END IF
*
*        If BETA is subnormal, it may lose relative accuracy
*
         DO 20 J = 1, KNT
            BETA = BETA*SMLNUM
 20      CONTINUE
         ALPHA = BETA
      END IF
*
      RETURN
*
*     End of CLARFGP
*
      END

max_stars_repo_name: lingochamp/OpenBLAS

max_issues_repo_name: lingochamp/OpenBLAS

max_forks_repo_name: lingochamp/OpenBLAS


if(allocated(obj%dNdgzi) ) then
    deallocate(obj%dNdgzi)
endif
allocate(obj%dNdgzi(obj%NumOfDim,obj%NumOfNode) )

if(obj%NumOfNode==1) then
    !#########################################################################
    !#######                                                           #######
    !#######                             +     (1)                     #######
    !#######                                                           #######
    !#######                                                           #######
    !#########################################################################
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==1 "
    obj%ierr=1

elseif(obj%NumOfNode==2) then

    if(obj%NumOfDim==1) then

        !#########################################################################
        !#######                                                           #######
        !#######                +>----------------------->+                #######
        !#######            (1)                               (2)          #######
        !#######                                                           #######
        !#########################################################################

        obj%dNdgzi(1,1)= -0.50d0
        obj%dNdgzi(1,2)=  0.50d0

        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction NumOfNode==2 NumOfOrder/=1 "
        obj%ierr=1
    endif
elseif(obj%NumOfNode==3) then
    if(obj%NumOfDim==1) then
        !#########################################################################
        !#######                                                           #######
        !#######                +-----------+-------------+                #######
        !#######            (1)             (2)            (3)             #######
        !#######                                                           #######
        !#########################################################################
        !allocate(obj%dNdgzi(obj%NumOfNode,obj%NumOfDim) )
        !
        !obj%dNdgzi(1,1)=0.50d0*( 1.0d0-gzi(1))
        !obj%dNdgzi(1,2)=0.50d0*(-1.0d0+gzi(1))
        !obj%dNdgzi(1,3)=0.50d0*( 1.0d0-gzi(1))
        !
        !obj%ErrorMsg="Succeed::GetShapeFunction "
        !obj%ierr=0
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    elseif(obj%NumOfDim==2) then
        !#########################################################################
        !#######                                                           #######
        !#######          (1)   +-------------------------+                #######
        !#######                 \                       /  (3)            #######
        !#######                   \                   /                   #######
        !#######                     \               /                     #######
        !#######                       \           /                       #######
        !#######                         \       /                         #######
        !#######                           \   /                           #######
        !#######                       (2)   +                             #######
        !#########################################################################

        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    else

        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif
elseif(obj%NumOfNode==4) then
    if(obj%NumOfDim==1) then
        !#########################################################################
        !#######                                                           #######
        !#######                +-------+---------+-------+                #######
        !#######          (1)          (2)        (3)       (4)            #######
        !#######                                                           #######
        !#########################################################################
        !obj%ErrorMsg="Succeed::GetShapeFunction "
        !obj%ierr=0
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    elseif(obj%NumOfDim==2) then
        !#########################################################################
        !#######                                                           #######
        !#######           (4)  +-------------------------+  (3)           #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######                |                         |                #######
        !#######           (1)  +-------------------------+  (2)           #######
        !#########################################################################

		obj%dNdgzi(1,1) = -1.0d0/4.0d0*(1.0d0-obj%gzi(2))
	    obj%dNdgzi(1,2) =  1.0d0/4.0d0*(1.0d0-obj%gzi(2))
		obj%dNdgzi(1,3) =  1.0d0/4.0d0*(1.0d0+obj%gzi(2))
		obj%dNdgzi(1,4) = -1.0d0/4.0d0*(1.0d0+obj%gzi(2))

		obj%dNdgzi(2,1) = -1.0d0/4.0d0*(1.0d0-obj%gzi(1))
	    obj%dNdgzi(2,2) = -1.0d0/4.0d0*(1.0d0+obj%gzi(1))
		obj%dNdgzi(2,3) =  1.0d0/4.0d0*(1.0d0+obj%gzi(1))
		obj%dNdgzi(2,4) =  1.0d0/4.0d0*(1.0d0-obj%gzi(1))

        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    elseif(obj%NumOfDim==3) then
        !#########################################################################
        !#######                    (4)   +                                #######
        !#######                        /   \                              #######
        !#######                      /       \                            #######
        !#######                    /           \                          #######
        !#######                  /           ___+    (3)                  #######
        !#######                /  ___----'''     \                        #######
        !#######          (1)  +  -                 \                      #######
        !#######                '''-->---____       \                      #######
        !#######                              -->--- +   (2)               #######
        !#########################################################################

        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif

elseif(obj%NumOfNode==5) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==6) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==7) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==8) then
    if(obj%NumOfDim==3) then
        !#########################################################################
        !#######           (8)   +-------------------------+ (7)           #######
        !#######                /|                        /|               #######
        !#######               / |                  (6)  / |               #######
        !#######          (5) +--|----------------------+  |               #######
        !#######              |  |                      |  |               #######
        !#######              |  +----------------------|--+ (3)           #######
        !#######              | / (4)                   | /                #######
        !#######              |/                        |/                 #######
        !#######          (1) +-------------------------+ (2)              #######
        !#########################################################################

        obj%dNdgzi(1,1)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
	    obj%dNdgzi(1,2)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(2))*(1.0d0 - obj%gzi(3))
		obj%dNdgzi(1,3)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
		obj%dNdgzi(1,4)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(2))*(1.0d0 - obj%gzi(3))
        obj%dNdgzi(1,5)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
	    obj%dNdgzi(1,6)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(2))*(1.0d0 + obj%gzi(3))
		obj%dNdgzi(1,7)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))
		obj%dNdgzi(1,8)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(2))*(1.0d0 + obj%gzi(3))

        obj%dNdgzi(2,1)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3))
	    obj%dNdgzi(2,2)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3))
		obj%dNdgzi(2,3)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(3))
		obj%dNdgzi(2,4)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(3))
        obj%dNdgzi(2,5)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3))
	    obj%dNdgzi(2,6)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3))
		obj%dNdgzi(2,7)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(3))
		obj%dNdgzi(2,8)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(3))

        obj%dNdgzi(3,1)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))
	    obj%dNdgzi(3,2)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))
		obj%dNdgzi(3,3)= - 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))
		obj%dNdgzi(3,4)= - 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))
        obj%dNdgzi(3,5)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 - obj%gzi(2))
	    obj%dNdgzi(3,6)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 - obj%gzi(2))
		obj%dNdgzi(3,7)= + 1.0d0/8.0d0*(1.0d0 + obj%gzi(1))*(1.0d0 + obj%gzi(2))
		obj%dNdgzi(3,8)= + 1.0d0/8.0d0*(1.0d0 - obj%gzi(1))*(1.0d0 + obj%gzi(2))

        obj%ErrorMsg="Succeed::GetShapeFunction "
        obj%ierr=0
    else
        print *, "ERROR::GetShapeFunction"
        obj%ErrorMsg="ERROR::GetShapeFunction "
        obj%ierr=1
    endif

elseif(obj%NumOfNode==9) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==10) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==11) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==12) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==13) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==14) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==15) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==16) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==17) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==18) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==19) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==20) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==21) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==22) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==23) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
elseif(obj%NumOfNode==24) then
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
else
    print *, "ERROR::GetShapeFunction"
    obj%ErrorMsg="ERROR::GetShapeFunction "
    obj%ierr=1
endif


text:
      DOUBLE PRECISION FUNCTION PHO_DZEROX(A0,B0,Eps,Maxf,F,Mode)
C**********************************************************************
C
C     Based on
C
C        J.C.P. Bus and T.J. Dekker, Two Efficient Algorithms with
C        Guaranteed Convergence for Finding a Zero of a Function,
C        ACM Trans. Math. Software 1 (1975) 330-345.
C
C        (MODE = 1: Algorithm M;    MODE = 2: Algorithm R)
C
C        CERNLIB C200
C
C***********************************************************************
      IMPLICIT NONE
      DOUBLE PRECISION a , A0 , atl , b , B0 , c , d , Eps , F , fa ,
     &                 fb , fc , fd , fda , fdb , h , HALF , hb , p , q
      DOUBLE PRECISION tol , w , Z1
      INTEGER ie , im1 , im2 , Maxf , mf , Mode
      SAVE

C  input/output channels
      INCLUDE 'inc/poinou'

      CHARACTER NAME*(*)
      PARAMETER (NAME='PHO_DZEROX')
      LOGICAL lmt
      DIMENSION im1(2) , im2(2) , lmt(2)
      EXTERNAL F

      PARAMETER (Z1=1,HALF=Z1/2)

      DATA im1/2 , 3/ , im2/ - 1 , 3/

      IF ( Mode.NE.1 .AND. Mode.NE.2 ) THEN
         c = -2D+10
         IF ( LPRi.GT.4 ) WRITE (LO,99010) NAME , Mode
99010    FORMAT (1X,A,': mode = ',I3,' illegal')
         GOTO 400
      END IF
      fa = F(B0)
      fb = F(A0)
      IF ( fa*fb.GT.0 ) THEN
         c = -3D+10
         IF ( LPRi.GT.4 ) WRITE (LO,99020) NAME
99020    FORMAT (1X,A,': F(A) and F(B) have the same sign')
         GOTO 400
      END IF
      atl = ABS(Eps)
      b = A0
      a = B0
      lmt(2) = .TRUE.
      mf = 2
 100  c = a
      fc = fa
 200  ie = 0
 300  IF ( ABS(fc).LT.ABS(fb) ) THEN
         IF ( c.NE.a ) THEN
            d = a
            fd = fa
         END IF
         a = b
         b = c
         c = a
         fa = fb
         fb = fc
         fc = fa
      END IF
      tol = atl*(1+ABS(c))
      h = HALF*(c+b)
      hb = h - b
      IF ( ABS(hb).GT.tol ) THEN
         IF ( ie.GT.im1(Mode) ) THEN
            w = hb
         ELSE
            tol = tol*SIGN(Z1,hb)
            p = (b-a)*fb
            lmt(1) = ie.LE.1
            IF ( lmt(Mode) ) THEN
               q = fa - fb
               lmt(2) = .FALSE.
            ELSE
               fdb = (fd-fb)/(d-b)
               fda = (fd-fa)/(d-a)
               p = fda*p
               q = fdb*fa - fda*fb
            END IF
            IF ( p.LT.0 ) THEN
               p = -p
               q = -q
            END IF
            IF ( ie.EQ.im2(Mode) ) p = p + p
            IF ( p.EQ.0 .OR. p.LE.q*tol ) THEN
               w = tol
            ELSE IF ( p.LT.hb*q ) THEN
               w = p/q
            ELSE
               w = hb
            END IF
         END IF
         d = a
         a = b
         fd = fa
         fa = fb
         b = b + w
         mf = mf + 1
         IF ( mf.GT.Maxf ) THEN
            IF ( LPRi.GT.4 ) WRITE (LO,99030) NAME
99030       FORMAT (1X,A,': too many function calls')
            GOTO 400
         END IF
         fb = F(b)
         IF ( fb.EQ.0 .OR. SIGN(Z1,fc).EQ.SIGN(Z1,fb) ) GOTO 100
         IF ( w.EQ.hb ) GOTO 200
         ie = ie + 1
         GOTO 300
      END IF
 400  PHO_DZEROX = c

      END FUNCTION


program constructs

    ! Load the OpenMP functions library
    use omp_lib

    ! Set variables
    implicit none
    integer :: num_threads, tnum, i, total

    ! Create parallel block
    !$OMP PARALLEL

    ! Create a section block
    !$OMP SECTIONS PRIVATE(tnum, i)

    ! Ask an available thread to print out the thread number.
    !$OMP SECTION
        tnum = OMP_GET_THREAD_NUM()
        print *, "I am thread number", tnum

    ! Ask another section to add up the thread numbers
    !$OMP SECTION
        num_threads = OMP_GET_NUM_THREADS()
        tnum = OMP_GET_THREAD_NUM()
        total = 0
        DO i = 1, num_threads
            total = total + i
        ENDDO
        print *, "thread number ", tnum, "says total =", total

    ! Close the section block. Normally this sets a barrier that requires all
    ! the threads to have completed processing by this point, but we've
    ! bypassed it with the NOWAIT argument.
    !$OMP END SECTIONS NOWAIT

    ! Print out the fact that the section block has finished. How many threads
    ! are still functional at this point?
    print *, "Finished sections block"

    ! We only want one thread to operate here.
    !$OMP SINGLE
        tnum = OMP_GET_THREAD_NUM()
        print *, "Single thread =", tnum
    !$OMP END SINGLE

    ! End parallel block
    !$OMP END PARALLEL

end program constructs


text:
Trails, Parks, Hiking Areas & Bird Watching

California Hiking and Backpacking Trails and Birds and Bird Watching Bird Watching Spots Organized Geographically by County.

• North Coast & Coast Range Mountain Hiking Trails.
wiki:mendocino:Hiking Del Norte, wiki:mendocino:Hiking Humboldt, wiki:mendocino:Hiking Trinity, wiki:mendocino:Hiking Mendocino, wiki:chico:Hiking_and_Backpacking Tehama, wiki:chico:Hiking_and_Backpacking Glenn, wiki:napa:Hiking Lake, wiki:napa:Hiking Napa, Hiking and Backpacking Colusa, wiki:sonomacounty:Hiking Sonoma, wiki:marin:Hiking Marin.

• Klamath Basin Hiking Trails.
wiki:mendocino:Hiking Trinity, wiki:redding:Hiking Siskiyou, wiki:redding:Hiking Modoc, wiki:redding:Hiking Shasta, wiki:redding:Hiking Lassen.

• Sierra Nevada Mountain Range Hiking Trails.
wiki:chico:Hiking_and_Backpacking Plumas, wiki:nevadacounty:Hiking Sierra, wiki:nevadacounty:Hiking Nevada, wiki:lincoln:Hiking Placer, wiki:ranchocordova:Hiking El Dorado, wiki:ranchocordova:Hiking Amador, wiki:ranchocordova:Hiking Alpine, wiki:modesto:Hiking Calaveras, wiki:modesto:Hiking Tuolumne, wiki:modesto:Hiking Mono, wiki:merced:Hiking Mariposa, wiki:merced:Hiking Madera, wiki:fresno:Hiking Fresno, wiki:visalia:Hiking Tulare, wiki:visalia:Hiking Kern, wiki:PCT:Hiking Pacific Crest Trail.

• Sacramento & San Joaquin Valley Hiking Trails.
wiki:redding:Hiking Shasta, wiki:chico:Hiking_and_Backpacking Tehama, wiki:chico:Hiking_and_Backpacking Glenn, wiki:chico:Hiking_and_Backpacking Butte, Hiking and Backpacking Colusa, Hiking and Backpacking Yolo, wiki:yubasutter:Hiking Sutter, wiki:yubasutter:Hiking Yuba, wiki:lincoln:Hiking Placer, wiki:vacaville:Hiking Solano, wiki:sacramento:Hiking Sacramento, wiki:sacramento:Hiking San Joaquin, wiki:modesto:Hiking Stanislaus, wiki:merced:Hiking Merced, wiki:merced:Hiking Madera, wiki:fresno:Hiking Fresno.

• San Francisco Bay Area Hiking Trails.
wiki:sf:Hiking San Francisco, wiki:marin:Hiking Marin, wiki:sonomacounty:Hiking Sonoma, wiki:napa:Hiking Napa, wiki:vacaville:Hiking Solano, wiki:walnutcreek:Hiking Contra Costa, http://berkeleywiki.org/Hiking Alameda, wiki:sanjose:Hiking Santa Clara, wiki:epa:Hiking San Mateo.

• South Coast & Coast Range Mountain Hiking Trails.
wiki:epa:Hiking San Mateo, wiki:santacruz:Hiking Santa Cruz, wiki:santacruz:Hiking Monterey, wiki:sanjose:Hiking San Benito, wiki:santabarbara:Hiking San Luis Obispo, wiki:santabarbara:Hiking Santa Barbara, wiki:ventura:hiking Ventura, wiki:la:Hiking Los Angeles, wiki:irvine:Hiking Orange, wiki:sandiego:hiking San Diego.

• Inland Empire & Mojave Desert Hiking Trails.
wiki:fresno:Hiking Kings, wiki:visalia:Hiking Tulare, wiki:visalia:Hiking Inyo, wiki:bakersfield:Hiking Kern, wiki:palmsprings:Hiking San Bernardino, wiki:palmsprings:Hiking Riverside, wiki:palmsprings:Hiking Imperial.

 California Hiking Trail Maps
  http://www.mapsofworld.com/usa/countymaps/california/ Maps of Individual Counties in California. http://www.mapsofworld.com/usa/states/california/californiacountymap.html California Map showing county boundaries.
  http://brtinsights.blogspot.com/2008/01/usingtopographicmapsfromacmemapper.html Acme Mapper includes online interactive road,  satellite & topographic maps. It provides geographic coordinates for any point on the map. Enter map coordinates or a landmark name and then jump to that location on the map.

 California Hiking Trail Directory
  http://brtinsights.blogspot.com/2010/02/californiahikingtrailsdirectory.html Hiking Blogs, Personal & Community Websites in California
  http://brtinsights.blogspot.com/2010/03/californiahikingtrailsdirectorypg2.html National, State & Regional Parks in California
  http://brtinsights.blogspot.com/2010/03/californiahikingtrailsdirectorypg3.html County Parks in California
  http://brtinsights.blogspot.com/2010/03/californiahikingtrailsdirectorypg4.html Hiking Trail Lists Organized by Counties in California
  http://brtinsights.blogspot.com/2010/04/californiahikingtrailsdirectorypage.html Hiking Clubs in California
  http://brtinsights.blogspot.com/2010/04/californiahikingtrailsdirectorypg6.html Photos from Hiking in California

 Other California Hiking Resources
  http://wildflowers.russellramblings.com California Wildflower Hikes

20100817 22:00:00 nbsp this page is sweet! Feather falls here I come Users/StevenDaubert


text:
      SUBROUTINE RWQICI
C
C CHANGE RECORD
C READ IN SPATIALLY AND/OR TEMPORALLY VARYING ICS (UNIT INWQICI).
C
      USE GLOBAL
      CHARACTER TITLE(3)*79
      REAL,SAVE,ALLOCATABLE,DIMENSION(:)::XWQV
      IF(.NOT.ALLOCATED(XWQV))THEN
		ALLOCATE(XWQV(NWQVM))
	    XWQV=0.0
	ENDIF
C
      OPEN(1,FILE=ICIFN,STATUS='UNKNOWN')
      OPEN(2,FILE='WQ3D'//ans(partid2)//'.OUT',STATUS='UNKNOWN',POSITION='APPEND')

      WRITE(2,60)'* READING INITIAL CONDITIONS'
      READ(1,50) (TITLE(M),M=1,3)
      WRITE(2,999)
      WRITE(2,50) (TITLE(M),M=1,3)

      READ(1,999)
      READ(1,50) TITLE(1)
      WRITE(2,50) TITLE(1)
      DO M=2,LA
        READ(1,84) I,J,(XWQV(NW),NW=1,NWQV)
        IF(IJCT(I,J).LT.1 .OR. IJCT(I,J).GT.8)THEN
          PRINT*, 'I, J, LINE# = ', I,J,M-1
          STOP 'ERROR!! INVALID (I,J) IN FILE 1'
        ENDIF
        L=LIJ(I,J)
        DO K=1,KC
          DO NW=1,NWQV
            WQV(L,K,NW)=XWQV(NW)
          ENDDO
        ENDDO
        WRITE(2,84) I,J,(WQV(L,1,NW),NW=1,NWQV)
      ENDDO
C
C: WQCHLX=1/WQCHLX
C
      DO L=2,LA
        DO K=1,KC
          WQCHL(L,K) = WQV(L,K,1)*WQCHLC + WQV(L,K,2)*WQCHLD
     &        + WQV(L,K,3)*WQCHLG
          IF(IWQSRP.EQ.1)THEN
            O2WQ_ = MAX(WQV(L,K,19), 0.0)
            WQTAMD = MIN( WQTAMDMX*EXP(-WQKDOTAM*O2WQ_), WQV(L,K,20) )
            WQTAMP(L,K) = WQV(L,K,20) - WQTAMD
            WQPO4D(L,K) = WQV(L,K,10) / (1.0 + WQKPO4P*WQTAMP(L,K))
            WQSAD(L,K)  = WQV(L,K,17) / (1.0 + WQKSAP*WQTAMP(L,K))
          ELSE IF(IWQSRP.EQ.2)THEN
            WQPO4D(L,K) = WQV(L,K,10) / (1.0 + WQKPO4P*SEDT(L,K))
            WQSAD(L,K)  = WQV(L,K,17) / (1.0 + WQKSAP*SEDT(L,K))
          ELSE
            WQPO4D(L,K) = WQV(L,K,10)
            WQSAD(L,K)  = WQV(L,K,17)
          ENDIF
        ENDDO
      ENDDO

      IWQICI = 0
      CLOSE(1)
      CLOSE(2)

  999 FORMAT(1X)
   50 FORMAT(A79)
   52 FORMAT(I7, 1X, A3)
   60 FORMAT(/, A24, I5, A24)
   84 FORMAT(3I5, 21E12.4)
      RETURN
      END

max_stars_repo_name: fearghalodonncha/DeepCurrent

max_issues_repo_name: fearghalodonncha/DeepCurrent

max_forks_repo_name: fearghalodonncha/DeepCurrent


text:
program main
  implicit none
  include 'mpif.h'
  integer i,n,ierr,mpisize,mpirank
  integer,dimension (128) :: iseed
  real(8) diff,norm
  real(8),parameter :: pi=3.14159265358979312d0
  real(8),allocatable,dimension(:) :: x,g,u,ud
  call mpi_init(ierr)
  call mpi_comm_size(mpi_comm_world,mpisize,ierr)
  call mpi_comm_rank(mpi_comm_world,mpirank,ierr)
  n = 10000
  allocate( x(3*n),g(3*n),u(3*n),ud(3*n) )
  do i = 1,128
     iseed(i) = mpirank
  enddo
  call random_seed(put=iseed)
  call random_number(x)
  call random_number(g)
  do i = 1,n
     x(3*i-2) = (x(3*i-2) - 0.5) * pi
     x(3*i-1) = (x(3*i-1) - 0.5) * pi
     x(3*i-0) = (x(3*i-0) - 0.5) * pi
     g(3*i-2) = (g(3*i-2) - 0.5) / n
     g(3*i-1) = (g(3*i-1) - 0.5) / n
     g(3*i-0) = (g(3*i-0) - 0.5) / n
     u(3*i-2) = 0
     u(3*i-1) = 0
     u(3*i-0) = 0
  enddo
  call fmm_init()
  call fmm_partition(n,x,g,u)
  call fmm_biot_savart(n,x,g,u)
  call direct_biot_savart(n,x,g,ud)
  diff = 0
  norm = 0
  do i = 1,n
     diff = diff + (u(3*i-2) - ud(3*i-2)) ** 2
     diff = diff + (u(3*i-1) - ud(3*i-1)) ** 2
     diff = diff + (u(3*i-0) - ud(3*i-0)) ** 2
     norm = norm + ud(3*i-2) ** 2
     norm = norm + ud(3*i-1) ** 2
     norm = norm + ud(3*i-0) ** 2
  enddo
!  call mpi_allreduce(diff,diff,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
!  call mpi_allreduce(norm,norm,1,mpi_double_precision,mpi_sum,0,mpi_comm_world,ierr)
  print '(a,es12.5)',"error         :",sqrt(diff/norm)
  call mpi_finalize(ierr)
  call fmm_finalize()
  deallocate( x,g,u,ud )
end program main

max_stars_repo_name: exafmm/exafmm-alpha

max_issues_repo_name: exafmm/exafmm-alpha

max_forks_repo_name: exafmm/exafmm-alpha


text:
c Routine AS 93 returns frequencies.   The following short routine calculates
c the distribution function from these frequencies (overwriting them).
c The calling arguments are as for AS 93.   The distribution function is
c returned in array A1.   The first element in A1 is F(ASTART).   N.B. ASTART
c is a real variable.
c
      subroutine wprob(test, other, astart, a1, l1, a2, a3, ifault)
      integer test, other, l1, ifault
      real astart, a1(l1), a2(l1), a3(l1)
c
c     Local variables
c
      real zero, sum
      data zero /0.0/
c
      call gscale(test, other, astart, a1, l1, a2, a3, ifault)
      if (ifault .ne. 0) return
c
c     Scale column of F
c
      nrows = 1 + (test * other)/2
      sum = zero
      do 10 i = 1, nrows
      sum = sum + a1(i)
      a1(i) = sum
   10 continue
      do 20 i = 1, nrows
   20 a1(i) = a1(i) / sum
c
      return
      end

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

      SUBROUTINE GSCALE(TEST, OTHER, ASTART, A1, L1, A2, A3, IFAULT)
C
C       ALGORITHM AS 93 APPL. STATIST. (1976) VOL.25, NO.1
C
C       FROM THE SIZES OF TWO SAMPLES THE DISTRIBUTION OF THE
C       ANSARI-BRADLEY TEST FOR SCALE IS GENERATED IN ARRAY A1.
C
      REAL ASTART, A1(L1), A2(L1), A3(L1), AI, ONE, FPOINT
      INTEGER TEST, OTHER
      LOGICAL SYMM
      DATA ONE /1.0/
C
C       TYPE CONVERSION (EFFECT DEPENDS ON TYPE STATEMENT ABOVE).
C
      FPOINT(I) = I
C
C       CHECK PROBLEM SIZE AND DEFINE BASE VALUE OF THE DISTRIBUTION.
C
      M = MIN0(TEST, OTHER)
      IFAULT = 2
      IF (M. LT. 0) RETURN
      ASTART = FPOINT((TEST + 1) / 2) * FPOINT(1 + TEST / 2)
      N = MAX0(TEST, OTHER)
C
C       CHECK SIZE OF RESULT ARRAY.
C
      IFAULT = 1
      LRES = 1 + (M * N) / 2
      IF (L1 .LT. LRES) RETURN
      SYMM = MOD(M + N, 2) .EQ. 0
C
C       TREAT SMALL SAMPLES SEPARATELY.
C
      MM1 = M - 1
      IF (M .GT. 2) GOTO 5
C
C       START-UP PROCEDURES ONLY NEEDED.
C
      IF (MM1) 1, 2, 3
C
C       ONE SAMPLE ONLY.
C
1     A1(1) = ONE
      GOTO 15
C
C       SMALLER SAMPLE SIZE = 1.
C
2      CALL START1(N, A1, L1, LN1)
      GOTO 4
C
C       SMALLER SAMPLE SIZE = 2.
C
3     CALL START2(N, A1, L1, LN1)
C
C     RETURN IF A1 IS NOT IN REVERSE ORDER.
C
4     IF (SYMM .OR. (OTHER .GT. TEST)) GOTO 15
      GOTO 13
C
C       FULL GENERATOR NEEDED
C       SET UP INITIAL CONDITIONS (DEPENDS ON MOD(N, 2)).
C
5     NM1 = N - 1
      NM2 = N - 2
      MNOW = 3
      NC = 3
      IF (MOD(N, 2) .EQ. 1) GOTO 6
C       SET UP FOR EVEN N.
C
      N2B1 = 3
      N2B2 = 2
      CALL START2(N, A1, L1, LN1)
      CALL START2(NM2, A3, L1, LN3)
      CALL START1(NM1, A2, L1, LN2)
      GOTO 8
C
C       SET UP FOR ODD N.
C
6     N2B1 = 2
      N2B2 = 3
      CALL START1(N, A1, L1, LN1)
      CALL START2(NM1, A2, L1, LN2)
C
C       INCREASE ORDER OF DISTRIBUTION IN A1 BY 2
C       (USING A2 AND IMPLYING A3).
C
7     CALL FRQADD(A1, LN1, L1OUT, L1, A2, LN2, N2B1)
      LN1 = LN1 + N
      CALL IMPLY(A1, L1OUT, LN1, A3, LN3, L1, NC)
      NC = NC + 1
      IF (MNOW .EQ. M) GOTO 9
      MNOW = MNOW + 1
C
C       INCREASE ORDER OF DISTRIBUTION IN A2 BY 2 (USING A3).
C
8     CALL FRQADD(A2, LN2, L2OUT, L1, A3, LN3, N2B2)
      LN2 = LN2 + NM1
      CALL IMPLY(A2, L2OUT, LN2, A3, J, L1, NC)
      NC = NC + 1
      IF (MNOW .EQ. M) GOTO 9
      MNOW = MNOW + 1
      GOTO 7
C
C       IF SYMMETRICAL, RESULTS IN A1 ARE COMPLETE.
C
9     IF (SYMM) GOTO 15
C
C       FOR A SKEW RESULT ADD A2 (OFFSET) INTO A1.
C
      KS = (M + 3) / 2
      J = 1
      DO 12 I = KS, LRES
      IF (I .GT. LN1) GOTO 10
      A1(I) = A1(I) + A2(J)
      GOTO 11
10    A1(I) = A2(J)
11    J = J + 1
12    CONTINUE
C
C       DISTRIBUTION IN A1 POSSIBLY IN REVERSE ORDER.
C
      IF (OTHER .LT. TEST) GOTO 15
C
C       REVERSE THE RESULTS IN A1.
C
13    J = LRES
      NDO = LRES / 2
      DO 14 I = 1, NDO
      AI = A1(I)
      A1(I) =A1(J)
      A1(J) = AI
      J = J - 1
14    CONTINUE
C
C       FINAL RESULTS NOW IN A1.
C
15    IFAULT = 0
      RETURN
      END

      SUBROUTINE START1(N, F, L, LOUT)
C
C       ALGORITHM AS 93.1 APPL. STATIST. (1976) VOL.25, NO.1
C
C       GENERATES A 1,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
      REAL F(L), ONE, TWO
      DATA ONE, TWO /1.0, 2.0/
      LOUT = 1 + N / 2
      DO 1 I = 1, LOUT
1     F(I) = TWO
      IF (MOD(N, 2) .EQ. 0) F(LOUT) = ONE
      RETURN
      END
C
      SUBROUTINE START2(N, F, L, LOUT)
C
C       ALGORITHM AS 93.2 APPL. STATIST. (1976) VOL.25, NO.1
C
C       GENERATES A 2,N ANSARI-BRADLEY DISTRIBUTION IN F.
C
      REAL F(L), ONE, TWO, THREE, FOUR
      DATA ONE, TWO, THREE, FOUR /1.0, 2.0, 3.0, 4.0/
C
C       DERIVE F FOR 2, NU, WHERE NU IS HIGHEST EVEN INTEGER
C       LESS THAN OR EQUAL TO N.
C       DEFINE NU AND ARRAY LIMITS.
C
      NU = N - MOD(N, 2)
      J = NU + 1
      LOUT = J
      LT1 = LOUT + 1
      NDO = LT1 / 2
      A = ONE
      B = THREE
C
C       GENERATE THE SYMMETRICAL 2,NU DISTRIBUTION.
C
      DO 1 I = 1, NDO
      F(I) = A
      F(J) = A
      J = J - 1
      A = A + B
      B = FOUR - B
1     CONTINUE
      IF (NU .EQ. N) RETURN
C
C       ADD AN OFFSET 1,N DISTRIBUTION INTO F TO GIVE 2,N RESULT.
C
      NU = NDO + 1
      DO 2 I = NU, LOUT
2     F(I) = F(I) + TWO
      F(LT1) = TWO
      LOUT = LT1
      RETURN
      END
C
      SUBROUTINE FRQADD(F1, L1IN, L1OUT, L1, F2, L2, NSTART)
C
C       ALGORITHM AS 93.3 APPL. STATIST. (1976) VOL.25, NO.1
C
C       ARRAY F1 HAS TWICE THE CONTENTS OF ARRAY F2 ADDED INTO IT
C       STARTING WITH ELEMENTS NSTART AND 1 IN F1 AND F2 RESPECTIVELY.
C
      REAL F1(L1), F2(L2), MUL2
      DATA MUL2 /2.0/
      I2 = 1
      DO 1 I1 = NSTART, L1IN
      F1(I1) = F1(I1) + MUL2 * F2(I2)
      I2 = I2 + 1
1     CONTINUE
      NXT = L1IN + 1
      L1OUT = L2 + NSTART - 1
      DO 2 I1 = NXT, L1OUT
      F1(I1) = MUL2 * F2(I2)
      I2 = I2 + 1
2     CONTINUE
      NSTART = NSTART + 1
      RETURN
      END
C
      SUBROUTINE IMPLY(F1, L1IN, L1OUT, F2, L2, L2MAX, NOFF)
C
C       ALGORITHM AS 93.4 APPL. STATIST. (1976) VOL.25, NO.1
C
C       GIVEN L1IN ELEMENTS OF AN ARRAY F1, A SYMMETRICAL
C       ARRAY F2 IS DERIVED AND ADDED ONTO F1, LEAVING THE
C       FIRST NOFF ELEMENTS OF F1 UNCHANGED AND GIVING A
C       SYMMETRICAL RESULT OF L1OUT ELEMENTS IN F1.
C
      REAL F1(L1OUT), F2(L2MAX), SUM, DIFF
C
C       SET-UP SUBSCRIPTS AND LOOP COUNTER.
C
      I2 = 1 - NOFF
      J1 = L1OUT
      J2 = L1OUT - NOFF
      L2 = J2
      J2MIN = (J2 + 1) / 2
      NDO = (L1OUT + 1) / 2
C
C       DERIVE AND IMPLY NEW VALUES FROM OUTSIDE INWARDS.
C
      DO 6 I1 = 1, NDO
C
C       GET NEW F1 VALUE FROM SUM OF L/H ELEMENTS OF
C       F1 + F2 (IF F2 IS IN RANGE).
C
      IF (I2 .GT. 0) GOTO 1
      SUM = F1(I1)
      GOTO 2
1     SUM = F1(I1) + F2(I2)
C
C       REVISE LEFT ELEMENT OF F1.
C
      F1(I1) = SUM
C
C       IF F2 NOT COMPLETE IMPLY AND ASSIGN F2 VALUES
C       AND REVISE SUBSCRIPTS.
C
2     I2 = I2 + 1
      IF (J2 .LT. J2MIN) GOTO 5
      IF (J1 .LE. L1IN) GOTO 3
      DIFF = SUM
      GOTO 4
3     DIFF = SUM - F1(J1)
4     F2(I1) = DIFF
      F2(J2) = DIFF
      J2 = J2 - 1
C
C       ASSIGN R/H ELEMENT OF F1 AND REVISE SUBSCRIPT.
C
5     F1(J1) = SUM
      J1 = J1 - 1
6     CONTINUE
      RETURN
      END

max_stars_repo_name: FRidh/scipy

max_issues_repo_name: FRidh/scipy

max_forks_repo_name: FRidh/scipy


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

! Modify hotstart.nc based on available data within certain regions

! Inputs:
! hotstart.nc
! elev_ic.gr3
! sal_ic.gr3
! sal_ic_reg.gr3
! tem_ic.gr3
! tem_ic_reg.gr3
! Inputs for Chesapeake Bay obs (come with the script, no need to change)

! Outputs:
! hotstart.nc (overwrite original values with *_ic.gr3 inside *.reg,
! and overwrite original values with Chesapeake Bay obs)

! ifort -O2 -mcmodel=medium -assume byterecl -o modify_hot compute_zcor.f90 modify_hot.f90 -I$NETCDF/include -I$NETCDF_FORTRAN/include -L$NETCDF_FORTRAN/lib -L$NETCDF/lib -lnetcdf -lnetcdff

  program gen_hot
  use netcdf
  use compute_zcor

  integer :: status ! netcdf local status variable
  integer :: sid, ncids(100) ! Netcdf file IDs
  integer :: latvid, lonvid, zmvid, hvid ! positional variables
  integer :: xdid, xvid ! longitude index
  integer :: ydid, yvid ! latitude index
  integer :: ldid, lvid ! vertical level, 1 is top
  integer :: svid, tvid ! salt & temp variable IDs
  integer :: uvid, vvid ! vel variable IDs
  integer :: evid ! SSH variable IDs
  integer, dimension(nf90_max_var_dims) :: dids

  real(4), allocatable :: eout(:),xl(:),yl(:),dp(:), &
    & sigma_lcl(:,:),z(:,:),ztot(:),sigma(:),x0(:),y0(:),val(:,:),z0(:,:),temp1(:), &
    & temp2(:),saltout(:,:),tempout(:,:),tr_nd(:,:,:),tsel(:,:,:)
  integer, allocatable :: i_ic_tem(:),i_ic_sal(:),elnode(:,:),i34(:),kbp2(:),nDepth(:),nn(:), &
    & iTest(:)
  character(Len = 1000),allocatable :: stName(:)

! -----------------Read grids------------------------
  open(17,file='sal_ic_reg.gr3',status='old')
  open(13,file='tem_ic_reg.gr3',status='old')
  open(16,file='hgrid.ll',status='old')
  open(14,file='hgrid.gr3',status='old') !only need depth info and connectivity
  open(19,file='vgrid.in',status='old')
  open(11,file='fort.11',status='replace')
  read(14,*)
  read(14,*)ne,np

  allocate(i_ic_tem(np),i_ic_sal(np),xl(np),yl(np),dp(np),stat=ierr)
  allocate(elnode(4,ne),i34(ne),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (1)'

  read(16,*); read(16,*); read(17,*); read(17,*); read(13,*); read(13,*)
  do i=1,np
    read(14,*)j,xtmp,ytmp,dp(i)
    read(16,*)j,xl(i),yl(i)

    read(17,*)j,xtmp,ytmp,tmp
    i_ic_sal(i)=nint(tmp)
    read(13,*)j,xtmp,ytmp,tmp
    i_ic_tem(i)=nint(tmp)

    if(i_ic_sal(i)<0.and.i_ic_sal(i)>2) then
      write(11,*)'sal ic flag wrong:',i,i_ic_sal(i)
      stop
    endif
    if(i_ic_tem(i)<0.and.i_ic_tem(i)>2) then
      write(11,*)'tem ic flag wrong:',i,i_ic_tem(i)
      stop
    endif
  enddo !i
  do i=1,ne
    read(14,*)j,i34(i),(elnode(l,i),l=1,i34(i))
  enddo !i

  close(14); close(16); close(17); close(13)

  print*, 'xl(1): ',xl(1)
  print*, 'xl(np): ',xl(np)
  print*, 'min/max xl: ',minval(xl),maxval(xl)
  print*, 'dp(1): ',dp(1)
  print*, 'dp(np): ',dp(np)
  print*, 'min/max dp: ',minval(dp),maxval(dp)

!     V-grid
  read(19,*)ivcor
  read(19,*)nvrt
  rewind(19)
  allocate(sigma_lcl(nvrt,np),z(np,nvrt),kbp2(np),stat=ierr)
  allocate(ztot(0:nvrt),sigma(nvrt),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (2)'
  call get_vgrid_single('vgrid.in',np,nvrt,ivcor,kz,h_s,h_c,theta_b,theta_f,ztot,sigma,sigma_lcl,kbp2)

!     Compute z-coord.
  do i=1,np
    if(ivcor==2) then
      call zcor_SZ_single(max(0.11,dp(i)),0.,0.1,h_s,h_c,theta_b,theta_f,kz,nvrt,ztot,sigma,z(i,:),idry,kbp2(i))
    else if(ivcor==1) then
      z(i,kbp2(i):nvrt)=max(0.11,dp(i))*sigma_lcl(kbp2(i):nvrt,i)
    else
      write(11,*)'Unknown ivcor:',ivcor
      stop
    endif

    !Extend below bottom for interpolation later
    z(i,1:kbp2(i)-1)=z(i,kbp2(i))
  enddo !i

  print*, 'z(1,1): ',z(1,1)
  print*, 'z(np,1): ',z(np,1)
  print*, 'min/max z: ',minval(z),maxval(z)

! -----------------Open old hotstart.nc------------------------
  status = nf90_open('hotstart.nc', nf90_nowrite, sid)
  call check(status)

  status = nf90_inq_varid(sid, "eta2", evid); call check(status)
  print*, 'Done reading eout variable ID'
  status = nf90_Inquire_Variable(sid, evid, dimids = dids); call check(status)
  print*, dids(1),dids(2),dids(3),dids(4)

  status = nf90_Inquire_Dimension(sid, dids(1), len = np0); call check(status)
  if (np0.ne.np) then
    print*, 'inconsistent np'
    stop
  endif
  print*, 'np: ',np

  allocate(eout(np),stat=ier)
  if(ierr/=0) stop 'Allocation failed (3)'
  status = nf90_get_var(sid, evid, eout); call check(status)
  print*, 'Done reading eout'
  print*, 'eout(1): ',eout(1)
  print*, 'eout(np): ',eout(np)
  print*, 'min/max eout: ',minval(eout),maxval(eout)

  status = nf90_inq_varid(sid, "tr_nd", svid); call check(status)
  print*, 'Done reading tr_nd variable ID'
  status = nf90_Inquire_Variable(sid, svid, dimids = dids); call check(status)
  print*, dids(1),dids(2),dids(3),dids(4)

  allocate(tr_nd(2,nvrt,np),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (4)'
  status = nf90_get_var(sid, svid, tr_nd); call check(status)
  print*, 'Done reading tr_nd'
  print*, 'tr_nd(1,1,1): ',tr_nd(1,1,1)
  print*, 'tr_nd(2,nvrt,np): ',tr_nd(2,nvrt,np)
  print*, 'min/max tem: ',minval(tr_nd(1,:,:)),maxval(tr_nd(1,:,:))
  print*, 'min/max sal: ',minval(tr_nd(2,:,:)),maxval(tr_nd(2,:,:))

  allocate(tempout(nvrt,np),saltout(nvrt,np),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (5)'
  tempout=tr_nd(1,:,:); saltout=tr_nd(2,:,:)

!--------------diagnostic outputs---------------------
  open(27,file='temp_surf.0.gr3',status='replace')
  open(28,file='salt_bot.0.gr3',status='replace')

  write(27,*)
  write(28,*)
  write(27,*) ne,np
  write(28,*) ne,np
  do i=1,np
    write(27,*)i,xl(i),yl(i),tempout(nvrt,i)
    write(28,*)j,xl(i),yl(i),saltout(1,i)
  enddo !i
  do i=1,ne
    write(27,*)j,i34(i),(elnode(l,i),l=1,i34(i))
    write(28,*)j,i34(i),(elnode(l,i),l=1,i34(i))
  enddo !i

  close(27)
  close(28)

!---------------subroutines--------------------

  iret=nf90_close(sid)

!-------------------modify CB's sal/tem------------------------------
  open(71,file='ChesBay_Data/hot.in',status='old')
  read (71,*) nSt, nSt_J

  allocate(temp1(4),temp2(nvrt),nn(2),iTest(nSt+nSt_J),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (6)'
  allocate(x0(nSt+nSt_J),y0(nSt+nSt_J),nDepth(nSt+nSt_J),val(nSt+nSt_J,200),stat=ierr)
  allocate(z0(nSt+nSt_J,200),stName(nSt+nSt_J),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (7)'

  do k=1,nSt+nSt_J
      read (71,*) x0(k), y0(k), stName(k)
  enddo
  close(71)

  do iter=1,2  !1: salinity; 2: temperature
    iTest=0
    !Read CBP casts
    do k=1,nSt+nSt_J
        if (iter==1) then
            open(72,file='ChesBay_Data/'//trim(adjustl(stName(k)))//'.sal.hot',status='old')
        else
            open(72,file='ChesBay_Data/'//trim(adjustl(stName(k)))//'.tem.hot',status='old')
        endif
        read (72,*) nDepth(k)
        print *, "station: ", k, stName(k)
        do m=1,nDepth(k)
            read(72,*) z0(k,m), val(k,m)
            !write(*,*) z0(k,m), val(k,m)
        enddo
    enddo
    close(72)

    !inside ChesBay, interpolate from CBP casts, FY
    do i=1,np
      !if(dp(i)<=0) then
      !    temp2(:)=0; !stop 'Dry'
      !    cycle
      !endif
      ! find lat between two CB*.* stations
      if(i_ic_sal(i).ne.1 .and. i_ic_tem(i).ne.1) then !outside Ches Bay
        cycle
      else !Chesapeake Bay main stem
          if (yl(i)>=y0(1)) then
              nn(1)=1; nn(2)=1
          elseif (yl(i)<y0(nSt)) then
              nn(1)=nSt; nn(2)=nSt
          else
              do k=1,nSt-1
                  if (yl(i)<y0(k) .and. yl(i)>=y0(k+1)) then
                      nn(1)=k; nn(2)=k+1
                  endif
              enddo
          endif
      endif

      !find z-level
      do k=1,nvrt
          thisS=abs(z(i,k)/z(i,1))
          do j=1,2
              n0=nn(j)
              if (thisS<=z0(n0,1)/z0(n0,nDepth(n0))) then
                  temp1(j)=val(n0,1)
              else
                  do m=1,nDepth(n0)-1
                      S1=z0(n0,m)/z0(n0,nDepth(n0)) !slayer of the lower
                      S2=z0(n0,m+1)/z0(n0,nDepth(n0)) !slayer of the upper
                      if (thisS>=S1 .and. thisS<=S2) then
                          rat1=abs(S2-thisS)/(S2-S1)
                          rat2=abs(S1-thisS)/(S2-S1)
                          temp1(j)=rat1*val(n0,m)+rat2*val(n0,m+1)
                          exit
                      endif
                  enddo
              endif
          enddo !j
          !horizontal
          if (nn(1)==nn(2)) then
              rat1=.5; rat2=.5
          else
              rat1=abs(y0(nn(2))-yl(i))/abs(y0(nn(2))-y0(nn(1)))
              rat2=abs(y0(nn(1))-yl(i))/abs(y0(nn(2))-y0(nn(1)))
          endif
          temp2(k)=rat1*temp1(1)+rat2*temp1(2)
          if (iter==1) then
              saltout(k,i)=rat1*temp1(1)+rat2*temp1(2)
          else
              tempout(k,i)=rat1*temp1(1)+rat2*temp1(2)
          endif
      enddo !k

      if (iTest(nn(1))==0) then
          print *, nn(1), nn(2)
          print *,  temp2(:)
          iTest(nn(1))=1
      endif

    enddo !i=1,np
  enddo ! iter sal tem

  ! other coastal regions
  iTest=0
  open(27,file='elev_ic.gr3',status='old')
  open(28,file='sal_ic.gr3',status='old')
  open(29,file='tem_ic.gr3',status='old')
  read(27,*); read(29,*); read(28,*)
  read(27,*); read(29,*); read(28,*)
  do i=1,np
    read(27,*) j,xtmp,ytmp,tmp1
    read(28,*) j,xtmp,ytmp,tmp2
    read(29,*) j,xtmp,ytmp,tmp3
    eout(i)=tmp1
    if (i_ic_sal(i)==2) then !in other coastal regions
      saltout(1,i)=tmp2
      do k=2,nvrt
        tmp=abs( (z(i,k)-z(i,1))/(z(i,nvrt)-z(i,1))  )
        saltout(k,i)=max(0.0,-2.5*tmp+saltout(1,i))
      enddo
    endif
    if (i_ic_tem(i)==2) then !in other coastal regions
      tempout(:,i)=tmp3
    endif

    if (i_ic_sal(i)==2.or.i_ic_tem(i)==2) then !in other coastal regions
      if (iTest(1)==0.and.dp(i)>15) then
        iTest(1)=1
        print *,'eout in DB at: ',i,xl(i),yl(i),eout(i)
        print *,'saltout in DB: ',i,saltout(:,i)
        print *,'tempout in DB: ',i,tempout(:,i)
        !pause
      endif
    else !outside DB
      if (iTest(2)==0.and.dp(i)>150) then
        iTest(2)=1
        print *,'eout outside DB at: ',i,xl(i),yl(i),eout(i)
        print *,'saltout outside DB: ',i,saltout(:,i)
        print *,'tempout outside DB: ',i,tempout(:,i)
      endif
    endif
  enddo !i
  close(27); close(29); close(28)

  !set salinity above 3m NAVD to be 0
  do i=1,np
    if (dp(i)<-3)then
      saltout(:,i)=0.0
    elseif (dp(i)<0)then
      rat=max(min(1.0, (dp(i)+3.0)/3.0), 0.0)
      saltout(:,i)=saltout(:,i) * rat
    endif
  enddo!i

  !node to element
  allocate(tsel(2,nvrt,ne),stat=ierr)
  if(ierr/=0) stop 'Allocation failed (8)'
  do i=1,ne
    do k=2,nvrt
      tsel(1,k,i)=(sum(tempout(k,elnode(1:i34(i),i)))+sum(tempout(k-1,elnode(1:i34(i),i))))/2/i34(i)
      tsel(2,k,i)=(sum(saltout(k,elnode(1:i34(i),i)))+sum(saltout(k-1,elnode(1:i34(i),i))))/2/i34(i)
    enddo !k
    tsel(1,1,i)=tsel(1,2,i) !mainly for hotstart format
    tsel(2,1,i)=tsel(2,2,i)
  enddo !i

!------------write new values for eta2, tr_nd, tr_nd0, tr_el-----------------
  status = nf90_open('hotstart.nc',nf90_write, sid); call check(status)

  status = nf90_inq_varid(sid, "eta2", evid); call check(status)
  iret=nf90_put_var(sid,evid,dble(eout(1:np)),(/1/),(/np/))

  status = nf90_inq_varid(sid, "tr_nd", svid); call check(status)
  iret=nf90_put_var(sid,svid,dble(tempout(1:nvrt,1:np)),(/1,1,1/),(/1,nvrt,np/))
  iret=nf90_put_var(sid,svid,dble(saltout(1:nvrt,1:np)),(/2,1,1/),(/1,nvrt,np/))
  status = nf90_inq_varid(sid, "tr_nd0", svid); call check(status)
  iret=nf90_put_var(sid,svid,dble(tempout(1:nvrt,1:np)),(/1,1,1/),(/1,nvrt,np/))
  iret=nf90_put_var(sid,svid,dble(saltout(1:nvrt,1:np)),(/2,1,1/),(/1,nvrt,np/))

  status = nf90_inq_varid(sid, "tr_el", svid); call check(status)
  iret=nf90_put_var(sid,svid,dble(tsel(1:2,1:nvrt,1:ne)))

  iret=nf90_close(sid)

!--------------diagnostic outputs---------------------
  open(27,file='temp_surf.gr3',status='replace')
  open(28,file='salt_bot.gr3',status='replace')

  write(27,*)
  write(28,*)
  write(27,*) ne,np
  write(28,*) ne,np
  do i=1,np
    write(27,*)i,xl(i),yl(i),tempout(nvrt,i)
    write(28,*)j,xl(i),yl(i),saltout(1,i)
  enddo !i
  do i=1,ne
    write(27,*)j,i34(i),(elnode(l,i),l=1,i34(i))
    write(28,*)j,i34(i),(elnode(l,i),l=1,i34(i))
  enddo !i

  close(27)
  close(28)

!---------------subroutines--------------------
  contains
!     Internal subroutine - checks error status after each netcdf,
!     prints out text message each time an error code is returned.
  subroutine check(status)
  integer, intent ( in) :: status

  if(status /= nf90_noerr) then
    print *, trim(nf90_strerror(status))
    print *, 'failed to open nc files'
    stop
  end if
  end subroutine check

  end program

max_stars_repo_name: rustychris/schism

max_issues_repo_name: rustychris/schism

max_forks_repo_name: rustychris/schism


text:
Shelter was one of the installations for the Art 154 in the Winter of 2006

20061207 19:28:38 nbsp So cool... it looks like an eagle, but one of the other eyes is all funky. Kind of odd to walk into late at night! Users/ChristyMarsden

20061208 13:08:40 nbsp Beautiful! Gee and Dalelong time no see! Users/EliseKane


text:
! { dg-do compile }
! PR 19929
! Deallocation of pointer components of derived type arrays
program der_ptr_component
    type :: t
        integer, pointer :: p
    end type t
    type(t) :: a(1)

    allocate(a(1)%p)
    deallocate(a(1)%p)

end program der_ptr_component

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab