text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C
c=======================================================================
      PROGRAM BLOTII
C=======================================================================

C                         *** BLOT ***

C   --*** BLOT *** (BLOT) Post-processing plot program
C   --
C   --BLOT is a graphics program for post-processing of finite element
C   --analyses output in the EXODUS II database format.  BLOT combines the
C   --plotting capabilities of DETOUR, TPLOT, and SPLOT.  It is command
C   --driven with free-format input.  BLOT can drive any graphics device
C   --supported by the Sandia Virtual Device Interface.
C   --
C   --DETOUR produces mesh plots with various representations of the
C   --analysis output variables. The major capabilities of DETOUR are
C   --deformed mesh plots, line contours, filled (painted) contours,
C   --vector plots of two variables (e.g., velocity vectors), and symbol
C   --plots of scalar variables (e.g., discrete cracks). DETOURs features
C   --include element selection by element block, element birth and death,
C   --multiple views for combining several displays on each plot, and
C   --symmetry mirroring.
C   --
C   --TPLOT generates time-versus-variable plots or variable-versus-variable
C   --plots where the variables are user-selected global, nodal or element
C   --EXODUS database variables.  TPLOTs features include multiple curve
C   --plots and neutral file output.
C   --
C   --SPLOT generates distance-versus-variable plots at selected time steps
C   --where the distance is the accumulated distance between pairs of nodes
C   --or element centers and the variable is a nodal or element EXODUS
C   --database variable.  SPLOTs features include multiple curve plots
C   --and neutral file output.
C   --
C   --Plot times from the database are processed in one of the following modes:
C   --   Uniform time interval - plot at every DELT time interval from
C   --      TMIN to TMAX.
C   --   All available times - plot at all available times from TMIN
C   --      to TMAX.
C   --   User selected times - plot at user selected times.
C   --The closest available time step is chosen.
C   --
C   --Expected input:
C   --   o The commands on the standard input device.
C   --   o The input EXODUS database on unit 11.
C   --
C   --Output:
C   --   o A listing of the input database information and any errors
C   --     found on the standard output device.
C   --   o The plots on the specified graphics device.
C   --   o A GRAFAID neutral file on unit 20.
C   --   o A GROPE listing file on unit 21.

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*                    ISSUED BY SANDIA LABORATORIES,                   *
*                      A PRIME CONTRACTOR TO THE                      *
*                  UNITED STATES DEPARTMENT OF ENERGY                 *
* * * * * * * * * * * * * *   N O T I C E   * * * * * * * * * * * * * *
* This program was prepared as an account of work sponsored by the    *
* United States Government.  Neither the United States nor the United *
* States Department of Energy nor any of their employees, nor any of  *
* their contractors, subcontractors, or 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 or process disclosed, or represents *
* that its use would not infringe privately owned rights.             *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

C   --Developed at Sandia National Laboratories.
C   --
C   --Current author and code sponsor: John Glick
C   --
C   --Revision History:
C   --   05/88  Added PATHLINE (Amy Gilkey)
C   --   03/88  Added Master/Slave logic (Amy Gilkey)
C   --   10/87  Added GROPE (Amy Gilkey)
C   --   10/87  Converted from SEACO to EXODUS database (Amy Gilkey)
C   --   07/87  Combined DETOUR, TPLOT, and SPLOT (Amy Gilkey)
C   --DETOUR:
C   --   05/86  Started to add 3D (Amy Gilkey)
C   --   03/85  New sponsor (Amy Gilkey)
C   --   11/82  Created (Dennis Flanagan)
C   --TPLOT:
C   --   xx/xx  Neutral file code (Greg Sjaardema)
C   --   01/86  New sponsor (Amy Gilkey)
C   --   xx/xx  Modified for the VAX11/780 (Johnny Biffle)
C   --   08/80  Created (Zelma Beisinger)
C   --SPLOT:
C   --   xx/xx  Neutral file code (Greg Sjaardema)
C   --   01/86  New sponsor (Amy Gilkey)
C   --   xx/xx  Created (Mary Sagartz)
C   --
C   --Source is in FORTRAN 77
C   --
C   --External software used:
C   --   SVDI graphics package
C   --   PLT graphics package
C   --   SUPES package (dynamic memory, free-field reader, FORTRAN extensions)
C   --
C   --Runs on VAX VMS, sun, Unicos, Ultrix
C   --
C   --Documentation: SAND86-0914, printed July 1987                       f
C   --   "DETOUR - A Deformed Mesh / Contour Plot Program"
C   --Documentation: SAND86-0883, printed August 1986
C   --   "TPLOT - A Time History or X-Y Plot Program for the Output
C   --   of a Finite Element Analysis"
C   --Documentation: SAND86-0882, printed August 1986
C   --   "SPLOT - A Distances-versus-Variable Plot Program for the Output
C   --   of a Finite Element Analysis"

C   --NOTE: All parameters and common areas are defined in BLKDAT

      include 'exodusII.inc'

C      --These parameters define the indices of 2D and 3D limit arrays
      PARAMETER (KLFT=1, KRGT=2, KBOT=3, KTOP=4, KNEA=5, KFAR=6)

C      --These parameters define the mesh display (see MSHLIN of /MSHOPT/)
      PARAMETER (MSHNON=0, MSHBOR=1, MSHDIV=2, MSHSEL=3, MSHALL=4)

      PARAMETER (MXQARC=MXSTLN)
      include 'progqa.blk'
      include 'dbase.blk'
      include 'dbname.blk'
      include 'dbtitl.blk'
      include 'dbnums.blk'
      include 'dbnumgq.blk'
      include 'd3nums.blk'
      include 'dbnams.blk'
      include 'times.blk'
      include 'layout.blk'
      include 'deform.blk'
      include 'mshlim.blk'

      include 'selne.blk'
      include 'neutr.blk'
      include 'csv.blk'
      include 'outfil.blk'
      include 'sizes.blk'

      include 'light.blk'
      include 'icrnbw.blk'
      include 'legopt.blk'

      include 'argparse.inc'

      common /debugc/ cdebug
      common /debugn/ idebug
      character*8 cdebug

C ... Used for hardcopy device drivers...
      common /blotans/ BLTANS
      character*2 BLTANS

      DIMENSION A(1),IA(1)
      EQUIVALENCE (A(1),IA(1))
      CHARACTER*1 C(1)
C      --A - the dynamic memory base array

      character*2048 scratch
      character*256  option, value

      LOGICAL MESHOK, DTOK, LNOK, SPOK, TPOK
      LOGICAL MAPND, MAPEL
      CHARACTER*(MXSTLN) CURPRO
      INTEGER NEWELB
      CHARACTER NEWPRO

      INTEGER IDUM
      REAL RDUM
      CHARACTER*8 CDUM
C     external sample_handler
      EXTERNAL TPREAD
      EXTERNAL BLKDAT

C     --The compute word size and I/O word size
      INTEGER CMPSIZ, IOWS

C ... Initialize Rainbow and Light (ICRNBW and LIGHT common blocks)
C     Block Data doesn't seem to work reliably on all systems

C ... Lights are stored as x, y, z, brightness, Vector is normalized in
C     shade.f.  NLIT is the number of lights.
      RMULT = 1.0
      GMULT = 1.0
      BMULT = 1.0
      LITE(1,1) = -1.0
      LITE(2,1) =  1.0
      LITE(3,1) =  1.0
      LITE(4,1) =  1.0
      LITE(5,1) = -0.57735
      LITE(6,1) =  0.57735
      LITE(7,1) =  0.57735
      LITE(8,1) =  1.0
      NLIT = 1
      AMBIENT = 0.2

      cdebug = '        '
      idebug = 0
      CMPSIZ = 0
      IOWS   = 0
      NEUTRL = 0

      NDB = 11
      NEU = 20
      NEUOPN = .FALSE.
      NEUGRF = 22
      GRFOPN = .FALSE.
      NCSV = 23
      CSVOPN = .FALSE.
      call version(qainfo)

C     Initialize variables
      VERS = 0.0
      IDUM = 0

      CALL STRTUP (QAINFO)

      draw(1) = '                                '
      draw(2) = '                                '
      draw(3) = '                                '
      draw(4) = '                                '
      DRAW(1) = QAINFO(1)
      DRAW(2) = '  ' // QAINFO(3)(3:8)
      DRAW(3)(1:8) = QAINFO(5)(1:8)
      DRAW(4)(1:8) = QAINFO(6)(1:8)

      WRITE (*, 10020)
      CALL BANNER (0, QAINFO,
     &  'A DEFORMED MESH / CONTOUR PLOT PROGRAM',
     &  'WITH X-Y PLOTTING CAPABILITIES',
     &  'FOR POST-PROCESSING OF FINITE ELEMENT ANALYSES'
     &  )
      CALL CPYRGT(0, '2009')

C   --Initialize dynamic memory

      CALL MDINIT (A)
      CALL MCINIT (C)
      CALL MDFILL (-999)
      CALL MCFILL ('Z')
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

C   --Open database file

      NARG = argument_count()
      if (narg .eq. 0) then
        CALL PRTERR ('FATAL', 'Filename not specified.')
        CALL PRTERR ('CMDSPEC',
     *    'Syntax is: "blot.dev [-basename basename] [-ps_option num]'//
     *    ' [-nomap node|element|all] filename"')
        GOTO 170
      end if

      CALL get_argument(narg,dbname, lfil)
      dbname(lfil+1:) = ' '
      ndb = exopen(dbname(:lfil), EXREAD, CMPSIZ, IOWS, vers, ierr)
      IF (IERR .NE. 0) THEN
        SCRATCH = 'Database "'//dbname(:lfil)//'" does not exist.'
        CALL PRTERR ('FATAL', SCRATCH(:LENSTR(SCRATCH)))
        CALL PRTERR ('CMDSPEC',
     *    'Syntax is: "blot.dev [-basename basename] [-ps_option num]'//
     *    ' [-nomap node|element|all] [-show_filename] filename"')
        GOTO 170
      END IF
      EXODUS = .FALSE.

      CALL INISTR (3, ' ', CAPTN(1,1))
      CALL INISTR (3, ' ', CAPTN(1,2))
      call exinq(ndb, EXDBMXUSNM, namlen, rdum, cdum, ierr)
      call exmxnm(ndb, namlen, ierr)

C ... Get basename of the database file to use for csv, neu and other output files
C     Assume first that the basename is from the database name.  It may be replaced
C     later if the user added a -hardcopy or -basename argument.
      last = indexr(dbname(:lfil), '.')
      if (last .gt. 2) then
        basenam = dbname(:last-1)
      else
        basenam = dbname
      end if

C ... Now see if a command-line option overrides this.
C     Options should all be of the form "-option arg"
      bltans = '7'

C ... By default, map both nodes and elements
      mapel = .true.
      mapnd = .true.

      if (narg .gt. 1) then
        i = 1
        do
          CALL get_argument(i,option, lo)
          i = i + 1
          if (option(:lo) .eq. '--show_filename' .or.
     *      option(:lo) .eq. '-show_filename') then
            captn(3,1) = dbname(:lfil)
            captn(3,2) = dbname(:lfil)

          else if (option(:lo) .eq. '-hardcopy' .or.
     *      option(:lo) .eq. '--hardcopy' .or.
     *      option(:lo) .eq. '-basename' .or.
     *      option(:lo) .eq. '--basename') then
            CALL get_argument(i,value,  lv)
            i = i + 1
            basenam = value(:lv)

          else if (option(:lo) .eq. '-ps_option' .or.
     *        option(:lo) .eq. '--ps_option') then
            CALL get_argument(i,value,  lv)
            i = i + 1
            if (lv .le. 2) then
              bltans = value(:lv)
            end if

          else if (option(:lo) .eq. '-nomap' .or.
     *        option(:lo) .eq. '--nomap') then
            CALL get_argument(i,value,  lv)
            i = i + 1
            if (value(1:1) .eq. 'n' .or. value(1:1) .eq. 'N')
     *        mapnd = .false.
            if (value(1:1) .eq. 'e' .or. value(1:1) .eq. 'E')
     *        mapel = .false.
            if (value(1:1) .eq. 'a' .or. value(1:1) .eq. 'A') then
              mapnd = .false.
              mapel = .false.
            end if
          end if
          if (i .gt. narg) exit
        end do
      end if

