text:
!----------------------Subroutine--------------------------------------!
!
      SUBROUTINE BANNER
!
!-------------------------Disclaimer-----------------------------------!
!
!     This material was prepared as an account of work sponsored by
!     an agency of the United States Government. Neither the
!     United States Government nor the United States Department of
!     Energy, nor Battelle, nor any of their employees, makes any
!     warranty, express or implied, or assumes any legal liability or
!     responsibility for the accuracy, completeness, or usefulness
!     of any information, apparatus, product, software or process
!     disclosed, or represents that its use would not infringe
!     privately owned rights.
!
!----------------------Acknowledgement---------------------------------!
!
!     This software and its documentation were produced with Government
!     support under Contract Number DE-AC06-76RLO-1830 awarded by the
!     United Department of Energy. The Government retains a paid-up
!     non-exclusive, irrevocable worldwide license to reproduce,
!     prepare derivative works, perform publicly and display publicly
!     by or for the Government, including the right to distribute to
!     other Government contractors.
!
!---------------------Copyright Notices--------------------------------!
!
!            Copyright Battelle Memorial Institute, 1996
!                    All Rights Reserved.
!
!----------------------Description-------------------------------------!
!
!     Display a banner on the screen and output file.
!
!----------------------Authors-----------------------------------------!
!
!     Written by Mark White, December 1992.
!     Last Modified by MD White, Battelle, PNL, December 31, 1992.
!     Last Modified by MD White, PNNL, 30 May 2002.

!     $Id: banner.F 1080 2017-03-14 16:22:02Z d3c002 $
!

!----------------------Fortran 90 Modules------------------------------!
!
      USE GLB_PAR
      USE SOLTN
      USE FILES
!

!----------------------Implicit Double Precision-----------------------!
!
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
!
!----------------------Parameter Statements----------------------------!
!

!
!----------------------Type Declarations-------------------------------!
!
!
!----------------------Common Blocks-----------------------------------!
!

!
!----------------------Executable Lines--------------------------------!
!
      ISUB_LOG = ISUB_LOG+1
      SUB_LOG(ISUB_LOG) = '/BANNER'
      IF( INDEX(SVN_ID(14)(1:1),'$').EQ.0 ) SVN_ID(14) =
     & '$Id: banner.F 1080 2017-03-14 16:22:02Z d3c002 $'
!
!---  Output File  ---
!
      WRITE(IWR,'(A,//)')' Welcome to ...'
      WRITE(IWR,'(A)')   '                           STOMP'
      WRITE(IWR,'(A,//)')'        Subsurface Transport Over Multiple '
     &  // 'Phases'
      WRITE(IWR,'(A)')   ' This file was produced by STOMP, a '
     &  // 'numerical simulator'
      WRITE(IWR,'(A)')   ' developed by the Pacific Northwest '
     &  // 'Laboratory, with'
      WRITE(IWR,'(A)')   ' support from the VOC-Arid Integrated '
     &  // 'Demonstration Project,'
      WRITE(IWR,'(A)')   ' Office of Technology Development, U.S. '
     &  // 'Department of Energy.'
      WRITE(IWR,'(A)')   ' Results from this version of STOMP should '
     &  // 'not be used for'
      WRITE(IWR,'(A,/)') ' license related applications.'
      WRITE(IWR,'(A)')   ' For support:  Tel: 509.372.6070'
      WRITE(IWR,'(A,/)') '               E-mail:  mark.white@pnnl.gov'
      WRITE(IWR,'(A,/)') '                      ---  OUTPUT  ---'

!
!---  Screen Echo ---
!
      WRITE(ISC,'(A,//)')' Welcome to ...'
      WRITE(ISC,'(A)')   '                           STOMP'
      WRITE(ISC,'(A,//)')'        Subsurface Transport Over Multiple '
     &  // 'Phases'
      WRITE(ISC,'(A)')   ' This file was produced by STOMP, a'
     &  // ' numerical simulator'
      WRITE(ISC,'(A)')   ' developed by the Pacific Northwest '
     &  // 'Laboratory, with'
      WRITE(ISC,'(A)')   ' support from the VOC-Arid Integrated '
     &  // 'Demonstration Project,'
      WRITE(ISC,'(A)')   ' Office of Technology Development, U.S. '
     &  // 'Department of Energy.'
      WRITE(ISC,'(A)')   ' Results from this version of STOMP should '
     &  // 'not be used for'
      WRITE(ISC,'(A,/)') ' license related applications.'
      WRITE(ISC,'(A)')   ' For support:  Tel: 509.372.6070'
      WRITE(ISC,'(A,/)') '               E-mail:  mark.white@pnnl.gov'

      ISUB_LOG = ISUB_LOG-1
!
!---  End of BANNER group ---
!
      RETURN
      END

max_stars_repo_name: pnnl/stomp-w

max_issues_repo_name: pnnl/stomp-w

max_forks_repo_name: pnnl/stomp-w


text:
KDVS resident Beat Poetry poet from the Midwest, David is the host of http://kdvs.org/shows/view/show_id/578 Upper Realm Shrieks: Music & Words. He is well loved at the station and on the airwaves for his shows idiosyncratic style, with his theatrical air breaks and the unconventional, eclectic logic of his song progressions. He occasionally publishes poetry in KDViationS.

Hes played music on numerous recordings, including triangles, bell tree, claves, ting, sha, tambourine, industrial kitchen equipment, and toys on Chip Handys P(MM). He often comes to music events around town.

He has many interesting stories to tell from his years spent in Berkeley and Chicago.

He works at the Silo on campus.

David is super cool for sure, trounced some cancer!   You should listen to his show!Users/StevenDaubert

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine fclaw2d_clawpatch5_fort_local_ghost_pack(mx,my,mbc,
     &      meqn, mint,qdata,area,qpack,psize,packmode,ierror)

      implicit none
      integer mx,my,mbc,meqn,psize, mint
      integer pack_area, packmode, ierror
      double precision qdata(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      double precision area(-mbc:mx+mbc+1,-mbc:my+mbc+1)
      double precision qpack(psize)

      integer packq, unpackq, packarea, unpackarea
      parameter(packq = 0, unpackq = 1, packarea = 2,
     &      unpackarea = 3)

      integer i,j,mq,k, ibc,jbc, kfinal
      integer nghost
      logical packdata

      ierror = 0
      if (packmode .ne. packq .and.
     &      packmode .ne. unpackq .and.
     &      packmode .ne. packarea .and.
     &      packmode .ne. unpackarea) then
         ierror = 1
         return
      endif
      packdata = packmode .eq. packq .or. packmode .eq. packarea

      nghost = 2
      k = 1
c     # Face 0
      do mq = 1,meqn
         do j = 1-nghost,my-mint
            do ibc = 1-nghost,mint
               if (packdata) then
                  qpack(k) = qdata(mq,ibc,j)
               else
                  qdata(mq,ibc,j) = qpack(k)
               endif
               k = k + 1
            enddo
         enddo

c        # Face 2
         do jbc = 1-nghost,mint
            do i = mint+1,mx+nghost
               if (packdata) then
                  qpack(k) = qdata(mq,i,jbc)
               else
                  qdata(mq,i,jbc) = qpack(k)
               endif
               k = k + 1
            enddo
         enddo

c        # Face 1
         do j = mint+1,my+nghost
            do ibc = mx-mint+1,mx+nghost
               if (packdata) then
                  qpack(k) = qdata(mq,ibc,j)
               else
                  qdata(mq,ibc,j) = qpack(k)
               endif
               k = k + 1
            enddo
         enddo

c        # Face 3
         do jbc = my-mint+1,my+nghost
            do i = 1-nghost,mx-mint
               if (packdata) then
                  qpack(k) = qdata(mq,i,jbc)
               else
                  qdata(mq,i,jbc) = qpack(k)
               endif
               k = k + 1
            enddo
         enddo
      enddo

      kfinal = k - 1
      if (packmode .ne. packarea .and. packmode .ne. unpackarea) then
         if (kfinal .ne. psize) then
            write(6,*) 'kfinal = ',kfinal
            write(6,*) 'psize = ',psize
            ierror = 2
         endif
         return
      endif

c     # Face 0
      do j = 1-nghost,my-mint
         do ibc = 1-nghost,mint
            if (packdata) then
               qpack(k) = area(ibc,j)
            else
               area(ibc,j) = qpack(k)
            endif
            k = k + 1
         enddo
      enddo

c     # Face 2
      do jbc = 1-nghost,mint
         do i = mint+1,mx+nghost
            if (packdata) then
               qpack(k) = area(i,jbc)
            else
               area(i,jbc) = qpack(k)
            endif
            k = k + 1
         enddo
      enddo

c     # Face 1
      do j = mint+1,my+nghost
         do ibc = mx-mint+1,mx+nghost
            if (packdata) then
               qpack(k) = area(ibc,j)
            else
               area(ibc,j) = qpack(k)
            endif
            k = k + 1
         enddo
      enddo

c     # Face 3
      do jbc = my-mint+1,my+nghost
         do i = 1-nghost,mx-mint
            if (packdata) then
               qpack(k) = area(i,jbc)
            else
               area(i,jbc) = qpack(k)
            endif
            k = k + 1
         enddo
      enddo

      kfinal = k-1
      if (kfinal .ne. psize) then
         ierror = 2
      endif

      end

      subroutine fclaw2d_clawpatch5_set_boundary_to_value(mx,my,mbc,
     &      meqn,q,val)
      implicit none

      integer mx,my,mbc,meqn
      double precision q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc)
      double precision val

      integer i,j,ibc,jbc,mq

c     # set boundary in strips, as with the packing routines

c     # Face 0
      do mq = 1,meqn
         do j = 1-mbc,my
            do ibc = 1,mbc
               q(mq,1-ibc,j) = val
            enddo
         enddo

c        # Face 2
         do jbc = 1,mbc
            do i = 1,mx+mbc
               q(mq,i,1-jbc) = val
            enddo
         enddo

c        # Face 1
         do j = 1,my+mbc
            do ibc = 1,mbc
               q(mq,mx+ibc,j) = val
            enddo
         enddo

c        # Face 3
         do jbc = 1,mbc
            do i = 1-mbc,mx
               q(mq,i,my+jbc) = val
            enddo
         enddo
      enddo

      end

      subroutine fclaw2d_clawpatch5_set_corners_to_value(mx,my,mbc,
     &      meqn, q,value)
      implicit none

      integer mx,my,mbc,meqn
      double precision value
      double precision q(1-mbc:mx+mbc,1-mbc:my+mbc,meqn)

      integer mq,ibc,jbc

      do mq = 1,meqn
         do ibc = mbc,mbc
            do jbc = mbc,mbc
               q(mq,1-ibc,1-jbc) = value
               q(mq,mx+ibc,1-jbc) = value
               q(mq,mx+ibc,my+jbc) = value
               q(mq,1-ibc,my+jbc) = value
            enddo
         enddo
      enddo

      end

max_stars_repo_name: scivision/forestclaw

max_issues_repo_name: scottaiton/forestclaw

max_forks_repo_name: scottaiton/forestclaw


text:
enum Punctuation {
    unknown,

    // Multiple characters punctuation
    // !!!!! Have to be sort by the number of characters
    // !!!!! to avoid bad forward detection (/** issue when last * can override the /* detection during the forward test)

    // TODO the escape of line return can be handle by the parser by checking if there is no more token after the \ on the line
    //    InhibitedLineReturn,    //    \\n or \\r\n          A backslash that preceed a line return (should be completely skipped)
    line_comment,           //    //
    open_block_comment,     //    /*
    close_block_comment,    //    */
    arrow,                  //    ->
    logical_and,            //    &&
    logical_or,             //    ||
    double_colon,           //    ::                 Used for namespaces
    equality_test,          //    ==
    difference_test,        //    !=
    // Not sure that must be detected as one token instead of multiples (especially <<, >>, <<= and >>=) because of templates
    //    LeftShift,              //    <<
    //    RightShift,             //    >>
    //    AdditionAssignment,     //    +=
    //    SubstractionAssignment, //    -=
    //    MultiplicationAssignment, //    *=
    //    DivisionAssignment,     //    /=
    //    DivisionAssignment,     //    %=
    //    DivisionAssignment,     //    |=
    //    DivisionAssignment,     //    &=
    //    DivisionAssignment,     //    ^=
    //    DivisionAssignment,     //    <<=
    //    DivisionAssignment,     //    >>=