C   --Set error reporting level
      CALL EXOPTS(EXABRT,IERR)

C   --Initialize graphics
      CALL GRINIT (DBORD0, CHLSIZ)

C   --Read and print database header

      call exgini(ndb, title,
     *  ndim, numnp, numel, nelblk,
     *  numnps, numess, ierr)
      if (numnps .gt. 0) then
        call exinq(ndb, EXNSNL, lnpsnl, rdum, cdum, ierr)
        call exinq(ndb, EXNSDF, lnpsdf, rdum, cdum, ierr)
      else
        lnpsnl = 0
        lnpsdf = 0
      end if
      if (numess .gt. 0) then
        call exinq(ndb, EXSSNL, lessnl, rdum, cdum, ierr)
        call exinq(ndb, EXSSEL, lessel, rdum, cdum, ierr)
        call exinq(ndb, EXSSDF, lessdf, rdum, cdum, ierr)
      else
        lessnl = 0
        lessel = 0
        lessdf = 0
      end if

      NUMNPF = NUMNP

      CALL PRINIT ('NTIS', -1, NDB, DBNAME, TITLE,
     &  NDIM, NUMNP, NUMEL, NELBLK,
     &  NUMNPS, LNPSNL, LNPSDF, NUMESS, LESSEL, LESSNL,
     &  NVARGL, NVARNP, NVAREL, NVARNS, NVARSS)

C   --Check whether mesh may be plotted

      CALL MSCHK (.FALSE., MESHOK)
      IF (.NOT. MESHOK)
     &  CALL PRTERR ('WARNING', 'No mesh is defined')

      IS3DIM = (NDIM .GE. 3)

C   --Read coordinates (full set)

C   --SCALER uses MDFIND to find XN, YN, ZN
      CALL MDGET (NDIM * NUMNP)
      CALL MDRSRV ('XN', KXN, NUMNP)
      CALL MDRSRV ('YN', KYN, NUMNP)
      IF (IS3DIM) THEN
        CALL MDRSRV ('ZN', KZN, NUMNP)
      ELSE
        KZN = 1
      END IF
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      call exgcor(ndb, a(kxn), a(kyn), a(kzn), ierr)

C     Element Block IDs
      CALL MDRSRV ('IDELB', KIDELB, NELBLK)
C     Element Block IDs of HEXSHELL element blocks
      CALL MDRSRV ('HEXID', KHEXID, NELBLK)
C     Element Block types
      CALL MCRSRV ('NAMELB', KNMLB, NELBLK * MXSTLN)

C     Element Block names
      CALL MCRSRV ('EBNAME', KNMEB, NELBLK * NAMLEN)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      CALL EXGEBI (NDB, IA(KIDELB), IERR)
      if (nelblk .gt. 0) then
        CALL MDLONG ('IDELB', KIDELB, NELBLK)
        CALL MCLONG ('NAMELB', KNMLB, NELBLK * MXSTLN)
        CALL MCLONG ('EBNAME', KNMEB, NELBLK * NAMLEN)
      end if
C     Number of elements in element block
      CALL MDRSRV ('NUMELB', KNELB, NELBLK)
C     Number of nodes per elements in element block
      CALL MDRSRV ('NUMLNK', KNLNKE, NELBLK)
C     Number of attributes in element block
      CALL MDRSRV ('NUMATR', KNATR, NELBLK)
C     Connectivity array
      CALL MDRSRV ('LINK', KLINKE, 0)
C     Connectivity pointer array
      CALL MDRSRV ('LPTR', KLPTR, NELBLK)
C     Attribute array
      CALL MDRSRV ('ATRIB', KATRIB, 0)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

C     Initialize element block integer arrays
      CALL INIINT(NELBLK, 0, IA(KIDELB))
      CALL INIINT(NELBLK, 0, IA(KNELB))
      CALL INIINT(NELBLK, 0, IA(KNLNKE))
      CALL INIINT(NELBLK, 0, IA(KNATR))

      CALL INISTR (NELBLK, ' ', C(KNMLB))

      CALL EXGEBI (NDB, IA(KIDELB), IERR)
C   --Read element block connectivity
      CALL DBIELB (NDB, '*', 1, NELBLK, IA(KIDELB), IA(KNELB),
     &  IA(KNLNKE), IA(KNATR), A, IA, KLINKE, KATRIB, C(KNMLB),
     &  C(KNMEB), IA(KLPTR), NAMLEN, *170)

C     Count the number of element blocks that contain HEXSHELLs. Store
C     the element block ID of HEXSHELL element blocks.
      ISHEX = 0
      CALL RDTYPE (NELBLK, C(KNMLB), IA(KIDELB), IA(KNELB),
     $  ISHEX, IA(KHEXID), NSHL)

      if (ishex .gt. 0) then
        nebsiz = nelblk + ishex
        CALL MDLONG ('IDELB', KIDELB, NEBSIZ)
        CALL MCLONG ('NAMELB', KNMLB, NEBSIZ * MXSTLN)
        CALL MCLONG ('EBNAME', KNMEB, NEBSIZ * NAMLEN)
C     Number of elements in element block
        CALL MDLONG ('NUMELB', KNELB, NEBSIZ)
C     Number of nodes per elements in element block
        CALL MDLONG ('NUMLNK', KNLNKE, NEBSIZ)
C     Number of attributes in element block
        CALL MDLONG ('NUMATR', KNATR, NEBSIZ)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

C     Initialize element block integer arrays
        CALL INIINT(ISHEX, 0, IA(KIDELB+NELBLK))
        CALL INIINT(ISHEX, 0, IA(KNELB+NELBLK))
        CALL INIINT(ISHEX, 0, IA(KNLNKE+NELBLK))
        CALL INIINT(ISHEX, 0, IA(KNATR+NELBLK))

C     Check for HEXSHELL - split HEXSHELL element block into a 'HEX'
C     element block and a 'SHELL' element block.
        CALL MDRSRV ('IDSCR', KIDSCR, NELBLK)
        do 350 j = 1, nelblk
          IA(KIDSCR+j-1) = IA(KIDELB+j-1)
 350    continue

        CALL PROCHS(A, IA, NELBLK, IA(KIDELB), IA(KIDSCR), IA(KNELB),
     &    IA(KNLNKE), IA(KNATR), KLINKE, KATRIB, C(KNMLB),
     &    IA(KLPTR), ISHEX, IA(KHEXID), *170)

C     Reset number of element blocks
        NELBLK = NELBLK + ISHEX
C     Reset number of elements
        INEL = 0
        do 20 I = 1, nelblk
          INEL = INEL + IA(KNELB+I-1)
 20     CONTINUE
        NUMEL = INEL

      ENDIF

C   --Scan element number map (global id)
      CALL MDRSRV ('MAPEL', KMAPEL, NUMEL)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160
      if (mapel) then
        call exgenm (ndb, ia(kmapel), ierr)
      else
        call iniseq(numel, ia(kmapel))
      end if

C   --Read node number map (global id)
      CALL MDRSRV ('MAPND', KMAPND, NUMNP)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160
      if (mapnd) then
        call exgnnm (ndb, ia(kmapnd), ierr)
      else
        call iniseq(numnp, ia(kmapnd))
      end if

C   --Change number of elements per element block to block index

C   --SCALER and MSMEMY and MSGEOM use MDFIND to find LENE
      CALL MDRSRV ('LENE', KLENE, 1+NELBLK)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      CALL NUM2IX (NELBLK, A(KNELB), A(KLENE))

C   --Special fix for 8-node 2D elements - order nodes consecutively
      IF (.NOT. IS3DIM) THEN
        CALL ORD8NP (NELBLK, A(KLENE), A(KNLNKE), A(KLINKE))
      END IF

C   --Make up the element to element block index

C   --MSGEOM and MSSTEP use MDFIND to find IE2ELB
      CALL MDRSRV ('IE2ELB', KE2ELB, NUMEL)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      CALL MAKE2B (NELBLK, A(KLENE), A(KE2ELB))

C   --Read node sets and side sets

      CALL MDRSRV ('IDNPS',  KIDNS, NUMNPS)
      CALL MDRSRV ('NNNPS', KNNNS, NUMNPS)
      CALL MDRSRV ('NDNPS', KNDNPS, NUMNPS)
C   --DBLIST and SSMAIN use MDFIND to find IXNNPS, LTNNPS, and FACNPS
      CALL MDRSRV ('IXNNPS', KIXNNS, NUMNPS)
C   --In exgcns call, must have space even if no dist factors defined
      CALL MDRSRV ('IXDNPS', KIXDNS, NUMNPS)
      CALL MCRSRV ('NSNAME', KNMNS,  NUMNPS * NAMLEN)
      CALL MDRSRV ('LTNNPS', KLTNNS, LNPSNL)
      CALL MDRSRV ('FACNPS', KFACNS, LNPSDF)
      CALL MDRSRV ('IDESS', KIDSS, NUMESS)
      CALL MDRSRV ('NEESS', KNESS, NUMESS)
      CALL MDRSRV ('NNESS', KNNSS, NUMESS)
      CALL MDRSRV ('NDESS', KNDSS, NUMESS)
C   --DBLIST and SSMAIN use MDFIND to find IXEESS, IXNESS, LTEESS, LTNESS,
C   --and FACESS
      CALL MDRSRV ('IXEESS', KIXESS, NUMESS)
      CALL MDRSRV ('IXNESS', KIXNSS, NUMESS)
      CALL MDRSRV ('IXDESS', KIXDSS, NUMESS)
      CALL MDRSRV ('LTNNSS', KLTNNN, LESSEL)
      CALL MDRSRV ('LTEESS', KLTESS, LESSEL)
      CALL MDRSRV ('LTNESS', KLTNSS, LESSNL)
      CALL MDRSRV ('LTSESS', KLTSSS, LESSEL)
      CALL EXINQ (NDB, EXSSDF, LESSDL, rdum, cdum, ierr)
      CALL MDRSRV ('FACESS', KFACSS, LESSDL)
      CALL MCRSRV ('SSNAME', KNMSS,  NUMESS * NAMLEN)
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      if (numnps .gt. 0) then
        call exgcns(ndb, a(kidns), a(knnns), a(kndnps), a(kixnns),
     *    a(kixdns), a(kltnns), a(kfacns), ierr)
        call getnams(ndb, EXNSET, numnps, c(knmns), namlen)
      end if
      if (numess .gt. 0) then
        call exgcssf(ndb, a(kidss), a(kness), a(kndss), a(kixess),
     *    a(kixdss), A(KLTeSS), a(kltsss), ierr)
        call getnams(ndb, EXSSET, numess, c(knmss), namlen)
        if (ierr .ne. 0) go to 170
      end if

C   --Read QA records

      call exinq(ndb, EXQA,   nqarec, rdum, cdum, ierr)
      call exinq(ndb, EXINFO, ninfo,  rdum, cdum, ierr)
      call mcrsrv('QAREC', kqarec, nqarec * 4 * MXSTLN)
      call mcrsrv('INFREC', kinfo, ninfo * MXLNLN)
      CALL MCSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160
      exodus = .true.
      if (nqarec .gt. 0) then
C ... Wrapper to get strings the right length
        call exgqaw(ndb, c(kqarec), ierr)
      end if
      if (ninfo .gt. 0) then
C ... Wrapper to get info record the right length
        call exginw(ndb, c(kinfo), ierr)
      end if

      CALL INISTR (4, ' ', CREATE)
      CALL INISTR (4, ' ', MODIFY)
C ... NOTE: cpyst8 must be called since cpystr knows that c() is only
C           a single character array.
      IF (NQAREC .GT. 0) CALL CPYST8 (4, C(KQAREC), CREATE)
      IF (NQAREC .GT. 1) THEN
        CALL CPYST8 (4, C(KQAREC+4*MXSTLN*(nqarec-1)), MODIFY)
      END IF

C   --Initialize (for GENESIS) and read database names

      NAMECO(1) = 'X'
      NAMECO(2) = 'Y'
      IF (IS3DIM) NAMECO(3) = 'Z'

      call exinq(ndb, EXTIMS, NTIMST, RDUM, CDUM, IERR)
      exodus = (ntimst .gt. 0)

      IF (EXODUS) THEN
C      --MSGEOM uses MDFIND to find ISEVOK (reserved in DBINAM)
        CALL DBINAM (NDB, 'CBVT', NDIM, NELBLK, NNDIM, NNELB,
     &    NVARHI, NVARGL, NVARNP, NVAREL, NVARNS, NVARSS,
     &    NAMECO, KNAMHV,
     &    KNAMGV, KNAMNV, KNAMEV, KNAMNS, KNAMSS,
     &    A, IA, KIEVOK, C, KNAMES,
     &    EXODUS, IA(KIDELB), ISHEX, KHEXID, NAMLEN, *110)

        GOTO 120

C      --Handle error reading variable names
 110    CONTINUE
        NVARHI = 0
        NVARGL = 0
        NVARNP = 0
        NVAREL = 0
        NVARNS = 0
        NVARSS = 0
        EXODUS = .FALSE.
 120    CONTINUE

      ELSE
        NNDIM = -999
        NNELB = -999
      END IF

      NVARHI = MAX (0, NVARHI)
      NVARGL = MAX (0, NVARGL)
      NVARNP = MAX (0, NVARNP)
      NVAREL = MAX (0, NVAREL)
      NVARNS = MAX (0, NVARNS)
      NVARSS = MAX (0, NVARSS)

      IF (EXODUS) THEN
        CALL PRNAME ('*', -1, namlen,
     *    NVARHI, NVARGL, NVARNP, NVAREL, NVARNS, NVARSS,
     &    C(KNAMES+NAMLEN*(KNAMHV-1)), C(KNAMES+NAMLEN*(KNAMGV-1)),
     &    C(KNAMES+NAMLEN*(KNAMNV-1)), C(KNAMES+NAMLEN*(KNAMEV-1)),
     &    C(KNAMES+NAMLEN*(KNAMNS-1)), C(KNAMES+NAMLEN*(KNAMSS-1)))
      END IF

C   --Get database time (allocates dynamic memory)
C   --SCALER uses MDFIND to find TIMES and WHOTIM

      if (EXODUS) then
        CALL DBITIM (NDB, '*', EXODUS,
     &    NVARNP, NELBLK, NVAREL, A(KIEVOK),
     &    NSTEPS, NSTEPW, A, KTIMES, KWHOLE, *170)
        CALL PRTIMS ('NM', -1, .TRUE., .TRUE.,
     &    NSTEPS, A(KTIMES), A(KWHOLE))
      else
C ... The 'TIMES' and 'WHOTIM' arrays are accessed even if there are
C     no timesteps on the model (See plcomd.f).
        call mdrsrv('TIMES', KTIMES, 1)
        call mdrsrv('WHOTIM', KWHOLE, 1)
      end if

      IF (.NOT. EXODUS) THEN
        WRITE (*, 10010) 'Database is in the GENESIS format'
      END IF

C   --Calculate the element centers

C      --SCALER uses MDFIND to find XE, YE, ZE
      CALL MDGET (NDIM * NUMEL)
      CALL MDRSRV ('XE', KXE, NUMEL)
      CALL MDRSRV ('YE', KYE, NUMEL)
      IF (IS3DIM) THEN
        CALL MDRSRV ('ZE', KZE, NUMEL)
      ELSE
        KZE = 1
      END IF
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      CALL ELECOR (NDIM, NELBLK, A(KLENE), A(KNLNKE), A(KLINKE),
     &  A(KXN), A(KYN), A(KZN), A(KXE), A(KYE), A(KZE))

C   --Break 3D elements into faces, and sort faces by element block

      IF (MESHOK) THEN
        IF (IS3DIM) THEN
          WRITE (*, 10010)
     &      'Please wait while the geometry is processed'
        END IF

        IF (.NOT. IS3DIM) THEN
          CALL MDRSRV ('LENF', KLENF, 1 + NELBLK + 2)
          CALL INIINT (1+NELBLK+2, 0, IA(KLENF))
        ELSE
          CALL MDRSRV ('LENF', KLENF, 1 + NELBLK + 4)
          CALL INIINT (1+NELBLK+4, 0, IA(KLENF))
        END IF
        CALL MDRSRV ('NLNKF', KNLNKF, NELBLK)
        CALL MDRSRV ('LINKF', KLINKF, 0)
        CALL MDRSRV ('IF2EL', KIF2EL, 0)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

        CALL MSSURF (A, A(KLENE), A(KNLNKE), A(KLINKE),
     &    A(KLENF), A(KNLNKF), KLINKF, KIF2EL, C(KNMLB))
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KNLNKF = 1
        KLENF = 1
        KLINKF = 1
        KIF2EL = 1
      END IF

C   --Group lines connecting nodes for efficiency

      IF (MESHOK) THEN
        CALL MDRSRV ('LENL', KLENL, 3 + NELBLK)
        CALL INIINT (3+NELBLK, 0, IA(KLENL))
        CALL MDRSRV ('LINSET', KLNSET, 0)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

        CALL MSLINS (A,
     &    A(KLENF), A(KNLNKF), A(KLINKF), A(KIF2EL),
     &    A(KLENL), KLNSET)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KLENL = 1
        KLNSET = 1
      END IF

C   --Calculate (and initialize) the deformed mesh limits

      RDMESH(1) = 0.0
      RDMESH(2) = 0.0
      RDMESH(3) = 0.0
      RDMESH(4) = 0.0

      IF (MESHOK) THEN

C      --Compute NPSURF nodes that determine the mesh limits

        NUMNPF = NUMNP

C      --QNPICK uses NPFIND to find NPSURF
        CALL MDRSRV ('NPSURF', KNPSUR, NUMNPF)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

        IF (.NOT. IS3DIM) THEN
          CALL MAKSU2 (A(KLENL), A(KLNSET), MSHBOR,
     &      .FALSE., IDUM, A(KNPSUR))
        ELSE
          CALL MAKSUR (A(KLENF), A(KNLNKF), A(KLINKF),
     &      .FALSE., IDUM, A(KNPSUR))
        END IF

        CALL MDLONG ('NPSURF', KNPSUR, NNPSUR)

C         --Find the displacement variables

        CALL FNDDIS (NAMECO, C(KNAMES+NAMLEN*(KNAMNV-1)),
     &    DEFOK, IXDEF, IYDEF, IZDEF, DEFFAC, NAMLEN)

C         --Find the undeformed mesh limits

        CALL MINMXS (NNPSUR, A(KNPSUR), A(KXN),
     &    UNMESH(KLFT), UNMESH(KRGT))
        CALL MINMXS (NNPSUR, A(KNPSUR), A(KYN),
     &    UNMESH(KBOT), UNMESH(KTOP))
        IF (IS3DIM) CALL MINMXS (NNPSUR, A(KNPSUR), A(KZN),
     &    UNMESH(KNEA), UNMESH(KFAR))

        CALL CPYREA (2*NDIM, UNMESH, ALMESH)

C         --Calculate the deformed mesh limits

        CALL DEFLIM (A, A(KWHOLE),
     &    A(KXN), A(KYN), A(KZN), A(KNPSUR))
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KNPSUR = 1
      END IF

C   --Reserve memory for all programs

      CALL MDGET (1+NUMNP + 1+NELBLK + 1+NUMEL
     &  + 1+NUMNPS + 1+NUMESS
     &  + 1+NVARHI + 1+NVARGL + 1+NVARNP + 1+NVAREL
     &  + 2+NVARHI+NVARGL)
      CALL MDRSRV ('LISNP', KLINP, 1+NUMNP)
      CALL MDRSRV ('NLISEL', KNLIEL, 1+NELBLK)
      CALL MDRSRV ('LISEL', KLIEL, 1+NUMEL)
      CALL MDRSRV ('LISNPS', KLINPS, 1+NUMNPS)
      CALL MDRSRV ('LISESS', KLIESS, 1+NUMESS)
      IF (EXODUS) THEN
        CALL MDRSRV ('LISHV', KLIHV, 1+NVARHI)
        CALL MDRSRV ('LISGV', KLIGV, 1+NVARGL)
        CALL MDRSRV ('LISNV', KLINV, 1+NVARNP)
        CALL MDRSRV ('LISEV', KLIEV, 1+NVAREL)
        CALL MDRSRV ('LIDSP', KLIDP, 2+NVARHI+NVARGL)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KLIHV = 1
        KLIGV = 1
        KLINV = 1
        KLIEV = 1
        KLIDP = 1
      END IF

      NPTIMS = 0
      NPTIMW = 0
      CALL MDRSRV ('IPTIMS', KPTIMS, MAX (NSTEPS, 1))
      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

C   --Reserve memory for mesh plots

      CALL MDRSRV ('BLKCOL', KBKCOL, 1+NELBLK)
      IF (MESHOK) THEN
        CALL MDRSRV ('IELBST', KELBST, NELBLK)
        CALL MDRSRV ('ISSNPS', KSSNPS, NUMNPS*4)
        CALL MDRSRV ('ISSESS', KSSESS, NUMESS*4)
        CALL MDRSRV ('SHDCOL', KSHDCL, NELBLK*7)
        CALL MDRSRV ('ISHDCL', KISHCL, NELBLK*3)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
        CALL INIREA (NELBLK*7, 0.0,  A(KSHDCL))
        CALL INIINT (NELBLK*3, 0,   IA(KISHCL))
      ELSE
        KELBST = 1
        KSSNPS = 1
        KSSESS = 1
        KBKCOL = 1
        KSHDCL = 1
        KISHCL = 1
      END IF

C   --Reserve memory for DETOUR, if able to run

      CALL DTCHK (.FALSE., DTOK)

      IF (DTOK) THEN
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      END IF

C   --Reserve memory for PATHLINE, if able to run

      CALL LNCHK (.FALSE., LNOK)

      IF (LNOK) THEN
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      END IF

C   --Reserve memory for TPLOT, if able to run

      CALL TPCHK (.FALSE., TPOK)

      IF (TPOK) THEN
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      END IF

C   --Reserve memory for SPLOT, if able to run

      CALL SPCHK (.FALSE., SPOK)

      IF (SPOK) THEN
        MAXNE = MAX (NUMNP, NUMEL)
        CALL MDRSRV ('NENUM', KNENUM, MAXNE)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KNENUM = 1
      END IF

      IF (LNOK .OR. SPOK) THEN
        CALL MDRSRV ('IPATH', KIPATH, 0)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160
      ELSE
        KIPATH = 1
      END IF

C        Initialize array containing list of display variables

      IF (EXODUS) CALL DISPV (.TRUE., ' ', IDUM, IDUM,
     &  ' ', C(KNAMES), A(KLIDP), NAMLEN)

C        Initialize BLKCOL array.

      CALL BCOLOR (.TRUE., ' ', IDUM, IDUM, IDUM,' ', A(KBKCOL))