    // Mostly in the order on a QWERTY keyboard (symbols making a pair are grouped)
    tilde,                  //    ~                  Should stay the first of single character symbols
    backquote,              //    `
    bang,                   //    !
    at,                     //    @
    hash,                   //    #
    dollar,                 //    $
    percent,                //    %
    caret,                  //    ^
    ampersand,              //    &					 bitwise and
    star,                   //    *
    open_parenthesis,       //    (
    close_parenthesis,      //    )
    underscore,             //    _
    dash,                   //    -
    plus,                   //    +
    equals,                 //    =
    open_brace,             //    {
    close_brace,            //    }
    open_bracket,           //    [
    close_bracket,          //    ]
    colon,                  //    :                  Used in ternaire expression
    semicolon,              //    ;
    single_quote,           //    '
    double_quote,           //    "
    pipe,                   //    |					 bitwise or
    slash,                  //    /
    backslash,              //    '\'
    less,                   //    <
    greater,                //    >
    comma,                  //    ,
    dot,                    //    .
    question_mark,          //    ?

    // White character at end to be able to handle correctly lines that terminate with a separator like semicolon just before a line return
    white_character,
    new_line_character
}

enum Keywords {
}

max_stars_repo_name: Flamaros/f-lang

max_issues_repo_name: Flamaros/f-lang

max_forks_repo_name: Flamaros/f-lang


text: This user only contributes http://daviswiki.org/Users/trma?actionuserinfo negative feedback. I have used the same spa and never had a problem. A woman has to think that only leaving negative feedback is the last resort of the desperate. Is there no business in Davis that deserves your praise? Why even live here.  Users/lrc


text:
module module1
    integer :: n
contains
    recursive subroutine sub1(x)
        integer, intent(inout) :: x
        integer :: y
        y = 0
        if (x < n) then
            x = x + 1
            y = x**2
            print *, 'x = ', x, ', y = ', y
            call sub1(x)
            print *, 'x = ', x, ', y = ', y
        end if
    end
end

program main
    use module1
    integer :: x = 0
    print *, 'Enter number of repeats'
    read (*, *) n
    call sub1(x)
end program main

max_stars_repo_name: leginee/netbeans

max_issues_repo_name: leginee/netbeans

max_forks_repo_name: leginee/netbeans


text:
      BLOCK DATA BLODAT
C     FOR THE EXTRA PARSER-GENERATOR
C     2005-03-29: demingle MEMSETS, PRESEMS, STASTRS
C     GF 19.07.1980: 'SYMPRO' INIT('PROHIB') FOR 'TRACCU', 'TRATTR'
C
      INCLUDE 'PARS.f'
      INCLUDE 'BUCS.f'
      INCLUDE 'CODS.f'
      INCLUDE 'DIRS.f'
      INCLUDE 'ITES.f'
      INCLUDE 'LINS.f'
      INCLUDE 'MEMS.f'
      INCLUDE 'PRES.f'
      INCLUDE 'PRIS.f'
      INCLUDE 'PROS.f'
      INCLUDE 'SEMS.f'
      INCLUDE 'SETS.f'
      INCLUDE 'SLOS.f'
      INTEGER*2 DIR(XXDIRH)
CEQU  EQUIVALENCE (SLOT(1,1),DIR(1))
      INCLUDE 'SPAS.f'
C----------------------------------------------------------------
      INCLUDE 'SPES.f'
      INCLUDE 'STAS.f'
      INCLUDE 'STAT1.f' ! STATICS FOR 'SEMANT'
C--------------------------------------------
      INTEGER*2 TRAFS(15)
      COMMON /STAT2/ TRAFS ! STATICS FOR 'TRAFOR'
      INCLUDE 'STKS.f'
      INCLUDE 'STRS.f'
      INCLUDE 'SYMS.f'
      INCLUDE 'TARS.f'
      INCLUDE 'TRAS.f'
      INCLUDE 'TRES.f'
C
C-------------------------------------------------------------
C     PARS
      DATA UPRI /XXUPRI/
      DATA ULIN /5/ ! /XXULIN/ 2005-04-07
      DATA UTAR /XXUTAR/
      DATA UPAR /XXUPAR/
      DATA UASS /XXUASS/
      DATA PARHIB /XXPARH/
      DATA PARLAS /1/
      DATA PARME(1) /1/
      DATA PARVAL(1) /0/
      DATA PARM(1,1) /'X'/
C----------------------------------------------------------------
C     BUCS
      DATA BUCHIB /XXBUCH/
      DATA BUCKET /XXBUCH * XXSYMH/
C------------------------------------------------------------
C     DIRS
      DATA UDIR/XXUDIR/
      DATA DIRHIB /XXDIRH/
C-----------------------------------------------------------
C     ITES
C     DATA ITEACT(XXITEH) /4/ ! ERROR
      DATA ITESYM(XXITEH) /XXSYMH/ ! FOR MERGING IN 'LACOPY'
      DATA ITEHIB/XXITEH/
      DATA SHIFT,REDUCE,ACCEPT,ERROR
     =     / 1  ,  2   ,  3   ,  4 /
C-------------------------------------------------------------
C     LINS
      DATA LINENO/0/
C-----------------------------------------------------------
C     MEMS
      DATA MEM /XXMEMH * 0/
      DATA FMEM /2/
      DATA MEMHIB/XXMEMH/
      DATA EOP/1/
C-----------------------------------------------------------
C     PRES
      DATA PREHIB/XXPREH/
C-------------------------------------------------------------
C     PRIS
      DATA FPRI /2/ ! POSITION 1 = CARRIAGE CONTROL
      INTEGER*2 XXPRIH2
      PARAMETER (XXPRIH2=XXPRIH*2-1)
      DATA PRIHIB /XXPRIH2/
      DATA PRIBUF /XXPRIH * '  '/
C-----------------------------------------------------------
C     PROS
      DATA PROHIB/XXPROH/
C-----------------------------------------------------------
C     SEMS
      DATA SEMHIB/XXSEMH/
      DATA
     =  ACCO,ACIN,ACMA,ACTA  ! ACCUMULATOR
     = ,ATCO,ATIN,ATMA,SYAT  ! ATTRIBUTE
     = ,SOCO,REIN,UNCH,SOTA  ! SON, MEMBER
     = ,SYCO,SYIN,SYMA,SEPR  ! NEW SYMBOL
     = ,ACCU,ATTR,EOSY,EOS   ! MISCELLANEOUS
     = /   1,   2,   3,   4
     = ,   5,   6,   7,   8
     = ,   9,  10,  11,  12
     = ,  13,  14,  15,  16
     = ,  17,  18,  19,  20
     = /
      INTEGER*2 XXSEMH5
      PARAMETER (XXSEMH5=XXSEMH-5)
C               DESL DESU PNOSEM      FSEM
C                 |    |          |     |
C                 V    V          V     V
      DATA SEM /  0 ,  0 ,  11  , 20 ,  20 , XXSEMH5 *20/
      DATA SEMSYM /XXSEMH * 0/
      DATA PNOSEM /4/ ! -> EOS
      DATA FSEM   /5/ ! -> THEREAFTER
C-----------------------------------------------------------
C     SETS
      DATA SETHIB/XXSETH/
C-----------------------------------------------------------
C     SLOS
      DATA SLOFUL /XXSLOH * 0/
      DATA FSLO /1/
      DATA SLOHIB /XXSLOH/
C----------------------------------------------------------
C     SPAS
      DATA SPAHIB/XXSPAH/
C-----------------------------------------------------------
C     STAS
      DATA STAITE /XXSTAH * 0/
      DATA FNUM /2/ ! = 'FSTA' AT BEGINNING
      DATA STAHIB/XXSTAH/
C----------------------------------------------------------
C     STAT1
      DATA DESL /1/
      DATA BUPROD /0/
      DATA INGRAM /1/ ! START WITH PARSING GRAMMAR
      DATA NEWNUM /7952/
C-----------------------------------------------------------
C     STKS
      DATA STKHIB/XXSTKH/
C-----------------------------------------------------------
C     STRS
      INTEGER*2 XXSTRH2
      PARAMETER (XXSTRH2=2*XXSTRH)
      DATA STRHIB/XXSTRH2/
      DATA FSTR /2/
C-----------------------------------------------------------
C     SYMS
      DATA FSYM /2/
      DATA SYMPOS(1) /1/
CEND  DATA SYMEND(1) /1/
      DATA SYMHIB/XXSYMH/
      DATA SYMPRO /XXSYMH * XXPROH/
      DATA SYMPOS(2) /2/ ! FOR 'HAMAP', 'HAPSE': SYMPOS(FSYM)=FSTR
CEND  DATA SYMEND(2) /2/ ! FOR 'HAMAP', 'HAPSE': SYMPOS(FSYM)=FSTR
C-------------------------------------------------------------
C     TARS
      DATA FTAR /1/ ! POSITION 1 = CARRIAGE CONTROL
      INTEGER*2 XXTARH2
      PARAMETER (XXTARH2=XXTARH*2-1)
      DATA TARHIB /XXTARH2/
      DATA TARBUF /XXTARH * '  '/
C-----------------------------------------------------------
C     TRAS
      DATA TRAHIB/XXTRAH/
      DATA TVOID,TKEYW,TSPEC,TCALL,TCOMT
     =    ,TGOTO,TIDEN,TNUMB,TSTRI
     =    /    1,    2,    3,    4,    5
     =    ,    6,    7,    8,    9
     =    /
      DATA TRASYM/XXTRAH * 0/  ! (NOT PRINTED BY 'TRADUM')
C-------------------------------------------------------------
C     TRES
      DATA UTRE /XXUTRE/
      DATA FTRE /1/
      DATA TREHIB /XXTREH/
      DATA TRECNO /1/
C-------------------------------------------------------------
      END

max_stars_repo_name: gfis/jextra

max_issues_repo_name: gfis/jextra

max_forks_repo_name: gfis/jextra


text:
The mission of the Friends of Adult Day Health Care is to help families to delay the placement of their lovedones in skilled Nursing Homes nursing facilities by providing funds and expertise to Yolo Adult Day Health Center and its diverse programming.

Meetings held 2nd Monday of the month and are held in Davis or Woodland.


text:
    program generate_bathymetry
    implicit none

    integer, parameter :: im = 500
    integer, parameter :: jm = 90
    integer :: i,j
    real(8) :: dx,dy,x,y,h0,h02,tht,cs,sn,pi,  &
               ytmp,xp,yp,ctmp1,ctmp2,depth(im,jm)
    character(len=80) :: test_case

    pi = 2.*acos(0.)

    test_case = 'landslide'

    if(trim(test_case)=='submerged_bar') then
      dx = 0.02
      do j = 1,jm
      do i = 1,im
        x = (i-1)*dx
        if(x<=6.0) then
          depth(i,j) = 0.4
        elseif(x>6.0.and.x<=12.0) then
          depth(i,j) = 0.4-(x-6.0)/20.0
        elseif(x>12.0.and.x<=14.0) then
          depth(i,j) = 0.1
        elseif(x>14.0.and.x<17.0) then
          depth(i,j) = 0.4-(17.0-x)/10.0
        else
          depth(i,j) = 0.4
        endif
      enddo
      enddo
    elseif(trim(test_case)=='landslide') then
      dx = 0.02
      do j = 1,jm
      do i = 1,im
        x = (i-1)*dx+0.5*dx
        if(x<=1.5/tan(15.*pi/180.)+0.1) then
          depth(i,j) = (x-0.1)*tan(15.*pi/180.)
        else
          depth(i,j) = 1.5
        endif
      enddo
      enddo
    elseif(trim(test_case)=='breaking_solitary_wave') then
      dx = 0.02
      do j = 1,jm
      do i = 1,im
        x = (i-1)*dx
        if(x<=4.0) then
          depth(i,j) = 0.21
        else
          depth(i,j) = 0.21-(x-4.0)/20.
        endif
      enddo
      enddo
    elseif(trim(test_case)=='dam_break') then
      dx = 0.05
      do j = 1,jm
      do i = 1,im
        x = (i-1)*dx
        if(x<=25.5) then
          depth(i,j) = 0.01
        elseif(x>25.5.and.x<=28.5) then
          depth(i,j) = -(x-25.5)*0.4/3.0
        elseif(x>28.5.and.x<=31.5) then
          depth(i,j) = -(31.5-x)*0.4/3.0
        elseif(x>31.5) then
          depth(i,j) = 0.01
        endif
      enddo
      enddo
    elseif(trim(test_case)=='tidal_wave') then
      dx = 7.5
      do j = 1,jm
      do i = 1,im
        x = (i-1)*dx
        if(abs(x-750)<=187.5) then
          depth(i,j) = 8.0
        else
          depth(i,j) = 16.0
        endif
      enddo
      enddo
    elseif(trim(test_case)=='Berkhoff_Shoal') then
      dx = 0.025
      dy = 0.05
      h0 = .45
      h02 = 0.07
      tht = 70.*3.14159/180
      cs    = cos(tht)
      sn    = sin(tht)
      ytmp  = 5.84-(h0-h02)/.02
      do i = 1, im
      do j = 1, jm
        x  = float(i-1)*dx- 12.0
        y  = float(j-1)*dy- 10.0
        xp =  x*cs+y*sn
        yp = -x*sn+y*cs
        ctmp1 = (xp/4.)**2+(yp/3.)**2
        ctmp2 = (xp/5.)**2+(yp/3.75)**2
        if(yp.ge.5.84) then
          depth(i,j) = h0
        elseif (yp.le.ytmp) then
          depth(i,j) = h02
        elseif (ctmp1.gt.1.) then
          depth(i,j) = h0-0.02*(5.84-yp)
        else
          depth(i,j) = .75-.02*(5.84-yp)-0.5*sqrt(1.-ctmp2)
        endif
      enddo
      enddo
    endif

    open(11,file='depth.txt')
    do j = 1,jm
      write(11,'(5000f8.4)') (depth(i,j),i=1,im)
    enddo
    close(11)

    end

max_stars_repo_name: NHWAVE/NHWAVE

max_issues_repo_name: NHWAVE/NHWAVE

max_forks_repo_name: NHWAVE/NHWAVE


text:
#include "scalar.fpp"

!=====================================================================!
! Parent class for solving n-th order differential equations. Specific
! integrators extend this class.
!
! Author: Komahan Boopathy (komahan@gatech.edu)
!=====================================================================!

module integrator_interface

  use dynamic_physics_interface, only : dynamics
  !use dynamic_analysis_interface, only : dynamic_analysis

  implicit none

  private
  public ::  integrator

  !-------------------------------------------------------------------!
  ! Define the type
  !-------------------------------------------------------------------!

  type, abstract :: integrator !, extends(dynamic_analysis)

     class(dynamics), allocatable :: system
     type(scalar)   , allocatable :: time(:)  ! time values (steps)
     type(scalar)   , allocatable :: U(:,:,:) ! state varibles (steps, deriv_ord, nvars)
     type(scalar)   , allocatable :: X(:,:)   ! xpoints
     type(scalar)                 :: tinit
     type(scalar)                 :: tfinal
     type(scalar)                 :: h
     type(logical)                :: implicit
     type(integer)                :: num_stages
     type(integer)                :: num_time_steps
     type(integer)                :: total_num_steps

   contains

     procedure                                    :: construct, destruct
     procedure                                    :: solve
     procedure(step_interface), deferred          :: step
     procedure(get_bandwidth_interface), deferred :: get_bandwidth
     procedure                                    :: write_solution
     procedure                                    :: to_string

  end type integrator

  !-------------------------------------------------------------------!
  ! Define interfaces to deferred procedures
  !-------------------------------------------------------------------!

  interface

     impure subroutine step_interface(this, t, u, h, p, ierr)

       import integrator

       class(integrator) , intent(inout) :: this
       type(scalar)      , intent(inout) :: t(:)
       type(scalar)      , intent(inout) :: u(:,:,:)
       type(integer)     , intent(in)    :: p
       type(scalar)      , intent(in)    :: h
       type(integer)     , intent(out)   :: ierr

     end subroutine step_interface

     pure type(integer) function get_bandwidth_interface(this, time_index) result(width)

       import integrator

       class(integrator), intent(in) :: this
       type(integer)    , intent(in) :: time_index

     end function get_bandwidth_interface

  end interface

contains

  !===================================================================!
  ! Base class constructor logic
  !===================================================================!

  subroutine construct(this, system, tinit, tfinal, h, implicit, num_stages)

    class(integrator) , intent(inout) :: this
    class(dynamics)   , intent(in)    :: system
    type(scalar)      , intent(in)    :: tinit, tfinal
    type(scalar)      , intent(in)    :: h
    type(logical)     , intent(in)    :: implicit
    type(integer)     , intent(in)    :: num_stages

    ! Set parameters
    allocate(this % system, source = system)
    this % tinit           = tinit
    this % tfinal          = tfinal
    this % h               = h
    this % num_time_steps  = floor((this % tfinal - this % tinit)/this % h) + 1
    this % num_stages      = num_stages
    this % total_num_steps = this % num_time_steps*(this % num_stages+1) - this % num_stages ! the initial step does not have stages
    this % implicit        = implicit

  end subroutine construct

  !===================================================================!
  ! Base class destructor
  !===================================================================!

  pure subroutine destruct(this)

    class(integrator), intent(inout) :: this

    ! Clear global states and time
    if(allocated(this % U)) deallocate(this % U)
    if(allocated(this % time)) deallocate(this % time)
    if(allocated(this % system)) deallocate(this % system)
    if(allocated(this % X)) deallocate(this % x)

  end subroutine destruct

  !===================================================================!
  ! Write solution to file
  !===================================================================!

  subroutine write_solution(this, filename)

    class(integrator)             :: this
    character(len=*), intent(in)  :: filename
    character(len=7), parameter   :: directory = "output/"
    character(len=:), allocatable :: path
    character(len=:), allocatable :: new_name
    integer                       :: k, j, i, ierr

    ! Open resource
    path = trim(filename)

    open(unit=90, file=trim(path), iostat= ierr)
    if (ierr .ne. 0) then
       write(*,'("  >> Opening file ", 39A, " failed")') path
       return
    end if

    write(90, *) "time ", "u ", "udot"
    ! Write data
    loop_time: do k = 1, this % total_num_steps
       if (mod(k-1,this % num_stages+1) .eq. 0) then
          ! Generalize this based on differential order
          write(90, *) this % time(k), this % U (k,1,:), this % U (k,2,:)!, this % U (k,3,:)
       end if
    end do loop_time

!!$    ! Write data
!!$    write(90, *) "time ", "x ", "y ", "z ",  "state "
!!$    loop_vars : do j = 1, this % system % get_num_state_vars()
!!$       loop_time: do k = 1, this % num_time_steps !total_num_steps !
!!$          write(90, *) this % time((this % num_stages+1)*k - this % num_stages), &
!!$               & this % system % x(1,j), &
!!$               & this % system % x(2,j), &
!!$               & this % system % x(3,j), &
!!$               & this % U ((this % num_stages+1)*k - this % num_stages, 1, j)
!!$       end do loop_time
!!$    end do loop_vars

    ! Close resource
    close(90)

  end subroutine write_solution

  !===================================================================!
  ! Time integration logic
  !===================================================================!

  impure subroutine solve( this )

    class(integrator), intent(inout) :: this

    integer :: k, p
    integer :: ierr

    ! State and time history
    if (allocated(this % time)) deallocate(this%time)
    if (allocated(this % U)) deallocate(this%U)

    allocate(this % time(this % total_num_steps))
    this % time = 0.0d0
    this % time(1) = this % tinit

    allocate( this % U( &
         & this % total_num_steps, &
         & this % system % get_differential_order() + 1, &
         & this % system % get_num_state_vars() &
         & ))
    this % U = 0.0d0

    ! Get the initial condition
    call this % system % get_initial_condition(this % U(1,:,:), this % X)

    ! March in time
    time: do k = 2, this % total_num_steps

       p = this % get_bandwidth(k)

       call this % step(this % time(k-p:k) , &
            & this % U(k-p:k,:,:), &
            & this % h, &
            & p, &
            & ierr)

       !! if (k.eq.50) stop

    end do time

  end subroutine solve

  !===================================================================!
  ! Prints important fields of the class
  !===================================================================!

  subroutine to_string(this)

    class(integrator), intent(in) :: this

    print '("  >> Physical System      : " ,A10)' , this % system % get_description()
    print '("  >> Num state variables  : " ,i8)'  , this % system % get_num_state_vars()
    print '("  >> Equation order       : " ,i4)'  , this % system % get_differential_order()
    print '("  >> Start time           : " ,F8.3)', this % tinit
    print '("  >> End time             : " ,F8.3)', this % tfinal
    print '("  >> Step size            : " ,E9.3)', this % h
    print '("  >> Number of time steps : " ,i10)' , this % num_time_steps
    print '("  >> Number of stages     : " ,i10)' , this % num_stages
    print '("  >> Tot. Number of steps : " ,i10)' , this % total_num_steps

  end subroutine to_string

end module integrator_interface

max_stars_repo_name: komahanb/time-marching-adjoint

max_issues_repo_name: komahanb/diffusion

max_forks_repo_name: komahanb/diffusion


text:
!**********************************************************************************************************************************
! LICENSING
! Copyright (C) 2021  National Renewable Energy Laboratory
!
!    This file is part of ROSCO.
!
! 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.
!**********************************************************************************************************************************
MODULE SysSubs

    USE ROSCO_Types
    USE Constants

    IMPLICIT NONE