C        Initialize line thicknesses for mesh plots
      CALL LINTHK (CDUM, IDUM, IDUM, IDUM, RDUM, CDUM, .TRUE.)

        write (*,9999)
 9999   FORMAT(/,
     *    10x,'NOTE: This version of blot uses global ids for both',
     *    ' node and element ids by default.',/,
     *    10x,'      To see the mapping from local to global, use',
     *    ' the commands:',/,
     *    10x,'          "LIST MAP" (element map), or ',
     *    '"LIST NODEMAP" (node map)',/,
     *    10x,'      To disable the maps and use local ids, restart',
     *    ' blot with "-nomap node|element|all"',//,
     *    10x,'      Notify gdsjaar@sandia.gov if bugs found')

        if (mapel .and. mapnd) then
          WRITE (*, 10010) 'Nodes and Elements using Global Ids'
        else if (mapel) then
          WRITE (*, 10010) 'Elements use Global Ids, Node Ids are Local'
        else if (mapnd) then
          WRITE (*, 10010) 'Element use Local Ids, Node Ids are Global'
        else
          WRITE (*, 10010) 'Nodes and Elements using Local Ids'
        end if

 130  CONTINUE
      IF (.TRUE.) THEN
        CALL MDLONG ('IPTIMS', KPTIMS, MAX (NSTEPS, 1))
        IF (SPOK) THEN
          CALL MDLONG ('NENUM', KNENUM, MAXNE)
        END IF
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

        CALL COMAND (A, CURPRO, C(KQAREC), C(KINFO),
     &    NAMECO, C(KNMLB), C(KNAMES), A(KTIMES), A(KWHOLE),
     *    A(KPTIMS), A(KMAPEL), A(KMAPND),
     &    A(KIDELB), NEWELB, A(KELBST), A(KE2ELB),
     &    A(KLENE), A(KNLNKE), A(KLINKE),
     &    A(KXN), A(KYN), A(KZN), A(KXE), A(KYE), A(KZE),
     &    A(KIEVOK),
     &    A(KSSNPS), A(KIDNS), A(KNNNS),
     &    A(KSSESS), A(KIDSS), A(KNESS), A(KNNSS),
     &    NCSTEP, A(KLINP), A(KNLIEL), A(KLIEL),
     &    A(KLINPS), A(KLIESS),
     &    A(KLIHV), A(KLIGV), A(KLINV), A(KLIEV), A(KLIDP),
     &    A(KBKCOL), A(KNENUM), NEUTRL, NEWPRO, A(KSHDCL),
     *    A(KISHCL), C(KNMEB), C(KNMNS), C(KNMSS), NAMLEN, *150)
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

C         --Calculate the undeformed and deformed mesh limits

        IF ((CURPRO .EQ. 'DETOUR') .AND. (DEFFAC .LT. 0.0)) THEN
          CALL DEFLIM (A, A(KWHOLE),
     &      A(KXN), A(KYN), A(KZN), A(KNPSUR))
          CALL MDSTAT (NERR, MEM)
          IF (NERR .GT. 0) GOTO 160
        END IF

        CALL MDLONG ('IPTIMS', KPTIMS, NPTIMS)
        IF (SPOK) THEN
          CALL MDLONG ('NENUM', KNENUM, NNENUM)
        END IF
        CALL MDSTAT (NERR, MEM)
        IF (NERR .GT. 0) GOTO 160

        IF (NEWPRO .EQ. 'P') THEN
          CALL PLOTER (A, CURPRO, NEUTRL,
     &      C(KNAMES), NPTIMS, A(KPTIMS), A(KTIMES), A(KWHOLE),
     &      A(KNLNKE), A(KIEVOK), A(KNENUM),
     &      A(KXN), A(KYN), A(KZN), A(KXE), A(KYE), A(KZE),
     &      NAMECO, A(KLENF), A(KNLNKF), KLINKF,
     &      A(KLENL), KLNSET,
     &      A(KE2ELB), NEWELB, A(KELBST), KNPSUR,
     &      A(KSSNPS), A(KIDNS), A(KSSESS), A(KIDSS),
     &      A(KLIDP), A(KBKCOL), A(KIDELB),C(KNMLB), NAMLEN,
     *      A(KMAPEL), A(KMAPND))
          CALL MDSTAT (NERR, MEM)
          IF (NERR .GT. 0) GOTO 160
        END IF

        GOTO 130
      END IF

 150  CONTINUE

C   --Finish graphics
      CALL GREXIT

      CALL MDSTAT (NERR, MEM)
      IF (NERR .GT. 0) GOTO 160

      GOTO 180

 160  CONTINUE
      WRITE (*, *) 'Total memory used = ', MEM, ' words.'
      CALL MEMERR
      GOTO 180

 170  CONTINUE
      GOTO 180

 180  CONTINUE

C   --Close files

      CALL EXCLOS(NDB, IERR)
      IF (NEUOPN) THEN
        WRITE (NEU, 10030) 'xmax', wxmax
        WRITE (NEU, 10030) 'xmin', wxmin
        WRITE (NEU, 10030) 'ymax', wymax
        WRITE (NEU, 10030) 'ymin', wymin
        CLOSE (NEU, IOSTAT=IDUM)
      END IF
      IF (ANYPRT) CLOSE (NPRT, IOSTAT=IDUM)
      IF (GRFOPN) CLOSE(NEUGRF)
      IF (CSVOPN) CLOSE(NCSV)
      call addlog (QAINFO(1)(:lenstr(QAINFO(1))))
      CALL WRAPUP (QAINFO(1))

10010 FORMAT (/, 1X, 5A)
10020 FORMAT (/
     &  21X,'BBBBBBBB    LL            OOOOOO    TTTTTTTTTT', /
     &  20X,'BBBBBBBBB   LL           OOOOOOOO   TTTTTTTTTT', /
     &  19X,'BB      BB  LL          OO      OO      TT    ', /
     &  18X,'BB      BB  LL          OO      OO      TT    ', /
     &  17X,'BBBBBBBBB   LL          OO      OO      TT    ', /
     &  16X,'BBBBBBBBB   LL          OO      OO      TT    ', /
     &  15X,'BB      BB  LL          OO      OO      TT    ', /
     &  14X,'BB      BB  LL          OO      OO      TT    ', /
     &  13X,'BBBBBBBBB   LLLLLLLLLL   OOOOOOOO       TT    ', /
     &  12X,'BBBBBBBB    LLLLLLLLLL    OOOOOO        TT    II-2')
10030 FORMAT ('@ world ',A4,1x,1pe15.7E3)
      END

      subroutine inimap(num, iar)
      dimension iar(*)
      do 10 i=1, num
        iar(i) = i
 10   continue
      return
      end

C=======================================================================
      SUBROUTINE CPYST8 (LEN, IFROM, ITO)
C=======================================================================
C   --Parameters:
C   --   LEN   - IN  - the number of strings in the list
C   --   IFROM - IN  - the input list
C   --   ITO   - OUT - the copied list
C   --   ISTLN - IN  - string length

      INTEGER LEN
      include 'params.blk'
      CHARACTER*(MXSTLN) IFROM(*), ITO(*)

      DO 100 I = 1, LEN
         ITO(I) = IFROM(I)
  100 CONTINUE

      RETURN
      END

C ... Wrapper to get strings the right length
      subroutine exgqaw(ndb, qarec, ierr)
      include 'params.blk'
      character*(mxstln) qarec(4, *)
      call exgqa(ndb, qarec, ierr)
      return
      end
      subroutine exginw(ndb, info, ierr)
      include 'params.blk'
      character*(mxlnln) info(*)
      call exginf(ndb, info, ierr)
      return
      end

C=======================================================================
      SUBROUTINE RDTYPE (NELBLK,NAMELB,IDELB,NUMELB,ISHEX,HEXID,NSHL)
C=======================================================================

      INTEGER NELBLK, ISHEX
      INTEGER IDELB(*), HEXID(*), numelb(*)
      CHARACTER*32 NAMELB(*)

C     Initialize HEX return variables
      ISHEX = 0
      NSHL = 0
C     Mark HEXID array with corresponding element blocks that have HEXSHELLs
      DO 310 I = 1, NELBLK
C        Check if element block contains HEXSHELLs
         IF (NAMELB(I)(1:8) .EQ. 'HEXSHELL') THEN
C           ISHEX - counter: how many HEXSHELL element blocks
C           index - HEXID(ISHEX) - stores HEXSHELL element block id
            ISHEX = ISHEX + 1
            HEXID(ISHEX) = IDELB(I)
         ELSE IF (NAMELB(I)(1:5) .EQ. 'SHELL') THEN
            NSHL = NSHL + NUMELB(I)
         ENDIF
 310  CONTINUE

      RETURN
      END

C=======================================================================
      SUBROUTINE GETNAMS (NDB, ITYPE, NUM, NAMES, NAMLEN)
C=======================================================================

      CHARACTER*(NAMLEN) NAMES(*)
      CALL EXGNAMS(NDB, ITYPE, NUM, names, ierr)
      RETURN
      END

      subroutine iniseq(icnt, map)
      integer map(*)
      do i=1, icnt
        map(i) = i
      end do
      return
      end

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


text:
module m1

end module m1

max_stars_repo_name: clementval/claw-compiler

max_issues_repo_name: skosukhin/claw-compiler

max_forks_repo_name: skosukhin/claw-compiler


text:
PROGRAM Assignment1
IMPLICIT NONE

INTEGER:: i, j
INTEGER, PARAMETER :: m = 1050
REAL :: Prior(0:m), P(0:m), TauArray(0:m), DecayMeasurments(10), l(1), dtau, NormalisationConstant, Y(1), X(1)

dtau = 0.1; Prior=1.0/99.0; P = Prior

OPEN(10, FILE = 'bayes_data.txt') !Opening decay measurment times
  READ(10,'(f6.3)') DecayMeasurments !Reading 10 decay measurments to array
CLOSE(10)

DO i = 1, 1050
  TauArray(i) = REAL(i)*dtau !Creating tauarray, the possible half life times, in intervals of 0.1s
END DO

DO j = 1, 10
  DO i = 1, 1050

    P(i)=P(i)*exp(-DecayMeasurments(j)/TauArray(i))/TauArray(i)

  END DO
END DO

NormalisationConstant = (SUM(P) - 0.5*(P(0)+P(m)))*dtau !Calculating normalisation constant

DO i = 1, 1050
  P(i) = P(i)/NormalisationConstant!Normalising probability array
END DO

OPEN(11, FILE = 'Results.dat') !Creating file to print results in
DO i = 1,1050
  WRITE(11,*) 'p ''Results.dat'' u 1:2 w l'
  WRITE(11,*) TauArray(i), P(i) !Printing half life times are their respective probabaility to file
END DO
CLOSE(11)

WRITE(6,*) 'Results printed to file ''Results.dat''', SUM(P)

l = Maxloc(P); Y = P(l); X = TauArray(l)
WRITE(6,*) L

WRITE(6,*) 'Most probably half life is', X, 'with a probability of:', Y

END PROGRAM

max_stars_repo_path: Monte Carlo Method and Bayesian Statistics/Assignment1.f90

max_stars_repo_name: DGrifferty/Fortran

max_issues_repo_name: DGrifferty/Fortran

max_forks_repo_name: DGrifferty/Fortran


text:
c----------------------------------------------------------------
c----------------------------------------------------------------
c--- The car axis problem of radau
c----------------------------------------------------------------
c----------------------------------------------------------------

c -------- radaudae.f -> radaudae.dll ------
c compile in R with: system("g77 -shared -o radaudae.dll radaudae.f")
c or with system("R CMD SHLIB radaudae.f")

c----------------------------------------------------------------
c Initialiser for parameter common block
c----------------------------------------------------------------
      subroutine initcaraxis(daeparms)

      external daeparms
      integer, parameter :: N = 8
      double precision parms(N)
      common /myparms/parms

      call daeparms(N, parms)
      return
      end

c----------------------------------------------------------------
c rate of change
c----------------------------------------------------------------
      subroutine caraxis(neq, t, y, ydot, out, ip)
      implicit none
      integer          neq, IP(*)
      double precision t, y(neq), ydot(neq), out(*)
      double precision eps, M, k, L, L0, r, w, g
      common /myparms/ eps, M, k, L, L0, r, w, g

      double precision xl, yl, xr, yr, ul, vl, ur, vr, lam1, lam2
      double precision yb, xb, Ll, Lr, dxl, dyl, dxr, dyr
      double precision dul, dvl, dur, dvr, c1, c2

c expand state variables
       xl = y(1)
       yl = y(2)
       xr = y(3)
       yr = y(4)
       ul = y(5)
       vl = y(6)
       ur = y(7)
       vr = y(8)
       lam1 = y(9)
       lam2 = y(10)

       yb = r * sin(w * t)
       xb = sqrt(L * L - yb * yb)
       Ll = sqrt(xl**2 + yl**2)
       Lr = sqrt((xr - xb)**2 + (yr - yb)**2)

       dxl = ul
       dyl = vl
       dxr = ur
       dyr = vr

       dul = (L0-Ll) * xl/Ll      + 2 * lam2 * (xl-xr) + lam1*xb
       dvl = (L0-Ll) * yl/Ll      + 2 * lam2 * (yl-yr) + lam1*yb - k*g
       dur = (L0-Lr) * (xr-xb)/Lr - 2 * lam2 * (xl-xr)
       dvr = (L0-Lr) * (yr-yb)/Lr - 2 * lam2 * (yl-yr) - k*g

       c1  = xb * xl + yb * yl
       c2  = (xl - xr)**2 + (yl - yr)**2 - L * L

c function values in ydot
       ydot(1)  = dxl
       ydot(2)  = dyl
       ydot(3)  = dxr
       ydot(4)  = dyr
       ydot(5)  = dul
       ydot(6)  = dvl
       ydot(7)  = dur
       ydot(8)  = dvr
       ydot(9)  = c1
       ydot(10) = c2
      return
      end

ext: f


text:
  subroutine dasray_intf( &
       n,                 nsmax,        ns_elas,          &
       maxintstr,         maxstr,       maxstr_internal,  &
       ngeomax,           mmax,         maxleg,           &
       maxlay,            nstokes,      accuracy,         &
       dnu,               tau,          omega_Ray,        &
       omega_Cab,         om_air,       scat_type,        &
       weightCabRam_air,  PCab_N2,      PRay_N2,          &
       PCab_O2,           PRay_O2,      PCab_Ar,          &
       PRay_Ar,           PCab_air,     PRay_air,         &
       PRam_air,          maxcoeffs,    maxcoeffstrunc,   &
       surf_type,         alb_Lambert,  gmu,              &
       gwt,               gmumin,       phi,              &
       imuv,              imu0,         plm_0,            &
       plm_m,             plm_p,        u,                &
       refl)

!--------------------------------------------------------------
! Subroutine:
!   Doubling-adding solver (elastic Rayleigh scattering)
!
! Purpose:
!   Solves radiative transfer in plane-parallel atmosphere
!   in the scalar approximation and includes inelastic
!   Raman scattering
!
! Author:
!   Rutger van Deelen
!   SRON - Netherlands Insititute for Space Research
!
! Reference:
!   R. van Deelen, J. Landgraf, and I. Aben --
!   "Multiple elastic and inelastic light-scattering in the
!   Earth's atmosphere: A doubling-adding method to include
!   rotational Raman scattering by air",
!   J. Quant. Spectrosc. Radiat. Transfer 95, 399-433, 2005
!------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

! Input ==========================================================

  integer,intent(in):: &
     n,          & !wavenr index, coarse grid
     nsmax,      & !nr of wavenrs 0:nsmax
     ns_elas,    & !wavenr index corresponding to elastically scattererd light
     maxintstr,  & !Nr of internal streams (must be even and >= 2)
     maxstr,     &
     maxstr_internal, &
     ngeomax,    & !Nr of viewing geometries
     mmax,       & !maximum Fourier moment (molecular scattering=2)
     maxleg,     & !maximum Legendre moment (molecular scatt.=2)
     maxlay,     & !Nr of atmospheric layers
     nstokes       !NR of stokes parameters

! wavenumber binwdith [cm-1]
  real(kind=dp),intent(in):: &
     dnu

! viewing geometries
  real(kind=dp)::  &
     phi(ngeomax)

! Gaussian quadrature
  real(kind=dp)::gmu(maxstr),gwt(maxstr),gmumin

! points to stream index for viewing angles and solar angles
  integer,intent(in):: &
     imuv(ngeomax), &
     imu0(ngeomax)

! needed for polarization and single scattering solution
  real(kind=dp),dimension(0:maxleg,0:maxleg,maxintstr+4*ngeomax),intent(in):: &
     plm_0, &
     plm_m, &
     plm_p
  real(kind=dp),dimension(maxintstr+4*ngeomax),intent(in):: &
     u

! doubling-adding accuracy
  real(kind=dp),intent(in):: &
      accuracy

! Internal ========================================================

  logical,parameter:: &
      single_scattering = .false.

  integer:: &
    m,   & !index of Fourier moment m
    k,   & !layer index
    ik,  & !layer index, reverse order
    ikk, & !layer index (fine layers for doubling)
    i,i1,i2,ns,ns1,ns2,lg,i_st,iv,is

!--------------------------------------------------------------
! Optical depth, single scattering albedo, etc.

  real(kind=dp),dimension(maxlay,0:nsmax),intent(in):: &
     tau,       & !total optical depth
     omega_Ray, & !Rayleigh single scattering albedo
     omega_Cab    !Cabannes single scattering albedo

  real(kind=dp),intent(in):: &
     om_air(maxlay,0:nsmax,0:nsmax)

! 0=not good; 1=Cab; 2=Ram; 3=mixCab&Ram; 4=Ray
  integer,intent(in):: &
     scat_type(0:nsmax,0:nsmax)

! weight CabRam for eff. phasematrix
  real(kind=dp),intent(in):: &
     weightCabRam_air(maxlay,0:nsmax)

!--------------------------------------------------------------
! phase matrices, depend on i,j, m, k and n
  real(kind=dp),intent(in), &
   dimension(nstokes,nstokes,0:maxleg,maxlay):: &
    PCab_N2,  & !Cabannes phase fct. moments
    PRay_N2,  & !Rayleigh phase fct. moments
    PCab_O2,  & !Cabannes phase fct. moments
    PRay_O2,  & !Rayleigh phase fct. moments
    PCab_Ar,  & !Cabannes phase fct. moments
    PRay_Ar,  & !Rayleigh phase fct. moments
    PCab_air, & !Cabannes phase fct. moments
    PRay_air, & !Rayleigh phase fct. moments
    PRam_air    !Raman    phase fct. moments

! Nr of Legendre terms for each layer
! Nr of coefficient per layer (clear sky = 0:2)
  integer,intent(in):: &
    maxcoeffs(maxlay), &    !all terms
    maxcoeffstrunc(maxlay)  !truncate sum for diffuse radiation

  real(kind=dp), &
   dimension(0:maxleg):: &
    wl_Cab_air, & !Cabannes phase fct. moments
    wl_Ray_air, & !Rayleigh phase fct. moments
    wl_Ram_air    !Raman    phase fct. moments

! phase matrices, depend on n, m
  real(kind=dp):: &
     Zm_Ray_min(maxstr,maxstr),   &
     Zm_Cab_min(maxstr,maxstr),   &
     Zm_Ram_min(maxstr,maxstr),   &
     Zm_Ray_plus(maxstr,maxstr),  &
     Zm_Cab_plus(maxstr,maxstr),  &
     Zm_Ram_plus(maxstr,maxstr)

! mixture of Cab and Ram, depends on n, m, k, mix for each ns
  real(kind=dp):: &
     Zm_CabRam_min(0:nsmax,maxstr,maxstr),  &
     Zm_CabRam_plus(0:nsmax,maxstr,maxstr)

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

  integer:: &
     ndl !nr of doublings

  real(kind=dp):: &
     dtau(0:nsmax)    !total (sca+abs) optical depth per layer

  real(kind=dp):: &
     dtau_A(0:nsmax),  &
     tau_atm(0:nsmax), &
     tau_AB(0:nsmax),  &
     tau_wl(0:nsmax)

!-----------------------------------------------------------
! Reflectance