! This module contains routines with system-specific logic and references, including all references to the console unit, CU.
! It also contains standard (but not system-specific) routines it uses.
! SysGnuLinux.f90 is specifically for the GNU Fortran (gfortran) compiler on Linux and macOS.

    INTEGER, PARAMETER            :: ConRecL     = 120                               ! The record length for console output.
    INTEGER, PARAMETER            :: CU          = 6                                 ! The I/O unit for the console.  Unit 6 causes ADAMS to crash.
    INTEGER, PARAMETER            :: MaxWrScrLen = 98                                ! The maximum number of characters allowed to be written to a line in WrScr
    LOGICAL, PARAMETER            :: KBInputOK   = .TRUE.                            ! A flag to tell the program that keyboard input is allowed in the environment.
    CHARACTER(*),  PARAMETER      :: NewLine     = ACHAR(10)                         ! The delimiter for New Lines [ Windows is CHAR(13)//CHAR(10); MAC is CHAR(13); Unix is CHAR(10) {CHAR(13)=\r is a line feed, CHAR(10)=\n is a new line}]
    CHARACTER(*),  PARAMETER      :: OS_Desc     = 'GNU Fortran for Windows'         ! Description of the language/OS
    CHARACTER( 1), PARAMETER      :: PathSep     = '\'                               ! The path separator.
    CHARACTER( 1), PARAMETER      :: SwChar      = '/'                               ! The switch character for command-line options.
    CHARACTER(11), PARAMETER      :: UnfForm     = 'UNFORMATTED'                     ! The string to specify unformatted I/O files.

    CONTAINS
!=======================================================================
    SUBROUTINE LoadDynamicLib ( DLL, ErrStat, ErrMsg )

        ! This SUBROUTINE is used to dynamically load a DLL.

        TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be loaded.
        INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
        CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

        INTERFACE  ! Definitions of Windows API routines

           !...........................
           !bjj: I have been unable to find a solution that works with both IVF and gfortran...
           !bjj: note that "Intel Fortran does not support use of STDCALL with BIND(C) at this time"
           !     See this link: http://software.intel.com/en-us/articles/replacing-intel-fortran-attributes-with-c-interoperability-features
           !bjj: Until this is fixed, Intel uses kernel32.f90 definitions instead of the interface below:
           !...........................

           FUNCTION LoadLibrary(lpFileName) BIND(C,NAME='LoadLibraryA')
              USE, INTRINSIC :: ISO_C_BINDING
              IMPLICIT NONE
              !GCC$ ATTRIBUTES STDCALL :: LoadLibrary
              INTEGER(C_INTPTR_T)        :: LoadLibrary
              CHARACTER(KIND=C_CHAR)     :: lpFileName(*)
           END FUNCTION LoadLibrary

        END INTERFACE

        ErrStat = ErrID_None
        ErrMsg = ''

        ! Load the DLL and get the file address:
        DLL%FileAddr = LoadLibrary( TRIM(DLL%FileName)//C_NULL_CHAR )  !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end)
        IF ( DLL%FileAddr == INT(0,C_INTPTR_T) ) THEN
           ErrStat = ErrID_Fatal
           WRITE(ErrMsg,'(I2)') BITS_IN_ADDR
           ErrMsg  = 'The dynamic library '//TRIM(DLL%FileName)//' could not be loaded. Check that the file '// &
                    'exists in the specified location and that it is compiled for '//TRIM(ErrMsg)//'-bit applications.'
           RETURN
        END IF

        ! Get the procedure address:
        CALL LoadDynamicLibProc ( DLL, ErrStat, ErrMsg )

        RETURN
     END SUBROUTINE LoadDynamicLib
     !=======================================================================
     SUBROUTINE LoadDynamicLibProc ( DLL, ErrStat, ErrMsg )

        ! This SUBROUTINE is used to dynamically load a procedure in a DLL.

        TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be loaded.
        INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
        CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
        INTEGER(IntKi)                            :: i

        INTERFACE  ! Definitions of Windows API routines

           !...........................
           !bjj: I have been unable to find a solution that works with both IVF and gfortran...
           !bjj: note that "Intel Fortran does not support use of STDCALL with BIND(C) at this time"
           !     See this link: http://software.intel.com/en-us/articles/replacing-intel-fortran-attributes-with-c-interoperability-features
           !bjj: Until this is fixed, Intel uses kernel32.f90 definitions instead of the interface below:
           !...........................

           FUNCTION GetProcAddress(hModule, lpProcName) BIND(C, NAME='GetProcAddress')
              USE, INTRINSIC :: ISO_C_BINDING
              IMPLICIT NONE
              !GCC$ ATTRIBUTES STDCALL :: GetProcAddress
              TYPE(C_FUNPTR)             :: GetProcAddress
              INTEGER(C_INTPTR_T),VALUE  :: hModule
              CHARACTER(KIND=C_CHAR)     :: lpProcName(*)
           END FUNCTION GetProcAddress

        END INTERFACE

        ErrStat = ErrID_None
        ErrMsg = ''

        ! Get the procedure addresses:
        do i=1,NWTC_MAX_DLL_PROC
           if ( len_trim( DLL%ProcName(i) ) > 0 ) then
              DLL%ProcAddr(i) = GetProcAddress( DLL%FileAddr, TRIM(DLL%ProcName(i))//C_NULL_CHAR )  !the "C_NULL_CHAR" converts the Fortran string to a C-type string (i.e., adds //CHAR(0) to the end)
              IF(.NOT. C_ASSOCIATED(DLL%ProcAddr(i))) THEN
                 ErrStat = ErrID_Fatal + i - 1
                 ErrMsg  = 'The procedure '//TRIM(DLL%ProcName(i))//' in file '//TRIM(DLL%FileName)//' could not be loaded.'
                 RETURN
              END IF
           end if
        end do

        RETURN
     END SUBROUTINE LoadDynamicLibProc
     !=======================================================================
     SUBROUTINE FreeDynamicLib ( DLL, ErrStat, ErrMsg )

        ! This SUBROUTINE is used to free a dynamically loaded DLL (loaded in LoadDynamicLib).

        TYPE (ExtDLL_Type),           INTENT(INOUT)  :: DLL         ! The DLL to be freed.
        INTEGER(IntKi),            INTENT(  OUT)  :: ErrStat     ! Error status of the operation
        CHARACTER(*),              INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None

        ! Local variable:
        INTEGER(C_INT)                            :: Success     ! Whether or not the call to FreeLibrary was successful
        INTEGER(C_INT), PARAMETER                 :: FALSE  = 0

        INTERFACE  ! Definitions of Windows API routines

           FUNCTION FreeLibrary(hLibModule) BIND(C, NAME='FreeLibrary')
              USE, INTRINSIC :: ISO_C_BINDING
              IMPLICIT NONE
              !GCC$ ATTRIBUTES STDCALL :: FreeLibrary
              INTEGER(C_INT)             :: FreeLibrary ! BOOL
              INTEGER(C_INTPTR_T),VALUE  :: hLibModule ! HMODULE hLibModule
           END FUNCTION

        END INTERFACE

        ! Free the DLL:
        IF ( DLL%FileAddr == INT(0,C_INTPTR_T) ) RETURN

        Success = FreeLibrary( DLL%FileAddr ) !If the function succeeds, the return value is nonzero. If the function fails, the return value is zero.

        IF ( Success == FALSE ) THEN !BJJ: note that this isn't the same as the Fortran LOGICAL .FALSE.
           ErrStat = ErrID_Fatal
           ErrMsg  = 'The dynamic library could not be freed.'
           RETURN
        ELSE
           ErrStat = ErrID_None
           ErrMsg = ''
           DLL%FileAddr = INT(0,C_INTPTR_T)
        END IF

        RETURN
     END SUBROUTINE FreeDynamicLib
     !=======================================================================

END MODULE SysSubs

max_stars_repo_name: ptrbortolotti/ROSCO

max_issues_repo_name: ptrbortolotti/ROSCO

max_forks_repo_name: mvanv/ROSCO


text:
(quantity)^2^ is an entertaining comic drawn by Katherine Harnack Katherine Harnack, and published in the California Aggie during the 200506 and 200607 school years.

About:
The name comes from the term used in mathematics to describe when everything inside a set of parentheses is squared.  It is also a play on words, where the characters and jokes in the comic are square. Miss Harnack is a selfproclaimed nerd and member of Davis Anime Club. Like any other comic artist for the Aggie, she does appreciate any email from those who enjoy her comic. Her comics, as well as other art, can be found on her deviantART page http://kathyvicki.deviantart.com/ here.  A compilation of her comics is available for sale at the UC Davis Bookstore.

Characters:
  Cowlick: The typical UC Davis student, apparently his real name. He is currently in possession of a VCard.
  1Up: A typical video game nerd whose real name is John Hochburg. He too is currently in possession of a VCard.
  Demi: The serious academic archetype character. VCard status is unknown.
  Stoic: A play on words, this character is sheerly unshakable. VCard status is unknown, and given her unreadable nature, may never be known.
  Seamus: An Irish exchange student. To help indicate the accent, Seamus speech bubbles contain a different font, and the spelling follows European English rules, i.e: honour.
  Blake the Freshman: A lost and patheticlooking freshman. Since UC Davis accepted far more students than it has beds for in 2006, Blake is living in a very crowded dorm.
  Lissa: 1Ups older sister. She is a foxy version of him, but with more social skills.
  Stan: Lissas fiancé and a lvl 12 Ranger.
  The Nameless Professor: A professor that the gang seems to have regardless of what the course is.
  Nudist Uncle Roy: 1Ups uncle who is also a nudist.

Art:
Everything but the speech bubbles are sketched, inked, and grayscaled on paper with Micron pens, Sharpies, and Prismacolor markers.  After being scanned into a computer, the speech bubbles are added with an inexpensive program called Comic Life Deluxe.

200708 Rejection:
On October 12th, 2007, Harnack received an email from Eddie Lee Eddie Lee, the editor of the California Aggie, telling her that (quantity)^2^ would not be published in the newspaper for its third year (Harnacks senior year at UC Davis).  Harnack emailed him back asking for a reason why the comic was rejected; as of yet no response has been given.  A copy of her application portfolio can be found http://kathyvicki.deviantart.com/art/The10ThatFailed67343729 here.

20070116 13:52:46 nbsp Something funky with lists. Users/DavidPoole

20070121 12:50:04 nbsp I worked on the list, and it seems okay now. Thanks for the adding the picture! Users/KatherineHarnack


text:
module ppum
use types, only: dp
use constants, only: pi
use bsplines, only: bspline, bspline_der, bspline_der2
use mesh, only: meshexp
use quadrature, only: gauss_pts, gauss_wts
use utils, only: stop_error, assert
use linalg, only: eigh
implicit none
private
public do_ppum_basis

contains

    elemental function h(r, rc)
    ! C^3 cutoff function h(r)
    real(dp), intent(in) :: r  ! Point at which to evaluate, must be r >= 0
    real(dp), intent(in) :: rc ! Cutoff radius; h(r) = 0 for r >= rc
    real(dp) :: h
    if (r < rc) then
        h = 1 + 20*(r/rc)**7-70*(r/rc)**6+84*(r/rc)**5-35*(r/rc)**4
    else
        h = 0
    end if
    end function

    elemental function hp(r, rc)
    ! C^3 cutoff function h'(r)
    real(dp), intent(in) :: r  ! Point at which to evaluate, must be r >= 0
    real(dp), intent(in) :: rc ! Cutoff radius; h(r) = 0 for r >= rc
    real(dp) :: hp
    if (r < rc) then
        hp = 140*r**6/rc**7 - 420*r**5/rc**6 + 420*r**4/rc**5 - 140*r**3/rc**4
    else
        hp = 0
    end if
    end function

    subroutine do_ppum_basis(p, xmin, xmax, Ne, penr, npenr, alpha, ortho, Nq, &
            eps, ft_only, xq, wq, B, Bp)
    integer, intent(in) :: p, Ne
    integer, intent(in) :: penr ! polynomial order of the poly enrichment, e.g.
        ! penr = 3 means cubic polynomials (i.e. 4 basis functions spanned by:
        ! 1, x, x^2, x^3)
    integer, intent(in) :: npenr ! number of non-polynomial enrichments, e.g.
        ! npenr = 2 means two non-polynomial enrichments will be added to the
        ! polynomial basis/enrichments of the order penr.
    integer, intent(in) :: Nq, ortho
    real(dp), intent(in) :: xmin, xmax, alpha, eps
    logical, intent(in) :: ft_only
    real(dp), allocatable, intent(out) :: xq(:), wq(:), B(:,:), Bp(:,:)
    logical, allocatable :: Bactive(:,:)
    integer ::  N_intervals, Nq_total
    real(dp) :: rmin, rmax, dx
    real(dp) :: xa, xb, jac, x0
    real(dp), allocatable :: t(:), xiq(:), wtq(:), x(:), hq(:)
    real(dp), allocatable :: W(:,:), Wp(:,:), Wpp(:,:), S(:), Sp(:), Spp(:)
    real(dp), allocatable :: wi(:,:), wip(:,:), wipp(:,:)
    real(dp), allocatable :: Sm(:,:), eigs(:), vecs(:,:)
    real(dp), allocatable :: enr(:,:,:), enrp(:,:,:), enrn(:,:,:), enrnp(:,:,:)
    real(dp), allocatable :: mesh(:)
    real(dp) :: rc
    integer :: i, j, u, bindex, Nb
    integer :: n, k, m, Nenr
    Nenr = (penr+1) + npenr
    k = p+1
    if (mod(p, 2) == 0) then
        n = k+2
    else
        if (p == 1) then
            n = k+1
        else
            n = k+3
        end if
    end if
    bindex = n/2+1
    print *, "Evaluating basis functions"
    print *, "Constructing B-spline n =", n, ", k =", k
    N_intervals = n-k+1

    allocate(mesh(2*Ne))
    dx = (xmax-xmin)/Ne*(alpha-1)/2
    mesh(1) = xmin
    do i = 1, Ne-1
        mesh(2*i)   = (xmax-xmin)/Ne * i + xmin - dx
        mesh(2*i+1) = (xmax-xmin)/Ne * i + xmin + dx
    end do
    mesh(2*Ne) = xmax

    Nq_total = Nq*(2*Ne-1)

    allocate(xq(Nq_total), wq(Nq_total))
    allocate(t(n+k), xiq(Nq), wtq(Nq), x(Nq))
    allocate(W(Nq_total,Ne), Wp(Nq_total,Ne), Wpp(Nq_total,Ne))
    allocate(S(Nq_total), wi(Nq_total,Ne))
    allocate(Sp(Nq_total), wip(Nq_total,Ne))
    allocate(Spp(Nq_total), wipp(Nq_total,Ne))
    allocate(enr(Nq_total,Nenr,Ne))
    allocate(enrp(Nq_total,Nenr,Ne))
    allocate(enrn(Nq_total,Nenr,Ne))
    allocate(enrnp(Nq_total,Nenr,Ne))
    allocate(Sm(Nenr,Nenr))
    allocate(eigs(Nenr))
    allocate(vecs(Nenr,Nenr))
    allocate(Bactive(Nenr,Ne))
    allocate(hq(Nq_total))

    ! Loop over the mesh, and constract a global quadrature rule. Integrals of a
    ! function hq evaluated at the points xq are calculated using: sum(wq*hq)
    xiq = gauss_pts(Nq)
    wtq = gauss_wts(Nq)
    do i = 1, 2*Ne-1
        xa = mesh(i)
        xb = mesh(i+1)
        jac = (xb-xa)/2
        x = (xiq(:)+1) * jac + xa
        xq((i-1)*Nq+1:i*Nq) = x
        wq((i-1)*Nq+1:i*Nq) = wtq*jac
    end do

    ! Construct weight functions Wi(x)
    do i = 1, Ne
        dx = (xmax-xmin)/Ne*(alpha-1)/2
        rmin = (xmax-xmin)/Ne * (i-1) + xmin - dx
        rmax = (xmax-xmin)/Ne * i + xmin + dx
        t(:k-1) = rmin
        t(k:n+1) = meshexp(rmin, rmax, 1._dp, N_intervals)
        t(n+2:) = rmax
        W(:,i)   = bspline     (t, bindex, k, xq)
        Wp(:,i)  = bspline_der (t, bindex, k, xq)
        Wpp(:,i) = bspline_der2(t, bindex, k, xq)
    end do

    ! Construct PU functions using Shepard constructions
    S = 0
    Sp = 0
    Spp = 0
    do i = 1, Ne
        S = S + W(:,i)
        Sp = Sp + Wp(:,i)
        Spp = Spp + Wpp(:,i)
    end do
    do i = 1, Ne
        wi(:,i) = W(:,i)/S
        wip(:,i) = (Wp(:,i)*S - W(:,i)*Sp)/S**2
        wipp(:,i) = (Wpp(:,i)*S**2 - 2*Wp(:,i)*S*Sp - W(:,i)*S*Spp &
            + 2*W(:,i)*Sp**2)/S**3
    end do

    ! Construct enrichment functions enr(x)
    enr = 0
    enrp = 0
    do i = 1, Ne
        dx = (xmax-xmin)/Ne*(alpha-1)/2
        rmin = (xmax-xmin)/Ne * (i-1) + xmin - dx
        rmax = (xmax-xmin)/Ne * i + xmin + dx
        x0 = (rmin+rmax)/2
        x0 = min(abs(rmin),abs(rmax),abs(x0))
        jac = (rmax-rmin)/2
        do m = 0, penr ! m is the polynomial order, e.g. m=2 is quadratic
            j = m+1 ! j is the basis function number, starts from 1
            enr(:,j,i) = legendre_p((xq-rmin)/jac-1, m) * &
                sqrt(m + 1._dp/2)/sqrt(jac)
            enrp(:,j,i) = legendre_p_der((xq-rmin)/jac-1, m)/jac * &
                sqrt(m + 1._dp/2)/sqrt(jac)
            where (xq < rmin .or. xq > rmax)
                enr(:,j,i) = 0
                enrp(:,j,i) = 0
            end where
        end do
        do m = 1, npenr ! loop over non-poly enrichments
            j = m + penr+1 ! total enrichment basis function index
            rc = 300._dp
            select case(m)
                case (1)
                    if (abs(x0) <= 5) then
                        enr(:,j,i) = exp(-xq**2/2) / pi**(1._dp/4) &
                            * h(abs(xq), rc)
                        enrp(:,j,i) = (hp(abs(xq), rc)*sign(1._dp, xq) &
                            - xq*h(abs(xq), rc)) * exp(-xq**2/2) / pi**(1._dp/4)
                    else
                        enr(:,j,i) = 0
                        enrp(:,j,i) = 0
                    end if
                case (2)
                    enr(:,j,i) = xq*exp(-xq**2/2) / pi**(1._dp/4)*sqrt(2._dp)
                    enrp(:,j,i) = (1-xq**2)*exp(-xq**2/2) &
                        / pi**(1._dp/4)*sqrt(2._dp)
                case default
                    call stop_error("non-poly enrichment index not implemented")
            end select
            where (xq < rmin .or. xq > rmax)
                enr(:,j,i) = 0
                enrp(:,j,i) = 0
            end where
        end do
    end do

    ! Orthogonalize the enrichment enr(x) -> enrn(x)
    print *, "Orthogonalization (eigs per element):"
    Bactive = .false.
    do i = 1, Ne
        dx = (xmax-xmin)/Ne*(alpha-1)/2
        rmin = (xmax-xmin)/Ne * (i-1) + xmin
        rmax = (xmax-xmin)/Ne * i + xmin
        if (.not. ft_only) then
            rmin = rmin - dx
            rmax = rmax + dx
        end if
        do j = 1, Nenr
        do k = 1, Nenr
            hq = enr(:,j,i)*enr(:,k,i)*wq
            where (xq < rmin) hq = 0
            where (xq > rmax) hq = 0
            Sm(j,k) = sum(hq)
        end do
        end do
        call eigh(Sm, eigs, vecs)
        print "(200es8.1)", eigs
        do j = 1, Nenr
            enrn(:,j,i) = 0
            enrnp(:,j,i) = 0
            do k = 1, Nenr
                enrn (:,j,i) = enrn (:,j,i) + vecs(k,j)*enr (:,k,i)
                enrnp(:,j,i) = enrnp(:,j,i) + vecs(k,j)*enrp(:,k,i)
            end do
            if (eigs(j) > eps) then
                Bactive(j,i) = .true.
                enrn (:,j,i) = enrn (:,j,i) / sqrt(eigs(j))
                enrnp(:,j,i) = enrnp(:,j,i) / sqrt(eigs(j))
            end if
        end do
    end do

    !do i = 1, Ne
    !    do j = 1, Nenr
    !    do k = 1, Nenr
    !        Sm(j,k) = sum(enrn(:,j,i)*enrn(:,k,i)*wq)
    !    end do
    !    print *, j, Sm(j,:)
    !    end do
    !    call eigh(Sm, eigs, vecs)
    !    print *, eigs
    !    print *
    !end do
    !stop "OK"

    ! Use the orthogonalized functions
    select case (ortho)
        case (0)
            ! Do not use orthogonalized functions
            Bactive = .true.
        case (1)
            ! use orthogonalized functions
            enr = enrn
            enrp = enrnp
        case default
            call stop_error("Invalid 'ortho'.")
    end select

    ! Construct basis functions B = wi*enr
    allocate(B(Nq_total,count(Bactive)))
    allocate(Bp(Nq_total,count(Bactive)))
    Nb = 0
    do i = 1, Ne
        do j = 1, Nenr
            if (Bactive(j,i)) then
                Nb = Nb+1
                B(:,Nb)  = wi(:,i)*enr(:,j,i)
                Bp(:,Nb) = wip(:,i)*enr(:,j,i)+wi(:,i)*enrp(:,j,i)
            end if
        end do
    end do
    call assert(Nb == count(Bactive))

    !open(newunit=u, file="pu.txt", status="replace")
    !write(u,*) xq
    !do i = 1, Ne
    !    write(u,*) W(:,i)
    !end do
    !do i = 1, Ne
    !    write(u,*) wi(:,i)
    !end do
    !do i = 1, Ne
    !    write(u,*) wip(:,i)
    !end do
    !do i = 1, Ne
    !    write(u,*) wipp(:,i)
    !end do
    !close(u)

    !open(newunit=u, file="enr.txt", status="replace")
    !do i = 1, Ne
    !    do j = 1, Nenr
    !        write(u,*) enr(:,j,i)
    !        write(u,*) enrp(:,j,i)
    !    end do
    !end do
    !close(u)

    !open(newunit=u, file="ppum.txt", status="replace")
    !do i = 1, Ne
    !    do j = 1, Nenr
    !        write(u,*) B(:,j,i)
    !        write(u,*) Bp(:,j,i)
    !    end do
    !end do
    !close(u)
    end subroutine

    function legendre_p(x, n) result(r)
    real(dp), intent(in) :: x(:)
    integer, intent(in) :: n
    real(dp) :: r(size(x))
    select case (n)
        case (0)
            r = 1
        case (1)
            r = x
        case (2)
            r = (3*x**2-1)/2
        case (3)
            r = (5*x**3-3*x)/2
        case (4)
            r = (35*x**4-30*x**2+3)/8
        case (5)
            r = (63*x**5-70*x**3+15*x)/8
        case (6)
            r = (231*x**6-315*x**4+105*x**2-5)/16
        case (7)
            r = (429*x**7-693*x**5+315*x**3-35*x)/16
        case (8)
            r = (6435*x**8-12012*x**6+6930*x**4-1260*x**2+35)/128
        case (9)
            r = (12155*x**9-25740*x**7+18018*x**5-4620*x**3+315*x)/128
        case (10)
            r = (46189*x**10-109395*x**8+90090*x**6-30030*x**4+3465*x**2-63)/256
        case default
            call stop_error("n is too high")
    end select
    end function

    function legendre_p_der(x, n) result(r)
    real(dp), intent(in) :: x(:)
    integer, intent(in) :: n
    real(dp) :: r(size(x))
    select case (n)
        case (0)
            r = 0
        case (1)
            r = 1
        case (2)
            r = (6*x)/2
        case (3)
            r = (15*x**2-3)/2
        case (4)
            r = (140*x**3-60*x)/8
        case (5)
            r = (315*x**4-210*x**2+15)/8
        case (6)
            r = (1386*x**5-1260*x**3+210*x)/16
        case (7)
            r = (3003*x**6-3465*x**4+945*x**2-35)/16
        case (8)
            r = (51480*x**7-72072*x**5+27720*x**3-2520*x)/128
        case (9)
            r = (109395*x**8-180180*x**6+90090*x**4-13860*x**2+315)/128
        case (10)
            r = (461890*x**9-875160*x**7+540540*x**5-120120*x**3+6930*x)/256
        case default
            call stop_error("n is too high")
    end select
    end function

end module

program test_ppum
use types, only: dp
use ppum, only: do_ppum_basis
use linalg, only: eigh
use utils, only: stop_error, assert
use schroed_util, only: lho
implicit none
integer :: ppu, Ne, penr, npenr, Nq, Nq_total, i, j, u, ortho, Nb
real(dp) :: alpha, xmin, xmax, eps, condA, condB
real(dp), allocatable :: xq(:), wq(:)
real(dp), allocatable :: B(:,:), Bp(:,:), eigs(:)
logical :: ft_only

ppu = 3
penr = 3
npenr = 1
Ne = 1
alpha = 1.5_dp
xmin = -10
xmax = 10
ortho = 1
eps = 1e-15_dp
Nq = 64
ft_only = .false.

open(newunit=u, file="ppum_conv.txt", status="replace")
do i = 1, 10
    call do_ppum_basis(ppu, xmin, xmax, Ne, penr, npenr, alpha, ortho, Nq, &
        eps, ft_only, xq, wq, B, Bp)
    Nb = size(B,2)
    print *, "Nb =", Nb
    call lho(xq, wq, B, Bp, eigs, condA, condB)
    call assert(size(eigs) >= 4)
    write(u,*) Nb, condA, condB, eigs(:4)

    ! Adjust this to ~250 for a longer convergence study
    if (Nb > 70) exit
    Ne = Ne*2
end do
close(u)

end program

max_stars_repo_name: certik/hfsolver

max_issues_repo_name: certik/hfsolver

max_forks_repo_name: certik/hfsolver


text:
! ==========================================================================
!
! Basic VMC/DMC code for educational purposes.
! This file: global constants and utility subroutines
!
! Copyright (C) 2012 Jindrich Kolorenc
!
! The software is released under MIT/X11 license.
!
! ==========================================================================

module types_const
  implicit none

  integer, parameter :: sp=kind(1.0)
  integer, parameter :: dp=kind(1.0d0)

  integer, parameter :: i2b=selected_int_kind(4)  ! 2 byte int < 2**15
  integer, parameter :: i4b=selected_int_kind(9)  ! 4 byte int < 2**31
  integer, parameter :: i8b=selected_int_kind(10) ! 8 byte int < 2**63

  real(dp), parameter :: pi_d=3.141592653589793238462643383279502884197_dp
  real(dp), parameter :: twopi_d=6.283185307179586476925286766559005768394_dp
  real(dp), parameter :: sqrt2_d=1.414213562373095048801688724209698078570_dp
  real(dp), parameter :: sqrt3_d=1.732050807568877293527446341505872366943_dp
  real(dp), parameter :: sqrt5_d=2.236067977499789696409173668731276235441_dp
  real(sp), parameter :: pi=3.141592653589793238462643383279502884197_sp
  real(sp), parameter :: twopi=6.283185307179586476925286766559005768394_sp
  real(sp), parameter :: sqrt2=1.414213562373095048801688724209698078570_sp
  real(dp), parameter :: sqrt3=1.732050807568877293527446341505872366943_sp
  real(dp), parameter :: sqrt5=2.236067977499789696409173668731276235441_sp

contains

  subroutine missing(var)
    ! {{{ error message for reading input files.
    character(len=*) var
    print *, "Missing variable '", trim(adjustl(var)), "' in the input file."
    stop
    ! }}}
  end subroutine missing

  subroutine finish_line(test,text)
    ! {{{ finish the line after a non-advancing write
    logical, intent(in) :: test
    character(len=*), intent(in) :: text
    if ( test ) then
       print *
    else
       print *, " ", text
    end if
    ! }}}
  end subroutine finish_line

end module types_const

! Local variables:
! folded-file: t
! End:

max_stars_repo_name: jirka-h/rnd-fortran

max_issues_repo_name: jirka-h/rnd-fortran

max_forks_repo_name: jirka-h/rnd-fortran


text:
Anderson Road is one of Davis main northsouth thoroughfares. It is named after JB Anderson John B. Anderson, the first president of the Board of Trustees, akin to the mayor, of Davis, and his brother, AG Anderson Gordon Anderson, who also served on the Board of Trustees, the precursor to the city council. Gordon Andersons descendants, his son Don Anderson and Dons daughter Jennifer Anderson, have continuously run Davis Lumber & Hardware Company, today known as Ace Hardware Davis Ace, for more than 100 years.

Anderson Road begins at its intersection with Russell Boulevard where it becomes La Rue Road upon entering Campus. Anderson runs North Davis north paralleling Highway 113 until it veers northeast, north of Peregrine Avenue and then turns due east, east of Tanager Avenue, ending where Anderson Road intersects F Street (aka County Road 101a) near the northern City Limits city limits. Bike lanes run along its length and there is a traffic circle at Alvarado Avenue north of the Covell Boulevard intersection.

 Gas Stations
   Arco AMPM Arco station at Russell Boulevard and Anderson (northwest corner of the intersection)
   Anderson Road Chevron Chevron station at Covell Boulevard and Anderson (southwest corner of the intersection
   Davis Shell Shell station at Covell Boulevard and Anderson (southeast corner of the intersection, connected to Anderson Plaza parking lot)

 Major Intersections
   Russell Boulevard
   8th Street West 8th Street
   Covell Boulevard
   Alvarado Avenue with Traffic Circle
   La Rue Road
   F Street

 Shopping
   University Mall
   Anderson Plaza
   AM PM

 Services
   Cesar Chavez School
   Congregation Bet Haverim
   Medical Vision Technology Optometry, ophthalmology & glasses
   Valley Oak Orthopaedics

 Houses
  1630 Anderson
  Bomb Shelter 720 Anderson


text:
      subroutine bjnser (x, l, jl, nl, ifl)

c-----------------------------------------------------------------------
c
c     subroutine: bjnser (x,l,jl,nl,ifl)
c
c     purpose:  to calculate the spherical bessel functions jl and nl
c
c     arguments:
c       x = argument of jl and nl
c       l = l value calculated (no offset)
c       jl = jl bessel function (abramowitz conventions)
c       nl = nl bessel function (abramowitz yl conventions)
c       ifl = 0 return both jl and nl
c             1 return jl only
c             2 return nl only
c
c     notes:  jl and nl are calculated by a series
c             expansion according to 10.1.2 and 10.1.3
c             in abramowitz and stegun (ninth printing),
c             page 437
c
c             error msgs written with PRINT statements.
c
c     first coded by r. c. albers on 26 jan 83
c
c     version 2
c
c     last modified: 27 jan 83 by r. c. albers
c
c-----------------------------------------------------------------------

      implicit double precision (a-h,o-z)

      complex*16 x,u,ux,del,pj,pn
      complex*16 jl,nl

      character*512 slog

      parameter (niter = 160, tol = 1.e-15)

      if (l .lt. 0) then
         call wlog(' l .lt. 0 in bjnser')
         stop 'bjnser 1'
      endif
      if (dble(x).lt. 0.) then
         write(slog,30) x
         call wlog(slog)
   30    format (' x = ', 1p, 2e14.6, ' is .le. 0 in bjnser')
         stop 'bjnser 2'
      endif

      lp1 = l+1
      u = x**2 / 2

c     make djl = 1 * 3 * 5 * ... * (2*l+1),
c          dnl = 1 * 3 * 5 * ... * (2*l-1)
      djl = 1
      fac = -1
      do 50 il = 1, lp1
         fac = fac + 2
         djl = fac * djl
   50 continue
      dnl = djl / (2*l+1)

      if (ifl .eq. 2)   goto 90
c     make jl
c     pj is term in { } in 10.1.2, del is last factor in the series
c     convergence test is (last factor)/(total term) <= tol
      pj = 1
      nf = 1
      nfac = 2*l + 3
      den = nfac
      sgn = -1
      ux = u
      do 60 il = 1, niter
         del = sgn*ux / den
         pj = pj + del
         trel = abs (del / pj)
         if (trel .le. tol)  goto 80
         sgn = -sgn
         ux = u*ux
         nf = nf+1
         nfac = nfac+2
         den = nf * nfac * den
   60 continue
      stop  'jl does not converge in bjnser'
   80 jl = pj * (x**l) / djl

   90 if (ifl.eq.1) return
c     make nl
c     pn is term in { } in 10.1.3, del is last factor in the series
c     convergence test is (last factor)/(total term) <= tol
      pn = 1
      nf = 1
      nfac = 1 - 2*l
      den = nfac
      sgn = -1
      ux = u
      do 100  il = 1, niter
         del = sgn * ux / den
         pn = pn + del
         trel = abs (del / pn)
         if (trel .le. tol) goto 120
         sgn = -sgn
         ux = u*ux
         nf = nf+1
         nfac = nfac+2
         den = nf * nfac * den
  100 continue
      stop  'nl does not converge in bjnser'
  120 nl = -pn * dnl / (x**lp1)

      return
      end


text:
!======================================================================!
SUBROUTINE UserCutLines_annulus()
    !----------------------------------------------------------------------!
    ! This program reads name.1D file created by NEU or GEN and averages
    ! the results in the homogeneous direction directions.
    ! The results are writen in files name_res.dat and name_res_plus.dat
    !----------------------------------------------------------------------!
    USE all_mod
    USE allp_mod
    USE les_mod
    USE pro_mod
    USE par_mod
    USE rans_mod
    !----------------------------------------------------------------------!
    IMPLICIT NONE
    !-----------------------------[Parameters]-----------------------------!
    REAL:: Ufric
    !------------------------------[Calling]-------------------------------!
    INTERFACE
        LOGICAL FUNCTION Approx(A,B,tol)
            REAL           :: A,B
            REAL, OPTIONAL :: tol
        END FUNCTION Approx
    END INTERFACE
    !-------------------------------[Locals]-------------------------------!
    INTEGER:: Nprob
    INTEGER:: pl
    INTEGER:: c
    INTEGER:: i
    INTEGER:: count
    INTEGER:: kk
    CHARACTER:: namCoo*80
    CHARACTER:: namPro*80
    CHARACTER:: namRes*80
    CHARACTER:: namRes_plus*80
    REAL,ALLOCATABLE    :: z_p(:)
    REAL,ALLOCATABLE    :: Ump(:)
    REAL,ALLOCATABLE    :: Vmp(:)
    REAL,ALLOCATABLE    :: Wmp(:)
    REAL,ALLOCATABLE    :: uup(:)
    REAL,ALLOCATABLE    :: vvp(:)
    REAL,ALLOCATABLE    :: wwp(:)
    REAL,ALLOCATABLE    :: uvp(:)
    REAL,ALLOCATABLE    :: uwp(:)
    REAL,ALLOCATABLE    :: vwp(:)
    REAL,ALLOCATABLE    :: Tmp(:)
    REAL,ALLOCATABLE    :: TTp(:)
    REAL,ALLOCATABLE    :: uTp(:)
    REAL,ALLOCATABLE    :: vTp(:)
    REAL,ALLOCATABLE    :: wTp(:)
    REAL,ALLOCATABLE    :: Ksgsp(:)
    REAL,ALLOCATABLE    :: ind(:)
    REAL,ALLOCATABLE    :: var_1(:)
    REAL,ALLOCATABLE    :: var_2(:)
    REAL,ALLOCATABLE    :: Wall_p(:)
    REAL,ALLOCATABLE    :: Ufric_p(:)
    INTEGER,ALLOCATABLE :: Np(:)
    INTEGER,ALLOCATABLE :: Ncount(:)
    REAL:: R
    REAL:: Lscale
    REAL:: R_max
    REAL:: b11
    REAL:: b22
    REAL:: b12
    REAL:: b21
    !--------------------------------[CVS]---------------------------------!
    !  $Id: UserProbe1D.f90,v 1.16 2002/11/25 10:33:17 niceno Exp $
    !  $Source: /home/muhamed/.CVSROOT/T-Rex/User/UserProbe1D.f90,v $
    !======================================================================!

    namPro = name

    namRes = name
    namRes(len_trim(name)+1:len_trim(name)+8) = '_res.dat'
    namRes_plus = name
    namRes_plus(len_trim(name)+1:len_trim(name)+13) = '_res_plus.dat'
    !>>>>>>>>>>>>>>>>>>>>>>!
    !     read 1D file     !
    !>>>>>>>>>>>>>>>>>>>>>>!
    namCoo = name
    namCoo(len_trim(name)+1:len_trim(name)+3) = '.1D'
    if(this < 2) write(6, *) '# Now reading the file:', namCoo
    open(9, FILE=namCoo)

    !---- write the number of searching intervals
    read(9,*) Nprob
    allocate(z_p(Nprob*2))
    allocate(ind(Nprob*2))

    !---- read the intervals positions
    do pl=1,Nprob
        read(9,*) ind(pl), z_p(pl)
    end do
    close(9)

    call SSORT (z_p, ind, Nprob, 0)

    call CalcShear(U % n, V % n, W % n, Shear)

    allocate(Np(Nprob));    Np=0
    allocate(Wall_p(Nprob));Wall_p=0.0
    allocate(Ump(Nprob));   Ump=0.0
    allocate(Vmp(Nprob));   Vmp=0.0
    allocate(Wmp(Nprob));   Wmp=0.0
    allocate(uup(Nprob));   uup=0.0
    allocate(vvp(Nprob));   vvp=0.0
    allocate(wwp(Nprob));   wwp=0.0
    allocate(uvp(Nprob));   uvp=0.0
    allocate(uwp(Nprob));   uwp=0.0
    allocate(vwp(Nprob));   vwp=0.0
    allocate(Ksgsp(Nprob)); Ksgsp=0.0
    allocate(var_1(Nprob));  var_1=0.0
    allocate(var_2(Nprob));  var_2=0.0
    allocate(Ufric_p(Nprob));  Ufric_p=0.0

    allocate(Ncount(Nprob)); Ncount=0
    count = 0

    allocate(Tmp(Nprob));   Tmp=0.0
    allocate(TTp(Nprob));   TTp=0.0
    allocate(uTp(Nprob));   uTp=0.0
    allocate(vTp(Nprob));   vTp=0.0
    allocate(wTp(Nprob));   wTp=0.0
    kk = 0

    Lscale = 0.0
    !+++++++++++++++++++++++++++++!
    !     average the results     !
    !+++++++++++++++++++++++++++++!
    !    R_max = 0.0

    !    do i = 1, Nprob
    !      R_max = max(z_p(i),R_max)
    !    end do
    !    open(12, FILE='inflow.dat')
    !    do c=1, NC
    !      if(zc(c) > 0.0.and.zc(c) < 0.16) then
    !        write(12,'(13E15.7)') xc(c), yc(c), U%n(c),V%n(c),W%n(c),uu%n(c),vv%n(c),ww%n(c),uv%n(c),uw%n(c),vw%n(c),Eps%n(c),f22%n(c)
    !      end if
    !    end do
    !    close(12)
    do i = 1, Nprob-1
        do c=1, NC
            Lscale = max(Lscale,WallDs(c))
            R = sqrt(xc(c)**2+yc(c)**2)
            if(R > abs(z_p(i+1)) .and. R < abs(z_p(i))) then
                if(zc(c) > 14.5.and.zc(c) < 15.0) then
                    Wall_p(i) = Wall_p(i) + R
                    b11 = xc(c)/R
                    b12 = yc(c)/R
                    b21 = -yc(c)/R
                    b22 = xc(c)/R
                    if(SIMULA==DNS) then

                        Ump(i)   = Ump(i) + U % n(c) * xc(c) / R  + V % n(c) * yc(c) / R
                        Vmp(i)   = Vmp(i) + (-U % n(c) * yc(c) / R  + V % n(c) * xc(c) / R)
                        !           Ump(i)   = Ump(i) + U % n(c)
                        !           Vmp(i)   = Vmp(i) + V % n(c)
                        Wmp(i)   = Wmp(i) + W % n(c)
                        uup(i)   = uup(i) + b11*b11*uu%n(c) + b11*b12*uv%n(c) &
                            + b12*b11*uv%n(c) + b12*b12*vv%n(c)
                        vvp(i)   = vvp(i) + b21*b21*uu%n(c) + b21*b22*uv%n(c) &
                            + b22*b21*uv%n(c) + b22*b22*vv%n(c)
                        wwp(i)   = wwp(i) + ww % n(c)
                        uvp(i)   = uvp(i) + b11*b21*uu%n(c) + b11*b22*uv%n(c) &
                            + b12*b21*uv%n(c) + b12*b22*vv%n(c)
                        uwp(i)   = uwp(i) + b11*uw%n(c) + b12*vw%n(c)
                        vwp(i)   = vwp(i) + b21*uw%n(c) + b22*vw%n(c)
                        var_1(i)   = var_1(i) + Kin % n(c)
                        var_2(i)   = var_2(i) + Eps % n(c)
                        if(IsNearWall(c)) then
                            Ufric_p(i) = Ufric_p(i) + (VISc * (U % n(c)**2 + V % n(c)**2 + W % n(c)**2)**0.5/WallDs(c))**0.5
                        end if

                    else if(SIMULA==LES .or. SIMULA == DNS) then
                        Ump(i)   = Ump(i) + U % mean(c) * xc(c) / R  + V % mean(c) * yc(c) / R
                        Vmp(i)   = Vmp(i) + (-U % mean(c) * yc(c) / R  + V % mean(c) * xc(c) / R)
                        Wmp(i)   = Wmp(i) + W % mean(c)
                        uup(i)   = uup(i) + (uu % mean(c)- (U % mean(c) * xc(c) / R  + V % mean(c) * yc(c) / R)**2.0)
                        vvp(i)   = vvp(i) + (vv % mean(c)- (-U % mean(c) * yc(c) / R  + V % mean(c) * xc(c) / R)**2.0)
                        wwp(i)   = wwp(i) + (ww % mean(c)- W % mean(c) * W % mean(c))
                        uvp(i)   = uvp(i) + (uv % mean(c)- (U % mean(c) * xc(c) / R  + V % mean(c) * yc(c) / R)* &
                            (-U % mean(c) * yc(c) / R  + V % mean(c) * xc(c) / R)   )
                        uwp(i)   = uwp(i) + (uw % mean(c)- (U % mean(c) * xc(c) / R  + V % mean(c) * yc(c) / R) * W % mean(c))
                        vwp(i)   = uwp(i) + (vw % mean(c)- (-U % mean(c) * yc(c) / R  + V % mean(c) * xc(c) / R) * W % mean(c))
                        if(IsNearWall(c)) then
                            Ufric_p(i) = Ufric_p(i) + (VISc * (U % mean(c)**2 + V % mean(c)**2 + W % mean(c)**2)**0.5/WallDs(c))**0.5
                        end if
                    end if

                    if(HOT==YES) then
                        if(SIMULA == LES .or. SIMULA == DNS) then
                            Tmp(i)   = Tmp(i) + T % mean(c)
                            TTp(i)   = TTp(i) + (TT % mean(c) - T % mean(c) * T % mean(c))
                            uTp(i)   = uTp(i) + (uT % mean(c) - (U % mean(c) * xc(c) / R  + V % mean(c) * yc(c) / R) * T % mean(c))
                            vTp(i)   = vTp(i) + (vT % mean(c) - (-U % mean(c) * yc(c) / R  + V % mean(c) * xc(c) / R) * T % mean(c))
                            wTp(i)   = wTp(i) + (wT % mean(c) - w % mean(c) * T % mean(c))
                        else
                            Tmp(i)   = Tmp(i) + T % n(c)
                        end if
                    end if
                    Ncount(i) = Ncount(i) + 1
                end if
            end if
        end do
    end do

    !---- average over all processors
    do pl=1, Nprob-1
        call IGlSum(Ncount(pl))

        call GloSum(Wall_p(pl))

        call GloSum(Ump(pl))
        call GloSum(Vmp(pl))
        call GloSum(Wmp(pl))

        call GloSum(uup(pl))
        call GloSum(vvp(pl))
        call GloSum(wwp(pl))

        call GloSum(uvp(pl))
        call GloSum(uwp(pl))
        call GloSum(vwp(pl))
        call GloSum(var_1(pl))
        call GloSum(var_2(pl))
        call GloSum(Ufric_p(pl))

        count =  count + Ncount(pl)

        if(HOT==YES) then
            call GloSum(Tmp(pl))
            call GloSum(TTp(pl))
            call GloSum(uTp(pl))
            call GloSum(vTp(pl))
            call GloSum(wTp(pl))
        end if
    end do

    do i = 1, Nprob-1
        if(Ncount(i) /= 0) then
            Wall_p(i) =  Wall_p(i)/Ncount(i)
            Ump(i)    =  Ump(i)/Ncount(i)
            Vmp(i)    =  Vmp(i)/Ncount(i)
            Wmp(i)    =  Wmp(i)/Ncount(i)
            uup(i)    =  uup(i)/Ncount(i)
            vvp(i)    =  vvp(i)/Ncount(i)
            wwp(i)    =  wwp(i)/Ncount(i)
            uvp(i)    =  uvp(i)/Ncount(i)
            uwp(i)    =  uwp(i)/Ncount(i)
            vwp(i)    =  vwp(i)/Ncount(i)
            var_1(i)  =  var_1(i)/Ncount(i)
            var_2(i)  =  var_2(i)/Ncount(i)
            Ufric_p(i)  =  Ufric_p(i)/Ncount(i)
            if(HOT == YES) then
                Tmp(i)    =  Tmp(i)/Ncount(i)
                TTp(i)    =  TTp(i)/Ncount(i)
                uTp(i)    =  uTp(i)/Ncount(i)
                vTp(i)    =  vTp(i)/Ncount(i)
                wTp(i)    =  wTp(i)/Ncount(i)
            end if
        end if
    end do

    Ufric = 0.0
    do i = 1, Nprob
        Ufric = max(Ufric_p(i),Ufric)
    end do

continue

open(3,FILE=namRes)
write(3,'(A1,2(A10, F13.5))') '#', 'Utau = ', Ufric, 'Re_tau = ', Ufric*R_max/VISc
if(SIMULA == LES .or. SIMULA == DNS) then
    write(3,'(A1,2X,A80)') '#', '1:Xrad, 2:U, 3:V, 4:W, 5:uu, 6:vv, 7:ww, 8:uv, 9:uw, 10:vw, 11:Kin'
    do i = 1, Nprob-1
        if(Ncount(i) /= 0) then
            write(3,'(11E15.7)') Wall_p(i), Ump(i), Vmp(i), Wmp(i),   &
                uup(i), vvp(i), wwp(i), uvp(i), uwp(i),      &
                vwp(i), 0.5*(uup(i) + vvp(i) + wwp(i))
        end if
    end do
else if(SIMULA == DNS) then
    write(3,'(A1,1X,A50)') '#', '1:Xrad, 2:U, 3:V, 4:W, 5:Kin, 6:Eps, 7:zeta, 8:f'
    do i = 1, Nprob-1
        if(Ncount(i) /= 0) then
            write(3,'(12E15.7)') Wall_p(i), Ump(i), Vmp(i), Wmp(i),   &
                uup(i), vvp(i), wwp(i), uvp(i), uwp(i), vwp(i), var_1(i), var_2(i)
        end if
    end do
end if
close(3)

open(3,FILE=namRes_plus)
write(3,'(A1,2(A10, F13.5))') '#', 'Utau = ', Ufric, 'Re_tau = ', Ufric*R_max/VISc
if(SIMULA == LES .or. SIMULA==DNS) then
    write(3,'(A1,2X,A80)') '#', '1:Xrad+, 2:U+, 3:V+, 4:W+, 5:uu+, 6:vv+, 7:ww+, 8:uv+, 9:uw+, 10:vw+, 11:Kin+'
    do i = 1, Nprob-1
        if(Ncount(i) /= 0) then
            write(3,'(11E15.7)') Wall_p(i)*Ufric/VISc, abs(Ump(i))/Ufric, abs(Vmp(i))/Ufric, abs(Wmp(i))/Ufric,   &
                uup(i)/Ufric**2.0, vvp(i)/Ufric**2.0, wwp(i)/Ufric**2.0, uvp(i)/Ufric**2.0, uwp(i)/Ufric**2.0,      &
                vwp(i)/Ufric**2.0, (uup(i) + vvp(i) + wwp(i))/Ufric**2.0
        end if
    end do
else if(SIMULA == DNS) then
    write(3,'(A1,1X,A80)') '#', '1:Xrad+, 2:U+, 3:V+, 4:W+, 5:Kin+, 6:Eps+, 7:zeta+, 8:f+                     '
    do i = 1, Nprob-1
        if(Ncount(i) /= 0) then
            write(3,'(11E15.7)') (Wall_p(i))*Ufric/VISc, Ump(i)/Ufric, Vmp(i)/Ufric, Wmp(i)/Ufric,   &
                uup(i)/Ufric**2.0, vvp(i)/Ufric**2.0, wwp(i)/Ufric**2.0, uvp(i)/Ufric**2.0, uwp(i)/Ufric**2.0, &
                vwp(i)/Ufric**2.0, 0.5*(uup(i)/Ufric**2.0+vvp(i)/Ufric**2.0+wwp(i)/Ufric**2.0)
        end if
    end do
end if
close(3)

deallocate(Np)
deallocate(z_p)
deallocate(Ump)
deallocate(Vmp)
deallocate(Wmp)
deallocate(uup)
deallocate(vvp)
deallocate(wwp)
deallocate(uvp)
deallocate(uwp)
deallocate(vwp)
deallocate(var_1)
deallocate(var_2)
deallocate(Ksgsp)
deallocate(Ufric_p)
deallocate(Wall_p)
if(HOT==YES) then
    deallocate(Tmp)
    deallocate(TTp)
    deallocate(uTp)
    deallocate(vTp)
    deallocate(wTp)
end if
END SUBROUTINE UserCutLines_annulus

max_stars_repo_name: palkinev/T-FlowS-old-compressible

max_issues_repo_name: palkinev/T-FlowS-old-compressible

max_forks_repo_name: palkinev/T-FlowS-old-compressible


text: #redirect Shelly Bailes and Ellen Pontac


! Mikayla Webber
! 4670 Numerical Analysis
! Homework Four Due 11/6/17

module secret

integer :: fcounter

end module secret

!-----------------------------------------------------------------------------------------------------

program WebberHomework4Question3
use secret
implicit none

stop
end program WebberHomework4Question3

!-----------------------------------------------------------------------------------------------------

implicit none

	double precision :: a,b,ci,d,f,simp2d
	integer :: m,n
	external f

	a = 0.0d0
	b = 2.0d0
	ci = 0.0d0
	d = 1.0d0
	simp2d = 0.0d0

	m = 20
	n = 20

	call CompSimpson2d(a,b,ci,d,m,n,f,simp2d)

	print*,"simpson 2d = ",simp2d

	stop
	end

	double precision function f(x,y)
	implicit none

	double precision :: x,y

	f = x*y**2.0d0

	return
	end


! { dg-do run }

  call test
contains
  subroutine check (x, y, l)
    integer :: x, y
    logical :: l
    l = l .or. x .ne. y
  end subroutine check

  subroutine foo (c, d, e, f, g, h, i, j, k, n)
    use omp_lib
    integer :: n
    character (len = *) :: c
    character (len = n) :: d
    integer, dimension (2, 3:5, n) :: e
    integer, dimension (2, 3:n, n) :: f
    character (len = *), dimension (5, 3:n) :: g
    character (len = n), dimension (5, 3:n) :: h
    real, dimension (:, :, :) :: i
    double precision, dimension (3:, 5:, 7:) :: j
    integer, dimension (:, :, :) :: k
    logical :: l
    integer :: p, q, r
    character (len = n) :: s
    integer, dimension (2, 3:5, n) :: t
    integer, dimension (2, 3:n, n) :: u
    character (len = n), dimension (5, 3:n) :: v
    character (len = 2 * n + 24) :: w
    integer :: x
    character (len = 1) :: y
    s = 'PQRSTUV'
    forall (p = 1:2, q = 3:5, r = 1:7) t(p, q, r) = -10 + p - q + 2 * r
    forall (p = 1:2, q = 3:7, r = 1:7) u(p, q, r) = 30 - p + q - 2 * r
    forall (p = 1:5, q = 3:7, p + q .le. 8) v(p, q) = '_+|/Oo_'
    forall (p = 1:5, q = 3:7, p + q .gt. 8) v(p, q) = '///|||!'
    l = .false.
!$omp parallel default (none) firstprivate (c, d, e, f, g, h, i, j, k) &
!$omp & firstprivate (s, t, u, v) reduction (.or.:l) num_threads (6) &
!$omp private (p, q, r, w, x, y)
    l = l .or. c .ne. 'abcdefghijkl'
    l = l .or. d .ne. 'ABCDEFG'
    l = l .or. s .ne. 'PQRSTUV'
    do 100, p = 1, 2
      do 100, q = 3, 7
	do 100, r = 1, 7
	  if (q .lt. 6) l = l .or. e(p, q, r) .ne. 5 + p + q + 2 * r
	  l = l .or. f(p, q, r) .ne. 25 + p + q + 2 * r
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. g(r, q) .ne. '0123456789AB'
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. g(r, q) .ne. '9876543210ZY'
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. h(r, q) .ne. '0123456'
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. h(r, q) .ne. '9876543'
	  if (q .lt. 6) l = l .or. t(p, q, r) .ne. -10 + p - q + 2 * r
	  l = l .or. u(p, q, r) .ne. 30 - p + q - 2 * r
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. v(r, q) .ne. '_+|/Oo_'
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. v(r, q) .ne. '///|||!'
100 continue
    do 101, p = 3, 5
      do 101, q = 2, 6
	do 101, r = 1, 7
	  l = l .or. i(p - 2, q - 1, r) .ne. 7.5 * p * q * r
	  l = l .or. j(p, q + 3, r + 6) .ne. 9.5 * p * q * r
101 continue
    do 102, p = 1, 5
      do 102, q = 4, 6
	l = l .or. k(p, 1, q - 3) .ne. 19 + p + 7 + 3 * q
102 continue
    x = omp_get_thread_num ()
    w = ''
    if (x .eq. 0) w = 'thread0thr_number_0THREAD0THR_NUMBER_0'
    if (x .eq. 1) w = 'thread1thr_number_1THREAD1THR_NUMBER_1'
    if (x .eq. 2) w = 'thread2thr_number_2THREAD2THR_NUMBER_2'
    if (x .eq. 3) w = 'thread3thr_number_3THREAD3THR_NUMBER_3'
    if (x .eq. 4) w = 'thread4thr_number_4THREAD4THR_NUMBER_4'
    if (x .eq. 5) w = 'thread5thr_number_5THREAD5THR_NUMBER_5'
    c = w(8:19)
    d = w(1:7)
    forall (p = 1:2, q = 3:5, r = 1:7) e(p, q, r) = 5 * x + p + q + 2 * r
    forall (p = 1:2, q = 3:7, r = 1:7) f(p, q, r) = 25 * x + p + q + 2 * r
    forall (p = 1:5, q = 3:7, p + q .le. 8) g(p, q) = w(8:19)
    forall (p = 1:5, q = 3:7, p + q .gt. 8) g(p, q) = w(27:38)
    forall (p = 1:5, q = 3:7, p + q .le. 8) h(p, q) = w(1:7)
    forall (p = 1:5, q = 3:7, p + q .gt. 8) h(p, q) = w(20:26)
    forall (p = 3:5, q = 2:6, r = 1:7) i(p - 2, q - 1, r) = (7.5 + x) * p * q * r
    forall (p = 3:5, q = 2:6, r = 1:7) j(p, q + 3, r + 6) = (9.5 + x) * p * q * r
    forall (p = 1:5, q = 7:7, r = 4:6) k(p, q - 6, r - 3) = 19 + x + p + q + 3 * r
    s = w(20:26)
    forall (p = 1:2, q = 3:5, r = 1:7) t(p, q, r) = -10 + x + p - q + 2 * r
    forall (p = 1:2, q = 3:7, r = 1:7) u(p, q, r) = 30 - x - p + q - 2 * r
    forall (p = 1:5, q = 3:7, p + q .le. 8) v(p, q) = w(1:7)
    forall (p = 1:5, q = 3:7, p + q .gt. 8) v(p, q) = w(20:26)
!$omp barrier
    y = ''
    if (x .eq. 0) y = '0'
    if (x .eq. 1) y = '1'
    if (x .eq. 2) y = '2'
    if (x .eq. 3) y = '3'
    if (x .eq. 4) y = '4'
    if (x .eq. 5) y = '5'
    l = l .or. w(7:7) .ne. y
    l = l .or. w(19:19) .ne. y
    l = l .or. w(26:26) .ne. y
    l = l .or. w(38:38) .ne. y
    l = l .or. c .ne. w(8:19)
    l = l .or. d .ne. w(1:7)
    l = l .or. s .ne. w(20:26)
    do 103, p = 1, 2
      do 103, q = 3, 7
	do 103, r = 1, 7
	  if (q .lt. 6) l = l .or. e(p, q, r) .ne. 5 * x + p + q + 2 * r
	  l = l .or. f(p, q, r) .ne. 25 * x + p + q + 2 * r
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. g(r, q) .ne. w(8:19)
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. g(r, q) .ne. w(27:38)
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. h(r, q) .ne. w(1:7)
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. h(r, q) .ne. w(20:26)
	  if (q .lt. 6) l = l .or. t(p, q, r) .ne. -10 + x + p - q + 2 * r
	  l = l .or. u(p, q, r) .ne. 30 - x - p + q - 2 * r
	  if (r .lt. 6 .and. q + r .le. 8) l = l .or. v(r, q) .ne. w(1:7)
	  if (r .lt. 6 .and. q + r .gt. 8) l = l .or. v(r, q) .ne. w(20:26)
103 continue
    do 104, p = 3, 5
      do 104, q = 2, 6
	do 104, r = 1, 7
	  l = l .or. i(p - 2, q - 1, r) .ne. (7.5 + x) * p * q * r
	  l = l .or. j(p, q + 3, r + 6) .ne. (9.5 + x) * p * q * r
104 continue
    do 105, p = 1, 5
      do 105, q = 4, 6
	l = l .or. k(p, 1, q - 3) .ne. 19 + x + p + 7 + 3 * q
105 continue
    call check (size (e, 1), 2, l)
    call check (size (e, 2), 3, l)
    call check (size (e, 3), 7, l)
    call check (size (e), 42, l)
    call check (size (f, 1), 2, l)
    call check (size (f, 2), 5, l)
    call check (size (f, 3), 7, l)
    call check (size (f), 70, l)
    call check (size (g, 1), 5, l)
    call check (size (g, 2), 5, l)
    call check (size (g), 25, l)
    call check (size (h, 1), 5, l)
    call check (size (h, 2), 5, l)
    call check (size (h), 25, l)
    call check (size (i, 1), 3, l)
    call check (size (i, 2), 5, l)
    call check (size (i, 3), 7, l)
    call check (size (i), 105, l)
    call check (size (j, 1), 4, l)
    call check (size (j, 2), 5, l)
    call check (size (j, 3), 7, l)
    call check (size (j), 140, l)
    call check (size (k, 1), 5, l)
    call check (size (k, 2), 1, l)
    call check (size (k, 3), 3, l)
    call check (size (k), 15, l)
!$omp end parallel
    if (l) call abort
  end subroutine foo

  subroutine test
    character (len = 12) :: c
    character (len = 7) :: d
    integer, dimension (2, 3:5, 7) :: e
    integer, dimension (2, 3:7, 7) :: f
    character (len = 12), dimension (5, 3:7) :: g
    character (len = 7), dimension (5, 3:7) :: h
    real, dimension (3:5, 2:6, 1:7) :: i
    double precision, dimension (3:6, 2:6, 1:7) :: j
    integer, dimension (1:5, 7:7, 4:6) :: k
    integer :: p, q, r
    c = 'abcdefghijkl'
    d = 'ABCDEFG'
    forall (p = 1:2, q = 3:5, r = 1:7) e(p, q, r) = 5 + p + q + 2 * r
    forall (p = 1:2, q = 3:7, r = 1:7) f(p, q, r) = 25 + p + q + 2 * r
    forall (p = 1:5, q = 3:7, p + q .le. 8) g(p, q) = '0123456789AB'
    forall (p = 1:5, q = 3:7, p + q .gt. 8) g(p, q) = '9876543210ZY'
    forall (p = 1:5, q = 3:7, p + q .le. 8) h(p, q) = '0123456'
    forall (p = 1:5, q = 3:7, p + q .gt. 8) h(p, q) = '9876543'
    forall (p = 3:5, q = 2:6, r = 1:7) i(p, q, r) = 7.5 * p * q * r
    forall (p = 3:6, q = 2:6, r = 1:7) j(p, q, r) = 9.5 * p * q * r
    forall (p = 1:5, q = 7:7, r = 4:6) k(p, q, r) = 19 + p + q + 3 * r
    call foo (c, d, e, f, g, h, i, j, k, 7)
  end subroutine test
end