! whole layer ,  single scattering approximation
  real(kind=dp):: &
     R1m_wl(0:nsmax,0:nsmax,ngeomax),     &
     T1m_wl(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_wl(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_wl(0:nsmax,0:nsmax,ngeomax)

! atm only
  real(kind=dp):: &
     R1m_atm(0:nsmax,0:nsmax,ngeomax),     &
     T1m_atm(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_atm(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_atm(0:nsmax,0:nsmax,ngeomax)

! atm+surf
  real(kind=dp):: &
     R1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     T1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     D1m_sys(0:nsmax,0:nsmax,ngeomax),     &
     R1mStar_sys(0:nsmax,0:nsmax,ngeomax), &
     T1mStar_sys(0:nsmax,0:nsmax,ngeomax)

! initial layer to start doubling
  real(kind=dp):: &
     R1m_bl(0:nsmax,0:nsmax,maxstr,maxstr), &
     T1m_bl(0:nsmax,0:nsmax,maxstr,maxstr)

  real(kind=dp):: &
     Rm_A(0:nsmax,0:nsmax,maxstr,maxstr),       &
     Tm_A(0:nsmax,0:nsmax,maxstr,maxstr),       &
     Rm_AB(0:nsmax,0:nsmax,maxstr,maxstr),      &
     Tm_AB(0:nsmax,0:nsmax,maxstr,maxstr),      &
     RmStar_AB(0:nsmax,0:nsmax,maxstr,maxstr),  &
     TmStar_AB(0:nsmax,0:nsmax,maxstr,maxstr),  &
     Rm_atm(0:nsmax,0:nsmax,maxstr,maxstr),     &
     Tm_atm(0:nsmax,0:nsmax,maxstr,maxstr),     &
     RmStar_atm(0:nsmax,0:nsmax,maxstr,maxstr), &
     TmStar_atm(0:nsmax,0:nsmax,maxstr,maxstr)

   real(kind=dp):: &
     Rm_sys(0:nsmax,0:nsmax,maxstr,maxstr),     &
     Tm_sys(0:nsmax,0:nsmax,maxstr,maxstr),     &
     RmStar_sys(0:nsmax,0:nsmax,maxstr,maxstr), &
     TmStar_sys(0:nsmax,0:nsmax,maxstr,maxstr)

! refl single scattering
  real(kind=dp):: &
    reflRay_ss(0:nsmax,nstokes,ngeomax),  &
    reflCab_ss(0:nsmax,nstokes,ngeomax),  &
    reflRam_ss(0:nsmax,nstokes,ngeomax),  &
    reflmRay_ss(0:nsmax,nstokes,ngeomax), &
    reflmCab_ss(0:nsmax,nstokes,ngeomax), &
    reflmRam_ss(0:nsmax,nstokes,ngeomax), &
    Zm_Ray_ss(nstokes,ngeomax), &
    Zm_Cab_ss(nstokes,ngeomax)

! aux
  real(kind=dp):: &
    tau_above(0:nsmax), &
    expuv(0:nsmax,ngeomax), &
    expu0(0:nsmax,ngeomax), &
    taut, &
    explay, &
    factor, &
    one_4pi, &
    uv, &
    u0

 real(kind=dp):: &
    deg2rad, &
    rad2deg, &
    pi

 real(kind=dp):: &
    one_d0m, &
    dfac, & !0.5 * (2-d0m)
    Bmplus(4,ngeomax), &
    Bmmin(4,ngeomax), &
    EpDelta(4), &
    EmDelta(4)

  real(kind=dp):: &
     R_sys_ns(0:nsmax,ngeomax)

  real(kind=dp),intent(out):: &
     refl(0:nsmax,nstokes,ngeomax)

! Surface reflection ==============================================

! Choose type of surface reflection
!   0 = no surface
!   1 = Lambertian surface reflection
!   2 = BDRF

  integer,intent(in):: &
     surf_type

! Lambertian surface albedo
  real(kind=dp),intent(in):: &
     alb_Lambert

!===========================================================

  pi = 2.d0*asin(1.d0)
  one_4pi = 1.d0/(4.d0*pi)

  rad2deg = 180.d0/pi
  deg2rad = pi/180.d0

!--------------------------------------------------------------------
! Reset
!--------------------------------------------------------------------

  do ns=0,nsmax
     do iv=1,ngeomax
        R_sys_ns(ns,iv) = 0.d0
     enddo
  enddo

!----------------------------------------------------------------------
! Set these to zero.
! (Have been filled above
!  in subroutines ss_elastic and ss_inelastic)
!---------------------------------------------------------------------

  do iv=1,ngeomax
     do i_st = 1,nstokes
        do ns=0,nsmax
           refl(ns,i_st,iv)       = 0.d0
           reflRay_ss(ns,i_st,iv) = 0.d0
           reflCab_ss(ns,i_st,iv) = 0.d0
           reflRam_ss(ns,i_st,iv) = 0.d0
        enddo
     enddo
  enddo

!********************************************************************
! Loop over Fourier moments
!********************************************************************

  do m=0,mmax

!--------------------------------------------------------------------
!    Setup quantities needed for Fourier expansion
!--------------------------------------------------------------------

     call Fourier_init( &
          m,        ngeomax,  phi,     &
          one_d0m,  dfac,     EpDelta, &
          EmDelta,  Bmplus,   Bmmin)

!********************************************************************
!    Loop over layers with index ik
!    from bottom to top (reverse order)
!
!    Optical properties are stored from top to bottom
!    corresponding layer index k = maxlay+1-ik
!********************************************************************

!    reverse order

     do ik=1,maxlay !reverse order

        k = maxlay+1-ik !corresponding k

!--------------------------------------------------------------------
!       Save phase matrix coefficients
!--------------------------------------------------------------------

        do lg=0,maxleg
           wl_Ray_air(lg) = pRay_air(1,1,lg,k)
           wl_Cab_air(lg) = pCab_air(1,1,lg,k)
           wl_Ram_air(lg) = pRam_air(1,1,lg,k)
        enddo

!-------------------------------------------------------------------
!       Determinine number of doublings,
!       the number of sublayers and their
!       optical thicknesses for each layer ik,
!       and for each Fourier component m,
!       and for each incident wavenumber n
!-------------------------------------------------------------------

        call nr_doublings( &
             accuracy,   ns_elas,  nsmax,      &
             k,          maxlay,   m,          &
             maxleg,     gmumin,   wl_Ray_air, &
             omega_Ray,  tau,      ndl)

!--------------------------------------------------------------------
!       Setup phasefunction moments Zm
!--------------------------------------------------------------------

        call Zm_setup( &
             nsmax,                        &
             k,             maxlay,            m,             &
             mmax,          maxstr,            maxleg,        &
             gmu,           weightCabram_air,  wl_Ray_air,    &
             wl_Cab_air,    wl_Ram_air,        Zm_Ray_min,    &
             Zm_Ray_plus,   Zm_Cab_min,        Zm_Cab_plus,   &
             Zm_Ram_min,    Zm_Ram_plus,       Zm_CabRam_min, &
             Zm_CabRam_plus)

!--------------------------------------------------------------------
!       Setup optical depth dtau for all wavenumbers around
!       wavenumber of incoming light
!--------------------------------------------------------------------

!       set for every layer ik, every m, and every wavenr ns
        do ns=0,nsmax
           dtau(ns) = tau(k,ns)/(2.d0**dble(ndl))
           tau_atm(ns) = 0.d0
        enddo

!********************************************************************
!       Single scattering solution of layer
!********************************************************************

!--------------------------------------------------------------------
!       ss, update
!--------------------------------------------------------------------

        call ss_Zm_setup( &
             m,         k,          nstokes,         &
             ngeomax,   maxstr,     maxintstr,       &
             maxleg,    maxlay,     maxcoeffstrunc,  &
             gmu,       imuv,       imu0,            &
             plm_0,     plm_m,      plm_p,           &
             u,         pRay_air,   Zm_Ray_ss)

        call ss_Zm_setup( &
             m,         k,          nstokes,         &
             ngeomax,   maxstr,     maxintstr,       &
             maxleg,    maxlay,     maxcoeffstrunc,  &
             gmu,       imuv,       imu0,            &
             plm_0,     plm_m,      plm_p,           &
             u,         pCab_air,   Zm_Cab_ss)

        do ns=0,nsmax
           tau_above(ns) = 0.d0
           if(ik.lt.maxlay)then
              do i=ik+1,maxlay
                 tau_above(ns) = tau_above(ns) + tau(maxlay+1-i,ns)
              enddo
           endif
        enddo

        do iv=1,ngeomax

           u0 = gmu(imu0(iv))
           uv = gmu(imuv(iv))

           do ns=0,nsmax

              expu0(ns,iv)  = exp(-tau_above(ns)/u0)
              expuv(ns,iv)  = exp(-tau_above(ns)/uv)

              taut   = tau(k,ns)
              explay = exp(-taut * (1.d0/uv + 1.d0/u0))
              factor = one_4pi *  u0/(u0+uv) * &
                          expu0(ns,iv)*expuv(ns,iv) * &
                           (1.d0-explay)

              do i_st=1,nstokes
                 reflmRay_ss(ns,i_st,iv) = &
                        reflmRay_ss(ns,i_st,iv) + &
                         omega_Ray(k,ns) * &
                          Zm_Ray_ss(i_st,iv) * &
                           factor
              enddo

           enddo

        enddo

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

        if(single_scattering)then

           call whole_sublayer( &
                ngeomax,        imuv,            imu0,         &
                maxstr,         maxlay,          nsmax,        &
                k,            &
                om_air,         scat_type,       tau,          &
                Zm_Ray_min,     Zm_Ray_plus,     Zm_Cab_min,   &
                Zm_Cab_plus,    Zm_Ram_min,      Zm_Ram_plus,  &
                Zm_CabRam_min,  Zm_CabRam_plus,  gmu,          &
                R1m_wl,         T1m_wl)

!--------------------------------------------------------------------
!          Adding layer A on top of layers below
!--------------------------------------------------------------------

!          Skip adding for first layer

           if(ik.eq.1)then

              do ns=0,nsmax
                 tau_atm(ns) = tau(k,ns)
              enddo

              do iv=1,ngeomax
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       R1mStar_wl(ns1,ns2,iv)  = R1m_wl(ns1,ns2,iv)
                       R1mStar_wl(ns1,ns2,iv)  = T1m_wl(ns1,ns2,iv)
                       R1m_atm(ns1,ns2,iv)     = R1m_wl(ns1,ns2,iv)
                       T1m_atm(ns1,ns2,iv)     = T1m_wl(ns1,ns2,iv)
                       R1mStar_atm(ns1,ns2,iv) = R1mStar_wl(ns1,ns2,iv)
                       T1mStar_atm(ns1,ns2,iv) = T1mStar_wl(ns1,ns2,iv)
                    enddo
                 enddo
              enddo

           else

              do ns=0,nsmax
                 tau_wl(ns)  = tau(k,ns)   !tau_A  layer ik
                 tau_atm(ns) = tau(k+1,ns) !tau_B  layer ik-1
              enddo

              do iv=1,ngeomax
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       R1mStar_wl(ns1,ns2,iv) = R1m_wl(ns1,ns2,iv)
                       T1mStar_wl(ns1,ns2,iv) = T1m_wl(ns1,ns2,iv)
                    enddo
                 enddo
              enddo

              call adder1(  &
                   ngeomax,     maxstr,       imuv,         &
                   imu0,        mmax,         nsmax,        &
                   ns_elas,     m,            gmu,          &
                   R1m_wl,      T1m_wl,       R1mStar_wl,   &
                   T1mStar_wl,  tau_wl,       R1m_atm,      &
                   T1m_atm,     R1mStar_atm,  T1mStar_atm,  &
                   tau_atm)

           endif

        else

!********************************************************************
!          Multiple scattering solution of layer
!********************************************************************

!--------------------------------------------------------------------
!          Setup initial sublayer for doubling
!--------------------------------------------------------------------

           call basic_sublayer( &
                maxstr,      maxlay,        nsmax,           &
                ns_elas,     n,               &
                k,           om_air,        scat_type,       &
                dtau,        Zm_Ray_min,    Zm_Ray_plus,     &
                Zm_Cab_min,  Zm_Cab_plus,   Zm_Ram_min,      &
                Zm_Ram_plus, Zm_CabRam_min, Zm_CabRam_plus,  &
                gmu,         R1m_bl,        T1m_bl)

           do i2=1,maxstr
              do i1=1,maxstr
                 do ns2=0,nsmax
                    do ns1=0,nsmax
                       Rm_A(ns1,ns2,i1,i2) = R1m_bl(ns1,ns2,i1,i2)
                       Tm_A(ns1,ns2,i1,i2) = T1m_bl(ns1,ns2,i1,i2)
                    enddo
                 enddo
              enddo
           enddo

           do ns=0,nsmax
              dtau_A(ns) = dtau(ns)
           enddo

!----------------------------------------------------------------
!          Start doubling
!--------------------------------------------------------------------

           do ikk = 1,ndl

                 call doublerray( &
                      maxstr,     maxstr_internal,  mmax,   &
                      nsmax,      ns_elas,          dnu,    &
                      m,          one_d0m,          ik,     &
                      ikk,        gmu,              gwt,    &
                      accuracy,   Rm_A,             Tm_A,   &
                      dtau_A,     Rm_AB,            Tm_AB)

                 do i2=1,maxstr
                    do i1=1,maxstr
                       do ns2=0,nsmax
                          do ns1=0,nsmax
                             Rm_A(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                             Tm_A(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                          enddo
                       enddo
                    enddo
                 enddo

                do ns=0,nsmax
                  dtau_A(ns) = 2.d0 * dtau_A(ns)
                enddo

           enddo

!--------------------------------------------------------------------
!          Adding layer on top of layer(s) below
!--------------------------------------------------------------------

!          Skip adding for first layer

           if(ik.eq.1)then

               do ns=0,nsmax
                  tau_atm(ns) = tau(k,ns)
               enddo

!              Total homogenous atmosphere layer
               do i2=1,maxstr
                  do i1=1,maxstr
                     do ns2=0,nsmax
                        do ns1=0,nsmax
                           Rm_atm(ns1,ns2,i1,i2)     = Rm_AB(ns1,ns2,i1,i2)
                           Tm_atm(ns1,ns2,i1,i2)     = Tm_AB(ns1,ns2,i1,i2)
                           RmStar_atm(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                           TmStar_atm(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                        enddo
                     enddo
                  enddo
               enddo

           else

              do ns=0,nsmax
                 tau_AB(ns)  = tau(k,ns)   !tau_A layer ik
                 tau_atm(ns) = tau(k+1,ns) !tau_B layer ik-1
              enddo

!             star
              do i2=1,maxstr
                 do i1=1,maxstr
                    do ns2=0,nsmax
                       do ns1=0,nsmax
                          RmStar_AB(ns1,ns2,i1,i2) = Rm_AB(ns1,ns2,i1,i2)
                          TmStar_AB(ns1,ns2,i1,i2) = Tm_AB(ns1,ns2,i1,i2)
                       enddo
                    enddo
                 enddo
              enddo

              call adderray( &
                   maxstr,     maxstr_internal, mmax,      &
                   nsmax,      ns_elas,         dnu,       &
                   m,          one_d0m,         gmu,       &
                   gwt,        accuracy,        Rm_AB,     &
                   Tm_AB,      RmStar_AB,       TmStar_AB, &
                   tau_AB,     Rm_atm,          Tm_atm,    &
                   RmStar_atm, TmStar_atm,      tau_atm)

          endif !first layer or not

        endif !multiple/single scat

     enddo !end loop over layers

!**********************************************************************
!    Adding surface
!**********************************************************************

!    total optical depth atmosphere above surface
     do ns=0,nsmax
        tau_atm(ns) = 0.d0
        do k=1,maxlay
           tau_atm(ns) = tau_atm(ns) + tau(k,ns)
        enddo
     enddo

     call add_surf(    &
          single_scattering,  surf_type,        alb_Lambert,   &
          m,                  mmax,             one_d0m,       &
          maxstr,             maxstr_internal,  nsmax,         &
          ns_elas,            dnu,              ngeomax,       &
          gmu,                gwt,              imuv,          &
          imu0,               accuracy,         tau_atm,       &
          R1m_atm,            T1m_atm,          R1mStar_atm,   &
          T1mStar_atm,        R1m_sys,          T1m_sys,       &
          D1m_sys,            R1mStar_sys,      T1mStar_sys,   &
          Rm_atm,             Tm_atm,           RmStar_atm,    &
          TmStar_atm,         Rm_sys,           Tm_sys,        &
          RmStar_sys,         TmStar_sys)

!    ss

     if(m.eq.0)then

        do iv=1,ngeomax

           u0 = gmu(imu0(iv))
           uv = gmu(imuv(iv))

           do i_st=1,nstokes
              do ns=0,nsmax

                 expu0(ns,iv) = exp(-tau_atm(ns)/u0)
                 expuv(ns,iv) = exp(-tau_atm(ns)/uv)

!                add Lambertian
                 if(i_st.eq.1)then
                    reflmRay_ss(ns,i_st,iv) = &
                      reflmRay_ss(ns,i_st,iv) + &
                          alb_Lambert * &
                           expu0(ns,iv) * expuv(ns,iv) * &
                            u0 / pi
                 endif

              enddo
           enddo

        enddo

     endif

!---------------------------------------------------------------------
!    Add Fourier moments
!
!    R_sys_ns = sum { Rm_sys_ns ...  }
!                m
!---------------------------------------------------------------------

     call add_Fourier_moment( &
            m,          nsmax,      ns_elas,           &
            ngeomax,    maxstr,     nstokes,           &
            single_scattering, &
            imuv,       imu0,       phi,               &
            deg2rad,    pi,         dfac,              &
            Bmplus,     Bmmin,      EpDelta,           &
            EmDelta,    R1m_sys,    Rm_sys,            &
            R_sys_ns)

!    ss
     do iv=1,ngeomax
        do i_st=1,nstokes
              do ns=0,nsmax
                 reflRay_ss(ns,i_st,iv) = reflRay_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmRay_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmRay_ss(ns,i_st,iv)

                 reflCab_ss(ns,i_st,iv) = reflCab_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmCab_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmCab_ss(ns,i_st,iv)

                 reflRam_ss(ns,i_st,iv) = reflRam_ss(ns,i_st,iv) + &
                                            EpDelta(i_st) * &
                                             dfac * Bmplus(i_st,iv) * &
                                              reflmRam_ss(ns,i_st,iv) + &
                                            EmDelta(i_st) * &
                                             dfac * Bmmin(i_st,iv) * &
                                              reflmRam_ss(ns,i_st,iv)

!                reflectivity
                 reflRay_ss(ns,i_st,iv) = reflRay_ss(ns,i_st,iv) * pi/u0
                 reflCab_ss(ns,i_st,iv) = reflCab_ss(ns,i_st,iv) * pi/u0
                 reflRam_ss(ns,i_st,iv) = reflRam_ss(ns,i_st,iv) * pi/u0

              enddo
           enddo
     enddo

  enddo !end loop over Fourier moments

!-----------------------------------------------------------------------
! for testing only
!-----------------------------------------------------------------------

  do iv=1,ngeomax
     do i_st=1,nstokes
        if(i_st.eq.1)then
           do ns=0,nsmax
              refl(ns,i_st,iv) = R_sys_ns(ns,iv)
           enddo
        else
           do ns=0,nsmax
              refl(ns,i_st,iv) = 0.d0
           enddo
        endif
     enddo
  enddo

  return
  end subroutine dasray_intf

max_stars_repo_name: rutgervandeelen/das

max_issues_repo_name: rutgervandeelen/das

max_forks_repo_name: rutgervandeelen/das


text: Eyad doesnt live in Davis anymore, but he has fond memories of the people and places associated with the town.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


! ::: -----------------------------------------------------------

      subroutine hypfill(adv,adv_l1,adv_l2,adv_l3,adv_h1,adv_h2,adv_h3, &
                         domlo,domhi,delta,xlo,time,bc) bind(C, name="hypfill")

      use amrex_fort_module, only : rt => amrex_real
      use meth_params_module

      implicit none
      include 'AMReX_bc_types.fi'
      integer adv_l1,adv_l2,adv_l3,adv_h1,adv_h2,adv_h3
      integer bc(3,2,*)
      integer domlo(3), domhi(3)
      real(rt) delta(3), xlo(3), time
      real(rt) adv(adv_l1:adv_h1,adv_l2:adv_h2,adv_l3:adv_h3,NVAR)

      integer n

      do n = 1,NVAR
          call filcc(adv(adv_l1,adv_l2,adv_l3,n), &
               adv_l1,adv_l2,adv_l3,adv_h1,adv_h2,adv_h3, &
               domlo,domhi,delta,xlo,bc(1,1,n))
      enddo

      end subroutine hypfill

! ::: -----------------------------------------------------------

      subroutine denfill(den,den_l1,den_l2,den_l3,den_h1,den_h2,den_h3, &
                         domlo,domhi,delta,xlo,time,bc) bind(C, name="denfill")

      use amrex_fort_module, only : rt => amrex_real
      implicit none
      include 'AMReX_bc_types.fi'
      integer den_l1,den_l2,den_l3,den_h1,den_h2,den_h3
      integer bc(3,2,*)
      integer domlo(3), domhi(3)
      real(rt) delta(3), xlo(3), time
      real(rt) den(den_l1:den_h1,den_l2:den_h2,den_l3:den_h3)

      call filcc(den,den_l1,den_l2,den_l3,den_h1,den_h2,den_h3,domlo,domhi,delta,xlo,bc)

      end subroutine denfill

! ::: -----------------------------------------------------------

      subroutine generic_fill(var,var_l1,var_l2,var_l3,var_h1,var_h2,var_h3, &
                              domlo,domhi,delta,xlo,time,bc) bind(C, name="generic_fill")

      use amrex_fort_module, only : rt => amrex_real
      implicit none
      include 'AMReX_bc_types.fi'
      integer var_l1,var_l2,var_l3,var_h1,var_h2,var_h3
      integer bc(3,2,*)
      integer domlo(3), domhi(3)
      real(rt) delta(3), xlo(3), time
      real(rt) var(var_l1:var_h1,var_l2:var_h2,var_l3:var_h3)

      call filcc(var,var_l1,var_l2,var_l3,var_h1,var_h2,var_h3,domlo,domhi,delta,xlo,bc)

      end subroutine generic_fill


text:
Ladle & Lotus offers holistic nutrition consultations, health education classes, and private yoga sessions.

Holistic Nutrition Consultations approach food as a healing modality.  Based on your interests and health goals this might take the form of: diet analysis and recommendations based on your metabolic assessment; assistance with meal planning; augmentation of recipes to make them healthier; individualized supplement recommendations; tips on label reading and grocery purchases; and strategies for implementing a seasonal and whole foods diet.

Health Education classes are held periodically at Kaya Yoga Kaya Yoga Studio.  Past classes have included: Detox & Digestive Tuneup; Immune Health; and  Weight Management.

Private Yoga sessions can be done with individuals, groups of friends, families, or at the workplace.  Private sessions honor the unique needs of the individual or group and tailors the session accordingly.

Background
Amy first began studying holistic health in 2004 while living in India.  Since then she has continued to expand her knowledge in the field by studying Ayurveda, massage therapy, yoga, nutrition, and integrative health more indepth.  In August 2011, she completed a 75credit Master’s degree in Holistic Health Education, with a specialization in nutrition, from John F. Kennedy University.

20110911 17:12:02 nbsp Ive taken Amys yoga classes, both at work (as part of the Wellness Program) and also at Kaya Yoga Studio.  Its been a great experience!  Shes knowledgeable about both yoga AND human anatomy (which isnt always true with yoga teachers).  She can also readily address concerns about past injuries affecting how one approaches yoga.  I have wrist and lower back issues but was able to feel comfortable in her yoga classes, knowing she could make suggestions that helped me adjust.   She is also a natural teacher who is friendly, accessible, and most importantly, truly passionate about what she teaches, someone who continues to be assiduous and humble in her own practice.  Having experienced a couple of yoga instructors who were (unconsciously) judgmental and/or bossy, I am really happy to have found Amy as my teacher.   It is my feeling that anyone can benefit from her instruction, but I particularly recommend her classes to beginners who may feel intimidated about starting yoga, or to those who feel they havent found a good instructor yet.                Users/ErikoMiura

20110912 14:59:37 nbsp Amy is a wonderful resource for personalized nutrition information. She is flexible and sincerely invested in her clients success in leading a healthy and fulfilling life through activity and a diet that works for them. She is incredibly knowledgeable in many disciplines and schools of thought within nutrition. She can speak to a practical, yet holistic approach to overall health. Despite moving out of town, Amy has still kept in contact and made sure that I am doing well in my transition. She is fantastic, full of wonderful information, and I happily have recommended her to many of my friends! Users/KCollins

20110917 18:12:09 nbsp Im fortunate enough to have had Amy as a yoga instructor at my workplace in Davis. Ive practiced yoga off & on for about 8 years (okay, more off than on), but Amys instruction has really given me a lot more enthusiasm for my practice. Ive learned more from her than from any of the other yoga instructors/studios Ive been to in town  shes got a strong understanding of anatomy, which is really helpful in adjusting poses or focusing practices to accommodate injuries or other physical limitations. But more importantly for me, Amy obviously takes a lot of enjoyment in teaching, and her classes are infused with humor. She also seems to be a good instructor for new or returning yoga students, based on the interactions Ive seen in classes. Users/TrayBiasiolli

20120105 14:41:59 nbsp I have learned so much about myself having Amy as my yoga instructor. Under her instruction I have been able to accomplish the unexpected and to accept and embrace my own physical (and mental) limitations. She has taught me to breathe. I can’t recommend her enough.

Her webpage also offers delicious and healthy recipes.
 Users/AngelaWhite


! { dg-do run }
! Test (re)allocation on assignment of scalars
!
! Contributed by Paul Thomas  <pault@gcc.gnu.org>
!
  call test_real
  call test_derived
  call test_char1
  call test_char4
  call test_deferred_char1
  call test_deferred_char4
contains
  subroutine test_real
    real, allocatable :: x
    real :: y = 42
    x = 42.0
    if (x .ne. y) call abort
    deallocate (x)
    x = y
    if (x .ne. y) call abort
  end subroutine
  subroutine test_derived
    type :: mytype
      real :: x
      character(4) :: c
    end type
    type (mytype), allocatable :: t
    t = mytype (99.0, "abcd")
    if (t%c .ne. "abcd") call abort
  end subroutine
  subroutine test_char1
    character(len = 8), allocatable :: c1
    character(len = 8) :: c2 = "abcd1234"
    c1 = "abcd1234"
    if (c1 .ne. c2) call abort
    deallocate (c1)
    c1 = c2
    if (c1 .ne. c2) call abort
  end subroutine
  subroutine test_char4
    character(len = 8, kind = 4), allocatable :: c1
    character(len = 8, kind = 4) :: c2 = 4_"abcd1234"
    c1 = 4_"abcd1234"
    if (c1 .ne. c2) call abort
    deallocate (c1)
    c1 = c2
    if (c1 .ne. c2) call abort
  end subroutine
  subroutine test_deferred_char1
    character(:), allocatable :: c
    c = "Hello"
    if (c .ne. "Hello") call abort
    if (len(c) .ne. 5) call abort
    c = "Goodbye"
    if (c .ne. "Goodbye") call abort
    if (len(c) .ne. 7) call abort
! Check that the hidden LEN dummy is passed by reference
    call test_pass_c1 (c)
    if (c .ne. "Made in test!") print *, c
    if (len(c) .ne. 13) call abort
  end subroutine
  subroutine test_pass_c1 (carg)
    character(:), allocatable :: carg
    if (carg .ne. "Goodbye") call abort
    if (len(carg) .ne. 7) call abort
    carg = "Made in test!"
  end subroutine
  subroutine test_deferred_char4
    character(:, kind = 4), allocatable :: c
    c = 4_"Hello"
    if (c .ne. 4_"Hello") call abort
    if (len(c) .ne. 5) call abort
    c = 4_"Goodbye"
    if (c .ne. 4_"Goodbye") call abort
    if (len(c) .ne. 7) call abort
! Check that the hidden LEN dummy is passed by reference
    call test_pass_c4 (c)
    if (c .ne. 4_"Made in test!") print *, c
    if (len(c) .ne. 13) call abort
  end subroutine
  subroutine test_pass_c4 (carg)
    character(:, kind = 4), allocatable :: carg
    if (carg .ne. 4_"Goodbye") call abort
    if (len(carg) .ne. 7) call abort
    carg = 4_"Made in test!"
  end subroutine
end


text:
Name: Dave Nguyen

About Me: Hi Everyone! I am a graduate student at UC Berkeley, but I currently live in Davis. As such, you can find my pretty much everyday riding on the morning to and evening back Davis Berkeley Shuttle. I usually sit in the back by the table so that I can read without falling asleep. Youll usually find me sleeping :\. I actually JUST moved here and wouldnt mind meeting some of you nice Davis folk, so please feel free to say hello!


!====================================================================
!
!  Define common block with all general parameters used by MadFKS
!  See their definitions in the file FKS_params.dat.
!
!====================================================================
module FKSParams
  character(len=64), parameter ::  paramFileName='FKS_params.dat'
  integer,parameter :: maxContribsSelected=100, &
                       maxCouplingsSelected=100, &
                       maxContribType=15, &
                       maxCouplingTypes=20
  real*8 :: IRPoleCheckThreshold,Virt_fraction, PrecisionVirtualAtRunTime,Min_virt_fraction
  integer  :: NHelForMCoverHels,VetoedContributionTypes(0:maxContribsSelected), &
              SelectedContributionTypes(0:maxContribsSelected),QED_squared_selected, &
              SelectedCouplingOrders(maxCouplingTypes,0:maxCouplingsSelected), &
              QCD_squared_selected
  logical :: separate_flavour_configs,IncludeBornContributions,use_poly_virtual

contains

  subroutine FKSParamReader(filename, printParam, force)
    ! Reads the file 'filename' and sets the parameters found in that file.
    implicit none
    logical, save :: HasReadOnce=.False.,paramPrinted=.false.
    logical :: force,couldRead,printParam
    character(*) :: filename
    CHARACTER(len=64) :: buff, buff2, mode
    include "orders.inc"
    integer :: i,j
    couldRead=.False.
    if (HasReadOnce.and..not.force) then
       goto 901
    endif
! Make sure to have default parameters if not set in the FKSParams.dat card
! (if it is an old one for instance)
    call DefaultFKSParam()
! Overwrite the default parameters from file:
    open(68, file=fileName, err=676, action='READ')
    do
       read(68,*,end=999) buff
       if(index(buff,'#').eq.1) then
          if (buff .eq. '#IRPoleCheckThreshold') then
             read(68,*,end=999) IRPoleCheckThreshold
             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
                stop 'IRPoleCheckThreshold must be >= -1.0d0.'
             endif
          elseif (buff .eq. '#PrecisionVirtualAtRunTime') then
             read(68,*,end=999) PrecisionVirtualAtRunTime
             if (IRPoleCheckThreshold .lt. -1.01d0 ) then
                stop 'PrecisionVirtualAtRunTime must be >= -1.0d0.'
             endif
          else if (buff .eq. '#NHelForMCoverHels') then
             read(68,*,end=999) NHelForMCoverHels
             if (NHelForMCoverHels .lt. -1) then
                stop 'NHelForMCoverHels must be >= -1.'
             endif
          else if (buff .eq. '#QCD^2==') then
             read(68,*,end=999) QCD_squared_selected
             if (QCD_squared_selected .lt. -1) then
                stop 'QCD_squared_selected must be >= -1.'
             endif
          else if (buff .eq. '#QED^2==') then
             read(68,*,end=999) QED_squared_selected
             if (QED_squared_selected .lt. -1) then
                stop 'QED_squared_selected must be >= -1.'
             endif
          else if (buff .eq. '#VirtualFraction') then
             read(68,*,end=999) Virt_fraction
             if (Virt_fraction .lt. 0 .or. virt_fraction .gt.1) then
                stop 'VirtualFraction should be a fraction between 0 and 1'
             endif
          else if (buff .eq. '#MinVirtualFraction') then
             read(68,*,end=999) Min_Virt_fraction
             if (min_virt_fraction .lt. 0 .or. min_virt_fraction .gt.1) then
                stop 'VirtualFraction should be a fraction between 0 and 1'
             endif
          else if (buff .eq. '#SeparateFlavourConfigurations') then
             read(68,*,end=999) separate_flavour_configs
          else if (buff .eq. '#UsePolyVirtual') then
             read(68,*,end=999) use_poly_virtual
          else if (buff .eq. '#VetoedContributionTypes') then
             read(68,*,end=999) VetoedContributionTypes(0)
             if (VetoedContributionTypes(0) .lt. 0.or. &
                  VetoedContributionTypes(0) .gt. maxContribsSelected) then
                write(*,*) 'VetoedContributionTypes length should be >= 0 and <=', &
                     maxContribsSelected
                stop 'Format error in FKS_params.dat.'
             endif
             read(68,*,end=999) (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
             do I=1,VetoedContributionTypes(0)
                if (VetoedContributionTypes(I).lt.1.or. &
                     VetoedContributionTypes(I).gt.maxContribType) then
                   write(*,*) 'VetoedContributionTypes must be >=1 and <=',maxContribType
                   stop 'Format error in FKS_params.dat.'
                endif
             enddo
             do I=VetoedContributionTypes(0)+1,maxContribsSelected
                VetoedContributionTypes(I)=-1
             enddo

          else if (buff .eq. '#SelectedContributionTypes') then
             read(68,*,end=999) SelectedContributionTypes(0)
             if (SelectedContributionTypes(0) .lt. 0 .or. &
                  SelectedContributionTypes(0) .gt. maxContribsSelected) then
                write(*,*) 'SelectedContributionTypes length should be >= 0 and <=', &
                     maxContribsSelected
                stop 'Format error in FKS_params.dat.'
             endif
             read(68,*,end=999) (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
             do I=1,SelectedContributionTypes(0)
                if (SelectedContributionTypes(I).lt.1.or. &
                     SelectedContributionTypes(I).gt.maxContribType) then
                   write(*,*) 'SelectedContributionTypes must be >=1 and <=',maxContribType
                   stop 'Format error in FKS_params.dat.'
                endif
             enddo
             do I=SelectedContributionTypes(0)+1,maxContribsSelected
                SelectedContributionTypes(I)=-1
             enddo
          else if (buff .eq. '#SelectedCouplingOrders') then
             read(68,*,end=999) SelectedCouplingOrders(1,0)
             if (SelectedCouplingOrders(1,0) .lt. 0 .or. &
                  SelectedCouplingOrders(1,0) .gt. maxCouplingsSelected) then
                write(*,*) 'SelectedCouplingOrders length should be >= 0 and <=', &
                     maxCouplingsSelected
                stop 'Format error in FKS_params.dat.'
             endif
             do j = 2, maxCouplingTypes
                SelectedCouplingOrders(j,0) = SelectedCouplingOrders(1,0)
             enddo
             do j = 1, SelectedCouplingOrders(1,0)
                read(68,*,end=999) (SelectedCouplingOrders(i,j),i=1,nsplitorders)
                do i=nsplitorders+1,maxCouplingTypes
                   SelectedCouplingOrders(i,j)=-1
                enddo
             enddo
          else
             write(*,*) 'The parameter name ',buff(2:),'is not reckognized.'
             stop 'Format error in FKS_params.dat.'
          endif
       endif
    enddo
999 continue
    couldRead=.True.
    goto 998

676 continue
    write(*,*) 'ERROR :: MadFKS parameter file ',fileName, &
         ' could not be found or is malformed. Please specify it.'
    stop
    !   Below is the code if one desires to let the code continue with
    !   a non existing or malformed parameter file
    write(*,*) 'WARNING :: The file ',fileName,' could not be ', &
         ' open or did not contain the necessary information. The ', &
         ' default MadFKS parameters will be used.'
    call DefaultFKSParam()
    goto 998
998 continue

    if(printParam.and..not.paramPrinted) then
       write(*,*) &
            '==============================================================='
       if (couldRead) then
          write(*,*) 'INFO: MadFKS read these parameters from ',filename
       else
          write(*,*) 'INFO: MadFKS used the default parameters.'
       endif
       write(*,*) &
            '==============================================================='
       write(*,*) ' > IRPoleCheckThreshold      = ',IRPoleCheckThreshold
       write(*,*) ' > PrecisionVirtualAtRunTime = ',PrecisionVirtualAtRunTime
       if (SelectedContributionTypes(0).gt.0) then
          write(*,*) ' > SelectedContributionTypes = ', &
               (SelectedContributionTypes(I),I=1,SelectedContributionTypes(0))
       else
          write(*,*) ' > SelectedContributionTypes = All'
       endif
       if (VetoedContributionTypes(0).gt.0) then
          write(*,*) ' > VetoedContributionTypes   = ', &
               (VetoedContributionTypes(I),I=1,VetoedContributionTypes(0))
       else
          write(*,*) ' > VetoedContributionTypes   = None'
       endif
       if (QCD_squared_selected.eq.-1) then
          write(*,*) ' > QCD_squared_selected      = All'
       else
          write(*,*) ' > QCD_squared_selected      = ',QCD_squared_selected
       endif
       if (QED_squared_selected.eq.-1) then
          write(*,*) ' > QED_squared_selected      = All'
       else
          write(*,*) ' > QED_squared_selected      = ',QED_squared_selected
       endif
       if (SelectedCouplingOrders(1,0).gt.0) then
          do j=1,SelectedCouplingOrders(1,0)
             write(*,*) ' > SelectedCouplingOrders(',j,') = ', &
                  (SelectedCouplingOrders(i,j),i=1,nsplitorders)
          enddo
       else
          write(*,*) ' > SelectedCouplingOrders    = All'
       endif
       write(*,*) ' > NHelForMCoverHels         = ',NHelForMCoverHels
       write(*,*) ' > VirtualFraction           = ',Virt_fraction
       write(*,*) ' > MinVirtualFraction        = ',Min_virt_fraction
       write(*,*) ' > SeparateFlavourConfigs    = ',separate_flavour_configs
       write(*,*) ' > UsePolyVirtual            = ',use_poly_virtual
       write(*,*) &
            '==============================================================='
       paramPrinted=.TRUE.
    endif

    close(68)
    HasReadOnce=.TRUE.
901 continue
  end subroutine FKSParamReader

  subroutine DefaultFKSParam()
    ! Sets the default parameters
    implicit none
    integer i,j
    IRPoleCheckThreshold=1.0d-5
    NHelForMCoverHels=5
    PrecisionVirtualAtRunTime=1d-3
    Virt_fraction=1d0
    QED_squared_selected=-1
    QCD_squared_selected=-1
    Min_virt_fraction=0.005d0
    separate_flavour_configs=.false.
    use_poly_virtual=.true.
    IncludeBornContributions=.true.
    SelectedContributionTypes(0)=0
    VetoedContributionTypes(0)=0
    do i=1, maxContribsSelected
       SelectedContributionTypes(I)=-1
       VetoedContributionTypes(I)=-1
    enddo
    do j=1,maxCouplingTypes
       SelectedCouplingOrders(j,0) = 0
    enddo
    do j=1,maxCouplingsSelected
       do i=1,maxCouplingTypes
          SelectedCouplingOrders(i,j) = -1
       enddo
    enddo
  end subroutine DefaultFKSParam

end module FKSParams


text: Noah is new to Davis this town and Front Page this wiki. He looks forward to getting to know people on or around campus, particularly people interested in Life, Ultimate Frisbee ultimate, or Plant Biology plants.


text:
program integration_test
    use, intrinsic :: iso_fortran_env, only : dp => REAL64, error_unit
    use :: gauss_mod, only : gen_gauss_type
    implicit none
    type(gen_gauss_type) :: quad
    integer :: n
    real(kind=dp) :: a = -1.0_dp, b = 1.0_dp, eps = 1.0e-9_dp
    character(len=1024) :: buffer

    if (command_argument_count() >= 3) then
        call get_command_argument(1, buffer)
        read (buffer, fmt='(I10)') n
        call get_command_argument(2, buffer)
        read (buffer, fmt='(E25.15)') a
        call get_command_argument(3, buffer)
        read (buffer, fmt='(E25.15)') b
    else
        write (unit=error_unit, fmt='(A)') &
            '#error: two command line argument expected, n, a, and b'
    end if
    call quad%init_gen_gauss(n, eps)
    print '(F25.15)', quad%integrate(f, a, b)

contains

    function f(x) result(r)
        implicit none
        real(kind=dp), intent(in) :: x
        real(kind=dp) :: r
        r = cos(x**2)*exp(0.37_dp*x)/(x**2 + 17.0_dp)
    end function f

end program integration_test

max_stars_repo_name: stijnvanhoey/training-material

max_issues_repo_name: stijnvanhoey/training-material

max_forks_repo_name: stijnvanhoey/training-material