text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains 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 NTESS 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 $Log: paintf.f,v $
C Revision 1.3  2009/03/25 12:36:46  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.2  1998/03/24 20:52:25  gdsjaar
C Fixed problem in finding average value of variable on 8-node quad
C face.
C Updated error message to send problem to SEACAS@sandia.gov instead of
C 'code sponsor'
C
C Revision 1.1  1994/04/07 20:06:43  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.4  1993/08/04  15:59:06  gdsjaar
c Performed some optimization to speed up blot - special-cased the most
c common situation.
c
c Revision 1.3  1993/04/29  17:27:34  gdsjaar
c Removed amy gilkeys name from error message.
c
c Revision 1.2  1990/12/14  08:54:34  gdsjaar
c Added RCS Id and Log to all files
c
C=======================================================================
      SUBROUTINE PAINTF (CNTR0, CNTR1, VARNP, NLNKF, LINKF1,
     &   XN, YN, ZN, XF, YF, ZF)
C=======================================================================

C   --*** PAINTF *** (DETOUR) Paint face contour for a face
C   --   Modified by John H. Glick - 10/26/88
C   --   Written by Amy Gilkey - revised 10/22/87
C   --
C   --PAINTF paints a contour section for a face.  The contour algorithm
C   --assumes that the elements do not contain internal nodes.
C   --The element interpolation field is approximated by logically drawing
C   --lines from each node to the element center and, thusly, dividing the
C   --element into triangles.  Contour sections are then drawn by connecting
C   --the intersection points of the sub-element edges and the contour
C   --plane.
C   --
C   --Parameters:
C   --   CNTR0, CNTR1 - the contour values delimiting the contour section
C   --   VARNP - IN - the contour function values
C   --   NLNKF - IN - the number of nodes per face
C   --   LINKF1 - IN - the connectivity for the face
C   --   XN, YN, ZN - IN - the nodal coordinates
C   --   XF, YF, ZF - IN - the face center coordinates
C   --
C   --Common Variables:
C   --   Uses IS3DIM of /D3NUMS/

      PARAMETER (KDONE=0, KERR1=-1, KERR2=-2, KERR3=-3,
     &   KCORNR=1, KISIDE=2, KMIDDL=3, KOSIDE=4, KFINDS=5)

      COMMON /D3NUMS/ IS3DIM, NNPSUR, NUMNPF, LLNSET
      LOGICAL IS3DIM

      REAL VARNP(*)
      INTEGER LINKF1(NLNKF)
      REAL XN(*), YN(*), ZN(*)

      LOGICAL INTERP_BL
      LOGICAL INTERC, INTEQC
      LOGICAL ISDONE(8)
      REAL XPTS(24), YPTS(24)

      ICNVS (I) = MOD (I+NNPF-1, NNPF) + 1
C      --ICNVS converts out of range side numbers (1..NLNKF)

      INTERC (CN, C0, C1) = (C0 .LT. CN) .AND. (CN .LT. C1)
C      --INTERC returns true iff corner is within two contours
      INTEQC (CN, CNN, C0, C1) = (CN .EQ. CNN)
     &   .AND. ((CN .EQ. C0) .OR. (CN .EQ. C1))
C      --INTEQC returns true iff corners are equal to each other and a contour

      if (nlnkf .eq. 4) then
C -- Most common case for 3D
C -- Compute center function value and coordinates

        FM = (varnp(linkf1(1)) + varnp(linkf1(2)) + varnp(linkf1(3)) +
     *    varnp(linkf1(4))) / 4.0
         XM = XF
         YM = YF
         NTRI = NLNKF
         NNPF = NLNKF
         MAXPT = 12

       else IF (NLNKF .EQ. 3) THEN

C      --Special case - triangle element

         N = LINKF1(2)
         FM = VARNP(N)
         XM = XN(N)
         YM = YN(N)
         IFL = 0
         NTRI = 1
         NNPF = NLNKF

      ELSE IF ((.NOT. IS3DIM) .AND. (NLNKF .EQ. 8)) THEN

C      --Compute center function value and coordinates

         FM = -0.25 * (varnp(linkf1(1)) + varnp(linkf1(3)) +
     *      varnp(linkf1(5)) + varnp(linkf1(7))) + 0.5 *
     *     (varnp(linkf1(2)) + varnp(linkf1(4)) +
     *      varnp(linkf1(6)) + varnp(linkf1(8)))
         FM = FM
         XM = XF
         YM = YF
         NTRI = NLNKF
         NNPF = NLNKF
         MAXPT = 24

      ELSE IF ((.NOT. IS3DIM) .AND. (NLNKF .EQ. 9)) THEN

C      --Compute center function value and coordinates

         FM = VARNP(LINKF1(9))
         XM = XN(LINKF1(9))
         YM = YN(LINKF1(9))
         NTRI = NLNKF - 1
         NNPF = NLNKF - 1
         MAXPT = 24

      ELSE

C      --Compute center function value and coordinates

         FM = 0.0
         DO 120 K = 1, NLNKF
            FM = FM + VARNP(LINKF1(K))
  120    CONTINUE
         FM = FM / NLNKF
         XM = XF
         YM = YF
         NTRI = NLNKF
         NNPF = NLNKF
         MAXPT = 12
      END IF

C   --Special check to paint the triangle formed by two adjacent nodes
C   --equal to each other, the middle value, and CNTR0

      IF ((NTRI .GT. 1) .AND. (CNTR0 .EQ. FM)) THEN
         IFL = 0
         FLAST = VARNP(LINKF1(NNPF))
         DO 130 K = 1, NNPF
            FTHIS = VARNP(LINKF1(K))
            IF ((FLAST .EQ. FTHIS) .AND. (FTHIS .EQ. FM)) IFL = K
            FLAST = FTHIS
  130    CONTINUE

         IF (IFL .GT. 0) THEN

            XPTS(1) = XM
            YPTS(1) = YM
            N1 = LINKF1(IFL)
            XPTS(2) = XN(N1)
            YPTS(2) = YN(N1)
            N = LINKF1(ICNVS(IFL-1))
            XPTS(3) = XN(N)
            YPTS(3) = YN(N)
            CALL MPD2PG (3, XPTS, YPTS, 'S')
         END IF
      END IF

      CALL INILOG (NTRI, .FALSE., ISDONE)

      NPTS = 0
      DO 150 ISIDE = 1, NTRI

C      --Check each side of the face, looking for a corner in the
C      --contour area or a dividing line in the contour area

         IF (ISDONE(ISIDE)) GOTO 150

         IS = ISIDE
         N = LINKF1(IS)
         N1 = LINKF1(ICNVS(IS+1))
         ISTART = IS

         ISTATE = KDONE
         IF (INTERC (VARNP(N), CNTR0, CNTR1)) THEN
            ICNTR = -1
            ISTATE = KCORNR
         ELSE IF ((INTERP_BL (CNTR0, VARNP(N), VARNP(N1), PSI))
     &      .AND. (VARNP(N1) .GT. CNTR0)) THEN
            ICNTR = 0
            CNTR = CNTR0
            ISTATE = KISIDE
         END IF

  140    CONTINUE
         IF ((ISTATE .GT. KDONE) .AND. (NPTS .LT. MAXPT)) THEN

            IF (ISTATE .EQ. KCORNR) THEN

C            --Corner is within the contour area, goto side-1

               NPTS = NPTS + 1
               XPTS(NPTS) = XN(N)
               YPTS(NPTS) = YN(N)
               ISDONE(IS) = .TRUE.
               IS = ICNVS(IS-1)
               ISDONE(IS) = .TRUE.
               N = LINKF1(IS)
               N1 = LINKF1(ICNVS(IS+1))
               ISTATE = KFINDS
               IF (IS .EQ. ISTART) ISTATE = KDONE

            ELSE IF (ISTATE .EQ. KISIDE) THEN

C            --Side crossing starts a contour area delimiter, find
C            --direction and follow contour line

               NPTS = NPTS + 1
               XPTS(NPTS) = XN(N) * (1.-PSI) + XN(N1) * PSI
               YPTS(NPTS) = YN(N) * (1.-PSI) + YN(N1) * PSI
               ISDONE(IS) = .TRUE.

C            --Find the direction to follow this contour line through
C            --the face center

               ISTATE = KERR3
               IF (INTERP_BL (CNTR, VARNP(N1), FM, PSI)) THEN
                  IDIR = +1
                  I = N1
                  ISTATE = KMIDDL
               ELSE IF (INTERP_BL (CNTR, VARNP(N), FM, PSI)) THEN
                  IDIR = -1
                  I = N
                  ISTATE = KMIDDL
               END IF

            ELSE IF (ISTATE .EQ. KMIDDL) THEN

C            --Crossing through the middle, continue following line
C            --to the side or through the middle

               NPTS = NPTS + 1
               XPTS(NPTS) = XN(I) * (1.-PSI) + XM * PSI
               YPTS(NPTS) = YN(I) * (1.-PSI) + YM * PSI
               ISDONE(IS) = .TRUE.
               IS = ICNVS(IS+IDIR)
               ISDONE(IS) = .TRUE.
               N = LINKF1(IS)
               N1 = LINKF1(ICNVS(IS+1))

               ISTATE = KERR1
               IF (IDIR .EQ. +1) THEN
                  I = N1
               ELSE
                  I = N
               END IF
               IF (INTERP_BL (CNTR, VARNP(N), VARNP(N1), PSI)) THEN
                  ISTATE = KOSIDE
C               --If the corner is exactly on the contour range, check
C               --if the line is moving away from the corner; if so, go
C               --through the middle to prevent cycling on this corner
                  IF (INTERP_BL (CNTR, VARNP(I), FM, PSIX)) THEN
                     IF (((IDIR .EQ. +1) .AND. (VARNP(N) .EQ. CNTR))
     &                  .OR.
     &                  ((IDIR .EQ. -1) .AND. (VARNP(N1) .EQ. CNTR)))
     &                  THEN
                        PSI = PSIX
                        ISTATE = KMIDDL
                     END IF
                  END IF
               ELSE IF (INTERP_BL (CNTR, VARNP(I), FM, PSI)) THEN
                  ISTATE = KMIDDL
               END IF

            ELSE IF (ISTATE .EQ. KOSIDE) THEN

C            --Contour line is intersecting the side, ending the line;
C            --search for enclosing corners/lines

               NPTS = NPTS + 1
               XPTS(NPTS) = XN(N) * (1.-PSI) + XN(N1) * PSI
               YPTS(NPTS) = YN(N) * (1.-PSI) + YN(N1) * PSI
               ISDONE(IS) = .TRUE.
               IF (VARNP(N) .EQ. CNTR) ISDONE(ICNVS(IS-1)) = .TRUE.
               IF (VARNP(N1) .EQ. CNTR) ISDONE(ICNVS(IS+1)) = .TRUE.

               ISTATE = KFINDS
               IF (IS .EQ. ISTART) ISTATE = KDONE
               IF (VARNP(N) .EQ. CNTR) THEN
                  IF ((IDIR .EQ. -1) .OR.
     &               ((ICNTR .EQ. 0) .AND. (IDIR .EQ. +1))) THEN
C                  --If the line exits at the corner, and enclosing area,
C                  --adjust side to be adjacent side
                     IS = ICNVS(IS-1)
                     ISDONE(IS) = .TRUE.
                     N = LINKF1(IS)
                     N1 = LINKF1(ICNVS(IS+1))
                     IF (IS .EQ. ISTART) ISTATE = KDONE
                  END IF
               END IF
            END IF

            IF (ISTATE .EQ. KFINDS) THEN

C            --Search for corner or first side crossing; the side can be
C            --a CNTR0 crossing if the last crossing was a corner or a
C            --CNTR1 crossing; it can be a CNTR1 crossing if the last
C            --crossing was a corner or a CNTR0 crossing

               ISTATE = KERR2
               IF (INTERC (VARNP(N), CNTR0, CNTR1) .OR.
     &            INTEQC (VARNP(N), VARNP(N1), CNTR0, CNTR1)) THEN
                  ICNTR = MIN (-1, ICNTR-1)
                  ISTATE = KCORNR
               ELSE IF ((ICNTR .EQ. -2) .AND.
     &            INTERP_BL (CNTR0, VARNP(N), VARNP(N1), PSI)) THEN
                  ICNTR = 0
                  CNTR = CNTR0
                  ISTATE = KISIDE
               ELSE IF ((ICNTR .NE. 1) .AND.
     &            INTERP_BL (CNTR1, VARNP(N), VARNP(N1), PSI)) THEN
                  ICNTR = 1
                  CNTR = CNTR1
                  ISTATE = KISIDE
               ELSE IF ((ICNTR .NE. 0) .AND.
     &            INTERP_BL (CNTR0, VARNP(N), VARNP(N1), PSI)) THEN
                  ICNTR = 0
                  CNTR = CNTR0
                  ISTATE = KISIDE
               END IF
            END IF

            GOTO 140
         END IF

         IF (ISTATE .NE. KDONE) THEN
C         --State will be <0 if error occurred, >0 if too many points
            CALL PLTFLU
            WRITE (*, 10000) ISTATE, IS, ICNTR, IDIR,
     &         (VARNP(LINKF1(I)), I=1,NLNKF), CNTR0, CNTR1
10000        FORMAT (' PAINT ERROR - please contact',
     &         ' SEACAS@sandia.gov with this information:',
     &         /, 3X, 3I3, SP, I3,
     &         2X, SP, 4 (1X, E9.2), 2X, 2 (1X, E9.2), /)
         END IF

C      --Paint enclosed area

         IF (NPTS .GT. 2) THEN
            CALL MPD2PG (NPTS, XPTS, YPTS, 'S')
            NPTS = 0
         END IF
  150 CONTINUE

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
C Segment-flux
C   calculate average flux in annuli bounded by a circular segment
C
C Document origin and update of this routine
C   P. Alexander, MRAO, Cambridge, 1/11/93
C
*-

C define work array
       integer       nm, nb
       integer        msize, bsize
       parameter    (msize = 256*256)
       parameter    (bsize=1024)
       parameter    (nm=16)
       parameter    (nb=3)
       real*4        map_array( nm*msize + nb*bsize )
C error status
       integer       status

C perform standard initialization
       call anm_start( 0,nm,msize,nb,bsize,status )

C read input map: imap is the catalogue identifier
       call do_segment_flux(map_array,status)

C check status value and report an error to the user
       call cmd_err(status,'SEGMENT-FLUX','Failed ')

C finally perform standard shut-down
       call anm_end( status )

       end
C
C
*+ do_segment_flux

       subroutine do_segment_flux(map_array,status)
C      -------------------------------------------------
C
C Add flux in elliptical annuli as defined by the user
C
C Given:
C   map data
       real*4        map_array(*)
C Returned:
C   error status word
       integer       status
C
C The flux is calculated in annuli bounded by a specified segment.
C The annuli are defined by an eccentricity, position-angle of the
C major axis and a width of the annulus.
C
C Output is to the standard output file which may be plotted or printed
C as required.
*-

       include '../include/anmap_sys_pars.inc'
       include '/mrao/include/chrlib_functions.inc'

C Local variables:
C   mini redtape
       integer    minirt(8)
C   used to hold the date (day, month, year)
       integer    idate(3)
C   input map and data pointer
       integer    imap, ip_map
C   results output variables
       integer    iout, n
       real*4     mean_radius
       integer    uv_range(4), number_annuli
       real*4     uv_centre(2)
C   radius of annulus (arcsec), position angle of ellipse (degrees),
C   eccentricity of ellipse
       real*4     radius_annulus, position_angle, eccentricity
C .. and in pixel coordintes / radians
       real*4     pix_rad_annulus, rad_pos_angle
C .. bounding angles of segment (degrees)
       real*4     seg_a1, seg_a2
C   maximum number of annuli the user may specify,  flux in each region,
C   error in each region, number of points found in each annulus
       integer    max_annuli
       parameter (max_annuli = 100)
       integer    number_points(max_annuli)
       real*4     average_flux(max_annuli), area_error(max_annuli)
C   redtape infortion used locally
       integer    equatorial_projection
       parameter (equatorial_projection = 1)
C   projection information
       integer    iproj, data_type
       real*8     usamp, skew, epoch, prang
C   source name and program name
       character  source*40, program*8

C check status
       if (status.ne.0) return

C read in the map
       call map_getmap('Map : ','Default_Map','READ',imap,status)
       call map_alloc_in(imap,'DIRECT',map_array,ip_map,status)
       call mapcat_enqsr(imap,source,program,status)
       if (status.ne.0) then
         call cmd_err(status,'SEGMENT-FLUX',' ')
         goto 999
       end if

C check that this image is a sky-plane image
       call enmapj(iproj,usamp,skew,epoch,prang,status)
       call enminirt(minirt,status)
       if (iproj.eq.equatorial_projection) then
         call cmd_wrerr('RING-FLUX',
     *             'Unable to work on EQUATORIAL-PLANE projection')
         call cmd_wrerr('SEGMENT-FLUX','Use REPROJECT-MAP')
         goto 999
       end if

C find user defined parameters
       call io_geti('Number of annuli : ','20',number_annuli,status)
       if (number_annuli.gt.max_annuli) then
         call cmd_wrerr('SEGMENT-FLUX','Too many annuli > 100')
         goto 999
       end if
       radius_annulus = usamp*1.0
       call io_getr('Radius of annulus (arcsec) : ','*',radius_annulus,
     *           status)
C .. convert to pixels
       pix_rad_annulus = radius_annulus/usamp

C find uv-range to analyse and centre of ellipse
       call enredt(uv_range,data_type,status)
       call plot_getuv('UV-range : ','*',uv_range,status)
       call plot_getpos2(minirt,'Centre of Ellipse : ','*',
     *                   uv_centre,status)

C read eccentricity and position angle of the major axis
       call io_getr('Eccentricity of ellipse: ','1.0',
     *               eccentricity,status)
       call io_getr('Position-angle (N thro'' E) : ','0.0',
     *           position_angle,status)
       call io_getr('Lower-segment-angle (+ve from X-axis, degrees) : ',
     *              '0.0',seg_a1,status)
       call io_getr('Upper-segment-angle (+ve from X-axis, degrees) : ',
     *              '360.0',seg_a2,status)
       if (seg_a1.lt.0.0) seg_a1 = seg_a1 + 360.0
       if (seg_a2.lt.0.0) seg_a2 = seg_a2 + 360.0
C .. convert to radians
       rad_pos_angle = position_angle*3.14159265/180.0

C do the work
      call image_addseg(minirt,map_array(ip_map),uv_range,uv_centre,
     *                   rad_pos_angle,eccentricity,
     *                   seg_a1, seg_a2,
     *                   number_annuli, pix_rad_annulus,
     *                   average_flux,number_points,area_error,status)
C
C write this out to the results file
       call io_opefil(iout,general_results_file,'WRITE',0,status)
       if (status.ne.0) then
         call cmd_err(status,'SEGMENT-FLUX',
     *                 'Error accessing results file')
         goto 999
       end if
       call util_enqdat(idate)
       write (iout,*)'%ndata ',number_annuli
       write (iout,*)'%ncols 4'
       write (iout,*)'%title RING FLUX output on ',idate
       write (iout,*)'%source ',source(1:chr_lenb(source))
       write (iout,*)'%title-1 Mean Radius of Annulus (arcsec)'
       write (iout,*)'%title-2 Average Flux (arbitrary units)'
       write (iout,*)'%title-3 No. pixels averaged in annulus'
       write (iout,*)'%title-4 Fractional Error in Area'
       do n = 1,number_annuli
         mean_radius = (float(n)-0.5)*radius_annulus
         write (iout,*) mean_radius, average_flux(n),
     *                  number_points(n), area_error(n)
       end do
       close (iout)

999    call map_end_alloc(imap,map_array,status)
       end
C
C
*+ image_addseg

      subroutine image_addseg(minirt,data,uv,ellipse_centre,
     *                         ellipse_theta, ellipse_eccen,
     *                         seg_a1, seg_a2,
     *                         n_annuli, r_annuli,
     *                         flux,ndata,frac_area,status   )
C     --------------------------------------------------------
C
C To calculate the average flux in elliptical annuli
C
C Given:
C   mini redtape of image
       integer     minirt(8)
C   image data
       real*4      data(*)
C   range on image to analyse
       integer     uv(4)
C   centre of ellipse
       real*4      ellipse_centre(2)
C   eccentricity and position-angle of ellipse
       real*4      ellipse_theta, ellipse_eccen
C   angles of acceptance for segment (degrees)
       real*4      seg_a1, seg_a2
C   number of annuli
       integer     n_annuli
C   spacing of annuli along major axis
       real*4      r_annuli
C   value above which averaging is performed
       real*4      gate
C Returned:
C   results for each annulus
C     flux
       real*4      flux(n_annuli)
C     number of data points
       integer     ndata(n_annuli)
C     fraction of correct annular arear foundd by cell averaging
       real*4      frac_area(n_annuli)
C   error return code
       integer     status
C
C The distribution of flux on an image centred at a given location
C within elliptical annuli is calculated according to the following
C algorithm.
C
C Divide image into elliptical bins and evaluate the mean flux density
C within each bin.  All the flux within a cell contributes to the bin
C enclosing the centre of that cell.  Two quantities are provided to
C estimate the usefulness of the output; the number of cells which have
C been averaged per bin; and the fraction of the correct area which has
C been evaluated by these points per bin.  Data values .le. the level
C 'gate' are not counted as valid data points
*-
C [A.J. Fitt, MRAO, Cambridge. June 1988.]
C [to NORD and ANMAP by P. Alexander, MRAO, Cambridge. June 1988.]
C [improved interface and standard coordinated, PA. 4/4/91]
C
C Local variables
C  coordinate along ellipse
       real*4      a, b
C  counters and sums
      integer      iu, iv, n, uv_point(2)
      equivalence (uv_point(1), iu)
      equivalence (uv_point(2), iv)
      real*4       x, y, value, rsq, areaconst, annarea
      real*4       sin_th, cos_th, seg_a
C  constant pi
      real*4       pi
      data         pi/3.14159265/

C check status on entry
      if (status.ne.0) return

C calculate geometric values
      sin_th = sin(ellipse_theta)
      cos_th = cos(ellipse_theta)

C initialise arrays
      do n = 1,n_annuli
        flux(n) = 0.0
        ndata(n) = 0
      end do

      do iv = uv(3),uv(4),-1
        do iu = uv(1),uv(2)
C .. find values
          call iuvval2(minirt,data,uv_point,value,status)
C .. find bin (n) containing point (x,y)
          x = float(iu) - ellipse_centre(1)
          y = float(iv) - ellipse_centre(2)
          seg_a = 180.0*atan2(y,x)/3.14159265
          if (seg_a.lt.0.0) seg_a = seg_a + 360.0
          a = -x*sin_th + y*cos_th
          b = x*cos_th + y*sin_th
          rsq = a**2 + (b*ellipse_eccen)**2
          n = sqrt(rsq)/r_annuli + 1.0
C .. add flux to correct bin
          if ((n.le.n_annuli) .and. (value.gt.gate) .and.
     *         (seg_a.ge.seg_a1 .and. seg_a.le.seg_a2) ) then
            flux(n) = flux(n) + value
            ndata(n) = ndata(n) + 1
          end if
        end do
      end do

C normalize and calculate fractional areas
      areaconst = pi*(r_annuli**2)/ellipse_eccen
      do n = 1,n_annuli
        if (ndata(n) .gt. 0 ) then
          flux(n) = flux(n)/float(ndata(n))
        else
          flux(n) = 0.0
        end if
        annarea = areaconst*float(2*n-1)
        frac_area(n) = ndata(n)/annarea
      end do
      call cmd_err(status,'image_addseg',' ')
      end
C

ext: f

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
      subroutine axibtc (s,d)
c     implicit double precision (a-h,o-z)                                    dp
       use mod_parameters
      common/range/mft,mlt,lft,llt,nftm1
      common/vect4/
     1 py1(nelemg),py2(nelemg),py3(nelemg),py4(nelemg),
     2 pz1(nelemg),pz2(nelemg),pz3(nelemg),pz4(nelemg),
     3 ph1(nelemg),ph2(nelemg),ph3(nelemg),ph4(nelemg)
      common/vect5/
     1 cb11(nelemg),cb21(nelemg),cb31(nelemg),cb41(nelemg),
     2 cb12(nelemg),cb22(nelemg),cb32(nelemg),cb42(nelemg),
     3 cb13(nelemg),cb23(nelemg),cb33(nelemg),cb43(nelemg),
     4 cb14(nelemg),cb24(nelemg),cb34(nelemg),cb44(nelemg),
     5 cb15(nelemg),cb25(nelemg),cb35(nelemg),cb45(nelemg),
     6 cb16(nelemg),cb26(nelemg),cb36(nelemg),cb46(nelemg),
     7 cb17(nelemg),cb27(nelemg),cb37(nelemg),cb47(nelemg),
     8 cb18(nelemg),cb28(nelemg),cb38(nelemg),cb48(nelemg)
      dimension s(44,1),d(4,4,1)

!$OMP PARALLEL DO
      do 10 i=lft,llt
      cb11(i)=d(1,3,i)*ph1(i)+d(1,1,i)*py1(i)+d(1,4,i)*pz1(i)
      cb21(i)=d(2,3,i)*ph1(i)+d(2,1,i)*py1(i)+d(2,4,i)*pz1(i)
      cb31(i)=d(3,3,i)*ph1(i)+d(3,1,i)*py1(i)+d(3,4,i)*pz1(i)
      cb41(i)=d(4,3,i)*ph1(i)+d(4,1,i)*py1(i)+d(4,4,i)*pz1(i)
      cb12(i)=d(1,2,i)*pz1(i)+d(1,4,i)*py1(i)
      cb22(i)=d(2,2,i)*pz1(i)+d(2,4,i)*py1(i)
      cb32(i)=d(3,2,i)*pz1(i)+d(3,4,i)*py1(i)
      cb42(i)=d(4,2,i)*pz1(i)+d(4,4,i)*py1(i)
   10 continue
!$OMP PARALLEL DO
      do 20 i=lft,llt
      cb13(i)=d(1,3,i)*ph2(i)+d(1,1,i)*py2(i)+d(1,4,i)*pz2(i)
      cb23(i)=d(2,3,i)*ph2(i)+d(2,1,i)*py2(i)+d(2,4,i)*pz2(i)
      cb33(i)=d(3,3,i)*ph2(i)+d(3,1,i)*py2(i)+d(3,4,i)*pz2(i)
      cb43(i)=d(4,3,i)*ph2(i)+d(4,1,i)*py2(i)+d(4,4,i)*pz2(i)
      cb14(i)=d(1,2,i)*pz2(i)+d(1,4,i)*py2(i)
      cb24(i)=d(2,2,i)*pz2(i)+d(2,4,i)*py2(i)
      cb34(i)=d(3,2,i)*pz2(i)+d(3,4,i)*py2(i)
      cb44(i)=d(4,2,i)*pz2(i)+d(4,4,i)*py2(i)
   20 continue
!$OMP PARALLEL DO
      do 30 i=lft,llt
      cb15(i)=d(1,3,i)*ph3(i)+d(1,1,i)*py3(i)+d(1,4,i)*pz3(i)
      cb25(i)=d(2,3,i)*ph3(i)+d(2,1,i)*py3(i)+d(2,4,i)*pz3(i)
      cb35(i)=d(3,3,i)*ph3(i)+d(3,1,i)*py3(i)+d(3,4,i)*pz3(i)
      cb45(i)=d(4,3,i)*ph3(i)+d(4,1,i)*py3(i)+d(4,4,i)*pz3(i)
      cb16(i)=d(1,2,i)*pz3(i)+d(1,4,i)*py3(i)
      cb26(i)=d(2,2,i)*pz3(i)+d(2,4,i)*py3(i)
      cb36(i)=d(3,2,i)*pz3(i)+d(3,4,i)*py3(i)
      cb46(i)=d(4,2,i)*pz3(i)+d(4,4,i)*py3(i)
   30 continue
!$OMP PARALLEL DO
      do 40 i=lft,llt
      cb17(i)=d(1,3,i)*ph4(i)+d(1,1,i)*py4(i)+d(1,4,i)*pz4(i)
      cb27(i)=d(2,3,i)*ph4(i)+d(2,1,i)*py4(i)+d(2,4,i)*pz4(i)
      cb37(i)=d(3,3,i)*ph4(i)+d(3,1,i)*py4(i)+d(3,4,i)*pz4(i)
      cb47(i)=d(4,3,i)*ph4(i)+d(4,1,i)*py4(i)+d(4,4,i)*pz4(i)
      cb18(i)=d(1,2,i)*pz4(i)+d(1,4,i)*py4(i)
      cb28(i)=d(2,2,i)*pz4(i)+d(2,4,i)*py4(i)
      cb38(i)=d(3,2,i)*pz4(i)+d(3,4,i)*py4(i)
      cb48(i)=d(4,2,i)*pz4(i)+d(4,4,i)*py4(i)
   40 continue
!$OMP PARALLEL DO
      do 50 i=lft,llt
      s(1,i)=s(1,i)+py1(i)*cb11(i)+pz1(i)*cb41(i)+ph1(i)*cb31(i)
      s(2,i)=s(2,i)+pz1(i)*cb21(i)+py1(i)*cb41(i)
      s(3,i)=s(3,i)+pz1(i)*cb22(i)+py1(i)*cb42(i)
      s(4,i)=s(4,i)+py2(i)*cb11(i)+pz2(i)*cb41(i)+ph2(i)*cb31(i)
      s(5,i)=s(5,i)+py2(i)*cb12(i)+pz2(i)*cb42(i)+ph2(i)*cb32(i)
      s(6,i)=s(6,i)+py2(i)*cb13(i)+pz2(i)*cb43(i)+ph2(i)*cb33(i)
      s(7,i)=s(7,i)+pz2(i)*cb21(i)+py2(i)*cb41(i)
      s(8,i)=s(8,i)+pz2(i)*cb22(i)+py2(i)*cb42(i)
      s(9,i)=s(9,i)+pz2(i)*cb23(i)+py2(i)*cb43(i)
      s(10,i)=s(10,i)+pz2(i)*cb24(i)+py2(i)*cb44(i)
   50 continue
!$OMP PARALLEL DO
      do 60 i=lft,llt
      s(11,i)=s(11,i)+py3(i)*cb11(i)+pz3(i)*cb41(i)+ph3(i)*cb31(i)
      s(12,i)=s(12,i)+py3(i)*cb12(i)+pz3(i)*cb42(i)+ph3(i)*cb32(i)
      s(13,i)=s(13,i)+py3(i)*cb13(i)+pz3(i)*cb43(i)+ph3(i)*cb33(i)
      s(14,i)=s(14,i)+py3(i)*cb14(i)+pz3(i)*cb44(i)+ph3(i)*cb34(i)
      s(15,i)=s(15,i)+py3(i)*cb15(i)+pz3(i)*cb45(i)+ph3(i)*cb35(i)
      s(16,i)=s(16,i)+pz3(i)*cb21(i)+py3(i)*cb41(i)
      s(17,i)=s(17,i)+pz3(i)*cb22(i)+py3(i)*cb42(i)
      s(18,i)=s(18,i)+pz3(i)*cb23(i)+py3(i)*cb43(i)
      s(19,i)=s(19,i)+pz3(i)*cb24(i)+py3(i)*cb44(i)
      s(20,i)=s(20,i)+pz3(i)*cb25(i)+py3(i)*cb45(i)
   60 continue
!$OMP PARALLEL DO
      do 70 i=lft,llt
      s(21,i)=s(21,i)+pz3(i)*cb26(i)+py3(i)*cb46(i)
      s(22,i)=s(22,i)+py4(i)*cb11(i)+pz4(i)*cb41(i)+ph4(i)*cb31(i)
      s(23,i)=s(23,i)+py4(i)*cb12(i)+pz4(i)*cb42(i)+ph4(i)*cb32(i)
      s(24,i)=s(24,i)+py4(i)*cb13(i)+pz4(i)*cb43(i)+ph4(i)*cb33(i)
      s(25,i)=s(25,i)+py4(i)*cb14(i)+pz4(i)*cb44(i)+ph4(i)*cb34(i)
      s(26,i)=s(26,i)+py4(i)*cb15(i)+pz4(i)*cb45(i)+ph4(i)*cb35(i)
      s(27,i)=s(27,i)+py4(i)*cb16(i)+pz4(i)*cb46(i)+ph4(i)*cb36(i)
      s(28,i)=s(28,i)+py4(i)*cb17(i)+pz4(i)*cb47(i)+ph4(i)*cb37(i)
      s(29,i)=s(29,i)+pz4(i)*cb21(i)+py4(i)*cb41(i)
      s(30,i)=s(30,i)+pz4(i)*cb22(i)+py4(i)*cb42(i)
   70 continue
!$OMP PARALLEL DO
      do 80 i=lft,llt
      s(31,i)=s(31,i)+pz4(i)*cb23(i)+py4(i)*cb43(i)
      s(32,i)=s(32,i)+pz4(i)*cb24(i)+py4(i)*cb44(i)
      s(33,i)=s(33,i)+pz4(i)*cb25(i)+py4(i)*cb45(i)
      s(34,i)=s(34,i)+pz4(i)*cb26(i)+py4(i)*cb46(i)
      s(35,i)=s(35,i)+pz4(i)*cb27(i)+py4(i)*cb47(i)
      s(36,i)=s(36,i)+pz4(i)*cb28(i)+py4(i)*cb48(i)
   80 continue
      return
      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
program test_distribution_PRNG
    use stdlib_error, only : check
    use stdlib_kinds, only: int8, int16, int32, int64
    use stdlib_random, only : random_seed, dist_rand

    implicit none
    logical ::  warn = .true.

    call test_random_seed
    call test_random_rand_iint8
    call test_random_rand_iint16
    call test_random_rand_iint32
    call test_random_rand_iint64

    contains

    subroutine test_random_seed
        integer :: put, get, res(5)
        integer :: ans(5) = [-1859553078, -1933696596, -642834430, &
                             1711399314,  1548311463]
        integer :: i

        print *, ""
        print *, "Test random_seed"
        put = 135792468
        do i = 1, 5
            call random_seed(put,get)
            res(i) = get
            put = get
        end do
        call check(all(res == ans), msg="random seed test failed.",warn=warn)
    end subroutine test_random_seed

    subroutine test_random_rand_iint8
        integer :: put, get, i

        integer(int8) :: res(5), ans(5)=[118, -15, -72, 101, 70]

        print *, ""
        print *, "Test random_rand with kind int8"
        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int8)
        end do
        call check(all(res == ans), msg="random_rand with kind int8 test" &
                                     //" failed.", warn=warn)
    end subroutine test_random_rand_iint8

    subroutine test_random_rand_iint16
        integer :: put, get, i

        integer(int16) :: res(5), ans(5)=[30286, -3799, -18204, 25947, 18148]

        print *, ""
        print *, "Test random_rand with kind int16"
        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int16)
        end do
        call check(all(res == ans), msg="random_rand with kind int16 test" &
                                     //" failed.", warn=warn)
    end subroutine test_random_rand_iint16

    subroutine test_random_rand_iint32
        integer :: put, get, i

        integer(int32) :: res(5), ans(5)=[1984865646, -248954393, -1192993267, &
                                  1700514835, 1189401802]

        print *, ""
        print *, "Test random_rand with kind int32"
        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int32)
        end do
        call check(all(res == ans), msg="random_rand with kind int32 test" &
                                     //" failed.", warn=warn)
    end subroutine test_random_rand_iint32

    subroutine test_random_rand_iint64
        integer :: put, get, i

        integer(int64) :: res(5), ans(5)=[8524933037632333570_int64, &
                                  -1069250973542918798_int64,        &
                                  -5123867065024149335_int64,        &
                                  7303655603304982073_int64,         &
                                  5108441843522503546_int64]

        print *, ""
        print *, "Test random_rand with kind int64"
        put = 12345678
        call random_seed(put, get)
        do i = 1, 5
            res(i) = dist_rand(1_int64)
        end do
        call check(all(res == ans), msg="random_rand with kind int64 test" &
                                     //" failed.", warn=warn)
    end subroutine test_random_rand_iint64

end program test_distribution_PRNG

max_stars_repo_name: wyphan/stdlib

max_issues_repo_name: alberteuler/stdlib


text:
module hipfort_types
  use iso_c_binding
  implicit none

  type,bind(c) :: dim3
     integer(c_int) :: x,y,z
  end type dim3

end module hipfort_types

max_stars_repo_name: afanfa/hipfort

max_issues_repo_name: afanfa/hipfort

max_forks_repo_name: afanfa/hipfort


text:
<ompts:test>
<ompts:testdescription>Test which checks the omp single directive by controlling how often a directive is called in an omp single region.</ompts:testdescription>
<ompts:ompversion>2.0</ompts:ompversion>
<ompts:directive>omp single</ompts:directive>
<ompts:dependences>omp parallel private,omp flush</ompts:dependences>
<ompts:testcode>
      INTEGER FUNCTION <ompts:testcode:functionname>omp_single</ompts:testcode:functionname>()
        IMPLICIT NONE
        INTEGER i
		<ompts:orphan:vars>
        INTEGER nr_threads_in_single,nr_iterations,result
        COMMON /orphvars/ nr_threads_in_single,nr_iterations,result
		</ompts:orphan:vars>
        INCLUDE "omp_testsuite.f"
        nr_threads_in_single=0
        result=0
        nr_iterations=0
!$omp parallel
        DO i=0, LOOPCOUNT-1
		<ompts:orphan>
!$omp single
!$omp flush
          nr_threads_in_single = nr_threads_in_single + 1
!$omp flush
          nr_iterations = nr_iterations + 1
          <ompts:check>nr_threads_in_single = nr_threads_in_single - 1</ompts:check>
          <ompts:crosscheck>nr_threads_in_single = nr_threads_in_single + 1</ompts:crosscheck>
          result = result + nr_threads_in_single
!$omp end single
		</ompts:orphan>
        END DO
!$omp end parallel
        IF ( result .EQ. 0 .AND. nr_iterations .EQ. LOOPCOUNT ) THEN
           <testfunctionname></testfunctionname> = 1
        ELSE
           <testfunctionname></testfunctionname> = 0
        END IF
      END FUNCTION
</ompts:testcode>
</ompts:test>

max_stars_repo_name: Alkaid-Benetnash/openmp


text:
Hey, I saw your comment you left over on Taqueria Guadalajara. Have you tried Taqueria Davis?  Users/WilliamLewis

20110518 19:57:27 nbsp I would concur wholeheartedly with the previous statement! Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C    %W% %G%
       subroutine mvdcx(ibus,iopt,iadr)
C * * *
C * * * THIS SUBROUTINE STORES THE REQUIRED DC OUTPUT QUANTITIES
C * * * FOR THE SPECIAL PLOTS OUTPUT.  IT IS CALLED BY NDCOUT.
C * * *
      include 'tspinc/params.inc'
      include 'tspinc/comn56.inc'                                              !dem
      include 'tspinc/mvn.inc'
*END  MVN
      include 'tspinc/link56.inc'
      include 'tspinc/room.inc'
      do 200 n = 1, 2*mvkt
      if(mvgen(n) .ne. ibus) go to 200
      if(mvcde(n) .ne. 'D') go to 200
      if(mvopt(n) .ne. iopt) go to 200
      mvadr(n) = madr1
      do 100 l = 1,icount+2
  100 vspce(madr1+l) = worksp(iadr+l)
      madr1 = madr1+icount+2
  200 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
! ******************************************************************************************************************************** !
! sedgem_lib.f90
!
! LIBRARY MODULE
! ******************************************************************************************************************************** !

MODULE sedgem_lib

  USE gem_cmn
  USE gem_util
  use gem_carbchem
  IMPLICIT NONE
  SAVE

  ! ****************************************************************************************************************************** !
  ! *** NAMELIST DEFINITIONS ***************************************************************************************************** !
  ! ****************************************************************************************************************************** !

  ! #### EDIT ADD AND/OR EXTEND NAME-LIST PARAMETER AND CONTROL OPTIONS ########################################################## !
  ! ------------------- RUN CONTROL ---------------------------------------------------------------------------------------------- !
  LOGICAL::ctrl_continuing                                       ! continuing run?
  NAMELIST /ini_sedgem_nml/ctrl_continuing
  REAL::start_year                                               ! Simulation start year [REAL]
  REAL::par_misc_t_runtime                                       ! run length (years)
  NAMELIST /ini_sedgem_nml/start_year,par_misc_t_runtime
  ! ------------------- PHYSICAL CONFIGURATION ----------------------------------------------------------------------------------- !
  REAL::par_sed_top_th                                           ! top ('well-mixed') sediment layer thickness (cm)
  REAL::par_sed_poros_det                                        ! detrital porosity (cm3(pore water) / cm3(sed))
  REAL::par_sed_poros_CaCO3                                      ! carbonate porosity in top layer (cm3(pore water) / cm3(sed))
  NAMELIST /ini_sedgem_nml/par_sed_top_th,par_sed_poros_det,par_sed_poros_CaCO3
  REAL::par_sed_Dmax_neritic                                     ! maximum depth of neritic (shallow water) sediments
  NAMELIST /ini_sedgem_nml/par_sed_Dmax_neritic
  LOGICAL::ctrl_sed_neritic_reef_force                           ! Force reef occurrence regardless of ocean depth?
  NAMELIST /ini_sedgem_nml/ctrl_sed_neritic_reef_force
  integer::par_n_sedcore_tot_min                                 ! Minimum (basic) number of sedcore layers
  integer::par_n_sedcore_tot_perky                               ! Number of dimensioned sedcore layers per ka of runtime
  NAMELIST /ini_sedgem_nml/par_n_sedcore_tot_min,par_n_sedcore_tot_perky
  integer::n_sed_tot                                             ! # sedimentary stack sub-layers
  integer::n_sed_tot_init                                        ! # initial sedimentary stack sub-layers filled
  integer::n_sed_tot_drop                                        ! # sedimentary stack sub-layers to drop off bottom
  NAMELIST /ini_sedgem_nml/n_sed_tot,n_sed_tot_init,n_sed_tot_drop
  ! ------------------- DIAGENESIS SCHEME: SELECTION ----------------------------------------------------------------------------- !
  character(len=63)::par_sed_diagen_CaCO3opt                     ! CaCO3 diagenesis scheme
  character(len=63)::par_sed_diagen_opalopt                      ! opal diagenesis scheme
  character(len=63)::par_sed_diagen_Corgopt                      ! Corg diagenesis scheme
  NAMELIST /ini_sedgem_nml/par_sed_diagen_CaCO3opt,par_sed_diagen_opalopt,par_sed_diagen_Corgopt
  ! ------------------- DIAGENESIS SCHEME: CONTROL ------------------------------------------------------------------------------- !
  LOGICAL::ctrl_sed_bioturb                                      ! Bioturbate sediment stack?
  logical::ctrl_sed_bioturb_Archer                               ! Use Archer et al. [2002] bioturbation scheme?
  NAMELIST /ini_sedgem_nml/ctrl_sed_bioturb,ctrl_sed_bioturb_Archer
  integer::par_n_sed_mix                                         ! maximum layer depth for bioturbation (# cm layers below surface)
  NAMELIST /ini_sedgem_nml/par_n_sed_mix
  REAL::par_sed_mix_k_sur_max                                    ! maximum surface bioturbation mixing rate (cm2 yr-1)
  REAL::par_sed_mix_k_sur_min                                    ! minimum surface bioturbation mixing rate (cm2 yr-1)
  NAMELIST /ini_sedgem_nml/par_sed_mix_k_sur_max,par_sed_mix_k_sur_min
  REAL::par_sed_fdet                                             ! prescribed (additional) flux of detrital material to the seds
  NAMELIST /ini_sedgem_nml/par_sed_fdet
  logical::ctrl_sed_noerosion
  NAMELIST /ini_sedgem_nml/ctrl_sed_noerosion
  logical::ctrl_sed_interface                                    ! CaCO3 interface dissolution?
  NAMELIST /ini_sedgem_nml/ctrl_sed_interface
  character(len=63)::opt_sed_CaCO3dislocation                    ! Location of CaCO3 dissolution
  NAMELIST /ini_sedgem_nml/opt_sed_CaCO3dislocation
  REAL::par_sed_CaCO3_fred                                       !
  REAL::par_sed_CaCO3_fblue                                      !
  NAMELIST /ini_sedgem_nml/par_sed_CaCO3_fred,par_sed_CaCO3_fblue
  logical::ctrl_sed_dyerestart                                   !
  NAMELIST /ini_sedgem_nml/ctrl_sed_dyerestart
  integer::par_sed_dyerestart_n                                  !
  NAMELIST /ini_sedgem_nml/par_sed_dyerestart_n
  ! ------------------- DIAGENESIS SCHEME: ORGANIC MATTER ------------------------------------------------------------------------ !
  LOGICAL::ctrl_sed_diagen_preserve_frac2                        ! Prevent frac2 from being remineralzied?
  NAMELIST /ini_sedgem_nml/ctrl_sed_diagen_preserve_frac2
  REAL::par_sed_diagen_fracCpres_ox                              ! Fractional POC burial -- oxic conditions
  REAL::par_sed_diagen_fracCpres_anox                            ! Fractional POC burial -- anoxic conditions
  REAL::par_sed_diagen_fracCpres_eux                             ! Fractional POC burial -- euxinic conditions
  NAMELIST /ini_sedgem_nml/par_sed_diagen_fracCpres_ox,par_sed_diagen_fracCpres_anox,par_sed_diagen_fracCpres_eux
  REAL::par_sed_diagen_fracC2Ppres_ox                            ! Fraction of P relative to C buried -- oxic
  REAL::par_sed_diagen_fracC2Ppres_anox                          ! Fraction of P relative to C buried -- anoxic
  REAL::par_sed_diagen_fracC2Ppres_eux                           ! Fraction of P relative to C buried -- euxinic
  NAMELIST /ini_sedgem_nml/par_sed_diagen_fracC2Ppres_ox,par_sed_diagen_fracC2Ppres_anox,par_sed_diagen_fracC2Ppres_eux
  ! ------------------- DIAGENESIS SCHEME: ARCHER 1991 --------------------------------------------------------------------------- !
  REAL::par_sed_archer1991_dissc                                 ! dissolution rate constant, units of 1/s
  REAL::par_sed_archer1991_dissn                                 ! dissolution rate order
  REAL::par_sed_archer1991_rc                                    ! organic degradation rate constant, 1/s
  NAMELIST /ini_sedgem_nml/par_sed_archer1991_dissc,par_sed_archer1991_dissn,par_sed_archer1991_rc
  integer::par_sed_archer1991_iterationmax                       ! loop limit in 'o2org' subroutine
  NAMELIST /ini_sedgem_nml/par_sed_archer1991_iterationmax
  NAMELIST /ini_sedgem_nml/par_sed_archer1991_iterationmax
  ! ------------------- CaCO3 PRODUCTION ----------------------------------------------------------------------------------------- !
  REAL::par_sed_CaCO3precip_sf                                   ! CaCO3 precipitation scale factor (abiotic)
  REAL::par_sed_CaCO3precip_exp                                  ! CaCO3 precipitation rate law lower (abiotic)
  NAMELIST /ini_sedgem_nml/par_sed_CaCO3precip_sf,par_sed_CaCO3precip_exp
  REAL::par_sed_reef_CaCO3precip_sf                              ! CaCO3 precipitation scale factor (corals)
  REAL::par_sed_reef_CaCO3precip_exp                             ! CaCO3 precipitation rate law power (corals)
  NAMELIST /ini_sedgem_nml/par_sed_reef_CaCO3precip_sf,par_sed_reef_CaCO3precip_exp
  logical::par_sed_reef_calcite                                  ! CaCO3 precipitation as calcite (otherwise aragonite)?
  NAMELIST /ini_sedgem_nml/par_sed_reef_calcite
  REAL::par_sed_CaCO3_abioticohm_min                             ! Min threshold for abiotic CaCO3 precipitation
  NAMELIST /ini_sedgem_nml/par_sed_CaCO3_abioticohm_min
  real::par_sed_poros_CaCO3_reef                                 ! reef CaCO3 porosity
  NAMELIST /ini_sedgem_nml/par_sed_poros_CaCO3_reef
  REAL::par_sed_CaCO3burial                                      ! prescribed CaCO3 production rate (mol cm-2 yr-1)
  NAMELIST /ini_sedgem_nml/par_sed_CaCO3burial
  ! ------------------- Corg PRODUCTION ------------------------------------------------------------------------------------------ !
  REAL::par_sed_Corgburial                                       ! prescribed Corg production rate (mol cm-2 yr-1)
  REAL::par_sed_Corgburial_Dd13C                                 ! POC d13C offset compared to the d13C of CaCO3
  NAMELIST /ini_sedgem_nml/par_sed_Corgburial,par_sed_Corgburial_Dd13C
  ! ------------------- TRACE METALS --------------------------------------------------------------------------------------------- !
  real::par_bio_red_CaCO3_LiCO3                                  ! Default CaCO3 Ca:Li ratio
  real::par_bio_red_CaCO3_LiCO3_alpha                            ! Partition coefficient (alpha)
  NAMELIST /ini_sedgem_nml/par_bio_red_CaCO3_LiCO3,par_bio_red_CaCO3_LiCO3_alpha
  ! ------------------- ISOTOPIC FRACTIONATION ----------------------------------------------------------------------------------- !
  CHARACTER(len=63)::opt_sed_foram_b_13C_delta                   ! Benthic foram 13C fractionation scheme ID string
  NAMELIST /ini_sedgem_nml/opt_sed_foram_b_13C_delta
  REAL::par_d7Li_LiCO3_epsilon                                   ! 7/6Li fractionation between Li and LiCO3
  NAMELIST /ini_sedgem_nml/par_d7Li_LiCO3_epsilon
  REAL::par_d44Ca_CaCO3_epsilon                                  ! 44/40Ca fractionation between Ca and CaCO3 -- corals
  NAMELIST /ini_sedgem_nml/par_d44Ca_CaCO3_epsilon
  REAL::par_d44Ca_abioticcal_epsilon0                            ! 44/40Ca fractionation between Ca and CaCO3 -- abiotic cal
  REAL::par_d44Ca_abioticarg_epsilon0                            ! 44/40Ca fractionation between Ca and CaCO3 -- abiotic arg
  NAMELIST /ini_sedgem_nml/par_d44Ca_abioticcal_epsilon0,par_d44Ca_abioticarg_epsilon0
  REAL::par_d44Ca_abioticcal_epsilondT                           ! T-dependence of 44/40Ca fractionation between Ca and CaCO3 (cal)
  REAL::par_d44Ca_abioticarg_epsilondT                           ! T-dependence of 44/40Ca fractionation between Ca and CaCO3 (arg)
  NAMELIST /ini_sedgem_nml/par_d44Ca_abioticcal_epsilondT,par_d44Ca_abioticarg_epsilondT
  ! ------------------- HYDROTHERMAL, OCEAN CRUSTAL WEATHERING, CLAY FORMATION --------------------------------------------------- !
  real::par_sed_hydroip_fLi                                      ! hydrothermal Li flux (mol yr-1)
  real::par_sed_hydroip_fLi_d7Li                                 ! hydrothermal Li flux d7Li (o/oo)
  NAMELIST /ini_sedgem_nml/par_sed_hydroip_fLi,par_sed_hydroip_fLi_d7Li
  real::par_sed_lowTalt_fLi_alpha                                ! Li low temperature alteration sink (mol yr-1) (Li/Ca normalized)
  real::par_sed_lowTalt_7Li_epsilon                              ! Li low temperature alteration sink 7Li epsilon (o/oo)
  NAMELIST /ini_sedgem_nml/par_sed_lowTalt_fLi_alpha,par_sed_lowTalt_7Li_epsilon
  real::par_sed_clay_fLi_alpha                                   ! Li clay formation sink (mol yr-1) (Li/Ca norm)
  real::par_sed_clay_7Li_epsilon                                 ! Li clay formation sink 7Li epsilon (o/oo)
  NAMELIST /ini_sedgem_nml/par_sed_clay_fLi_alpha,par_sed_clay_7Li_epsilon
  real::par_sed_hydroip_fCa                                      ! hydrothermal Ca flux (mol yr-1)
  real::par_sed_hydroip_fCa_d44Ca                                ! hydrothermal Ca flux d44Ca (o/oo)
  NAMELIST /ini_sedgem_nml/par_sed_hydroip_fCa,par_sed_hydroip_fCa_d44Ca
  real::par_sed_lowTalt_fCa_alpha                                ! Ca low-T alteration sink (mol yr-1) (Ca/Mg norm)
  real::par_sed_lowTalt_44Ca_epsilon                             ! Ca low-T alteration sink 44Ca epsilon (o/oo)
  NAMELIST /ini_sedgem_nml/par_sed_lowTalt_fCa_alpha,par_sed_lowTalt_44Ca_epsilon
  ! ------------------- MISC CONTROLS -------------------------------------------------------------------------------------------- !
  logical::ctrl_sed_forcedohmega_ca                              ! Ca-only adjustment for forced ocean saturation?
  NAMELIST /ini_sedgem_nml/ctrl_sed_forcedohmega_ca
  real::par_sed_ohmegamin                                        ! forced minimum saturation (calcite ohmega) anywhere
  real::par_sed_ohmegamin_flux                                   ! imposed sed->ocn flux (mol Ca cm-2 per time-step) for saturation
  NAMELIST /ini_sedgem_nml/par_sed_ohmegamin,par_sed_ohmegamin_flux
  ! ------------------- I/O: DIRECTORY DEFINITIONS ------------------------------------------------------------------------------- !
  CHARACTER(len=255)::par_indir_name                             !
  CHARACTER(len=255)::par_outdir_name                            !
  CHARACTER(len=255)::par_rstdir_name                            !
  NAMELIST /ini_sedgem_nml/par_indir_name,par_outdir_name,par_rstdir_name
  CHARACTER(len=127)::par_infile_name,par_outfile_name           !
  NAMELIST /ini_sedgem_nml/par_infile_name,par_outfile_name
  CHARACTER(len=127)::par_sed_topo_D_name                        ! Sediment water depth grid name
  CHARACTER(len=127)::par_sed_reef_mask_name                     ! Shallow water sediment (coral reef) mask name
  CHARACTER(len=127)::par_sedcore_save_mask_name                 ! Sediment core save mask name
  CHARACTER(len=127)::par_sed_mix_k_name                         ! Biodiffusion profile name
  NAMELIST /ini_sedgem_nml/par_sed_topo_D_name,par_sed_reef_mask_name,par_sedcore_save_mask_name,par_sed_mix_k_name
  CHARACTER(len=127)::par_output_years_file_0d                   ! file containing years for 0D output (to summary file)
  CHARACTER(len=127)::par_output_years_file_2d                   ! file containing years for 2D output (to netcdf/ascii)
  NAMELIST /ini_sedgem_nml/par_output_years_file_0d,par_output_years_file_2d
  ! ------------------- I/O: MISC ------------------------------------------------------------------------------------------------ !
  logical::ctrl_append_data                                      ! append data to output files on restart
  logical::ctrl_timeseries_output                                ! save timeseries output
  NAMELIST /ini_sedgem_nml/ctrl_append_data,ctrl_timeseries_output
  logical::ctrl_data_save_ascii                                  ! Save (sedcore) output in ascii format?
  logical::ctrl_data_save_sedcorenv                              ! Save sedcorenv output (ascii)?
  logical::ctrl_data_save_wtfrac                                 ! Report sediment data as a mass fraction?
  NAMELIST /ini_sedgem_nml/ctrl_data_save_ascii,ctrl_data_save_sedcorenv,ctrl_data_save_wtfrac
  logical::ctrl_misc_debug1                                      ! Debug level #1?
  logical::ctrl_misc_debug2                                      ! Debug level #2?
  logical::ctrl_misc_debug3                                      ! Debug level #3?
  logical::ctrl_misc_debug4                                      ! Debug level #4?
  NAMELIST /ini_sedgem_nml/ctrl_misc_debug1,ctrl_misc_debug2,ctrl_misc_debug3,ctrl_misc_debug4
  logical::ctrl_misc_report_err                                  ! report errors?
  NAMELIST /ini_sedgem_nml/ctrl_misc_report_err
  integer::par_misc_debug_i                                      ! i sediment coordinate for debug reporting
  integer::par_misc_debug_j                                      ! j sediment coordinate for debug reporting
  NAMELIST /ini_sedgem_nml/par_misc_debug_i,par_misc_debug_j
  real::par_sed_age_save_dt                                      ! threshold of accumulated time for sedcorenv save (yr)
  NAMELIST /ini_sedgem_nml/par_sed_age_save_dt
  ! ------------------- DATA SAVING: MISC ---------------------------------------------------------------------------------------- !
  LOGICAL::ctrl_ncrst                                          ! restart as netCDF format?
  NAMELIST /ini_sedgem_nml/ctrl_ncrst
  CHARACTER(len=127)::par_ncrst_name                           !
  NAMELIST /ini_sedgem_nml/par_ncrst_name
  CHARACTER(len=127)::par_ncsedcore_name
  NAMELIST /ini_sedgem_nml/par_ncsedcore_name
  ! ############################################################################################################################## !

  ! ****************************************************************************************************************************** !
  ! MODEL CONFIGURATION CONSTANTS - ARRAY DIMENSIONS
  ! ****************************************************************************************************************************** !

  ! grid dimensions
  INTEGER :: n_i, n_j
  ! grid properties array dimensions
  INTEGER :: n_phys_sed, n_opt_sed

  ! *** array index values ***
  ! sediment grid properties array indices
  INTEGER,PARAMETER::ips_lat                              = 01 ! latitude (degrees) [mid-point]
  INTEGER,PARAMETER::ips_lon                              = 02 ! longitude (degrees) [mid-point]
  INTEGER,PARAMETER::ips_dlat                             = 03 ! latitude (degrees) [width]
  INTEGER,PARAMETER::ips_dlon                             = 04 ! longitude (degrees) [width]
  INTEGER,PARAMETER::ips_latn                             = 05 ! latitude (degrees) [north edge]
  INTEGER,PARAMETER::ips_lone                             = 06 ! longitude (degrees) [east edge]
  INTEGER,PARAMETER::ips_D                                = 07 ! depth (m)
  INTEGER,PARAMETER::ips_A                                = 08 ! area (m2)
  INTEGER,PARAMETER::ips_rA                               = 09 ! reciprocal area (to speed up numerics)
  INTEGER,PARAMETER::ips_mask_sed                         = 10 ! sediment grid point mask (sediment = 1.0)
  INTEGER,PARAMETER::ips_mask_sed_reef                    = 11 ! reef grid point mask (reef = 1.0)
  INTEGER,PARAMETER::ips_mask_sed_muds                    = 12 ! shallow sediment grid point mask (muds = 1.0)
  INTEGER,PARAMETER::ips_poros                            = 13 ! sediment surface porosity
  INTEGER,PARAMETER::ips_mix_k0                           = 14 ! maximum (surface) sediment bioturbation mixing rate (cm2 yr-1)
  ! options - sediements
  INTEGER,PARAMETER::iopt_sed_save_diag_final             = 20 ! save final sediment data?
  INTEGER,PARAMETER::iopt_sed_save_diag                   = 21 ! save sediment diagnostics time-slice data?
  INTEGER,PARAMETER::iopt_sed_diagen_AltoasymSi           = 23 ! asymptotic [Si] dependence on %refrac/%opal?
  INTEGER,PARAMETER::iopt_sed_diagen_AltoKSi              = 24 ! KSi dependence on %refrac/%opal?

  ! *** look-up table constants ***
  ! CaCO3 (calcite)
  ! NOTE: following Ridgwell [2001]
  INTEGER,PARAMETER::lookup_i_D_min      = 0                   !
  INTEGER,PARAMETER::lookup_i_D_max      = 10                  !
  INTEGER,PARAMETER::lookup_i_dCO3_min   = -100                !
  INTEGER,PARAMETER::lookup_i_dCO3_max   = 100                 !
  INTEGER,PARAMETER::lookup_i_concO2_min = 4                   ! (equivalent to 200 umol kg-1)
  INTEGER,PARAMETER::lookup_i_concO2_max = 4                   ! (equivalent to 200 umol kg-1)
  INTEGER,PARAMETER::lookup_i_frac_min   = 1                   !
  INTEGER,PARAMETER::lookup_i_frac_max   = 10                  !
  INTEGER,PARAMETER::lookup_i_fCorg_min  = 0                   !
  INTEGER,PARAMETER::lookup_i_fCorg_max  = 50                  !
  REAL,PARAMETER::lookup_D_max      = 10000.0                  !
  REAL,PARAMETER::lookup_dCO3_max   = 100.0 * 1.0E-06          !
  REAL,PARAMETER::lookup_concO2_max = 200.0 * 1.0E-06          ! D(concO2) = 50 umol kg-1
  REAL,PARAMETER::lookup_frac_max   = 1.0                      !
  REAL,PARAMETER::lookup_fCorg_max  = 50.0 * 1.0E-06           !
  ! opal
  ! NOTE: following Ridgwell [2001]
  INTEGER,PARAMETER::lookup_i_opalpc_min       = 1             ! (2%)
  INTEGER,PARAMETER::lookup_i_opalpc_max       = 50            ! (100%)
  INTEGER,PARAMETER::lookup_i_concSi_min       = 0             ! (0 umol kg-1)
  INTEGER,PARAMETER::lookup_i_concSi_max       = 25            ! (250 umol kg-1)
  INTEGER,PARAMETER::lookup_i_T_min            = 270           ! (270 K)
  INTEGER,PARAMETER::lookup_i_T_max            = 280           ! (280 K)
  INTEGER,PARAMETER::lookup_i_KSi0_min         = 1             ! (0.010 yr-1 == 0.0 s-1)
  INTEGER,PARAMETER::lookup_i_KSi0_max         = 100           ! (1.000 yr-1 == 3.16E-08)
  INTEGER,PARAMETER::lookup_i_opaltorefrac_min = 0             ! (0.0 [refrac%/opal%])
  INTEGER,PARAMETER::lookup_i_opaltorefrac_max = 10            ! (10.0 [refrac%/opal%])
  REAL,PARAMETER::lookup_opalpc_max       = 1.0                ! D(opalpc) = 2%
  REAL,PARAMETER::lookup_concSi_max       = 250.0 * 1.0E-06    ! D(concSi) = 10 umol kg-1
  REAL,PARAMETER::lookup_T_max            = 280.0              ! D(T)      = 1 K
  REAL,PARAMETER::lookup_KSi0_max         = 1.000 / conv_yr_s  ! D(KSi0)   = 0.010 yr-1
  REAL,PARAMETER::lookup_opaltorefrac_max = 10.0               ! D(opaltorefrac) = 1.0

  ! ****************************************************************************************************************************** !
  ! GLOBAL VARIABLE AND RUN-TIME SET PARAMETER ARRAYS
  ! ****************************************************************************************************************************** !

  ! *** GRid parameters ***
  ! I/O - strings
  CHARACTER(len=63)::string_runid
  CHARACTER(len=127)::string_ncout2d                           ! name for netcdf output file
  ! netCDF and netCDF restart parameters
  CHARACTER(len=31)::string_rstid                              !
  CHARACTER(len=7) ::string_ncrunid                            !
  CHARACTER(len=254) ::string_ncrst                            !
  integer::ncrst_ntrec                                         ! count for netcdf datasets
  integer::ncrst_iou                                           ! io for netcdf restart
  !
  CHARACTER(len=254)::string_ncsedcorein                       !
  CHARACTER(len=254)::string_ncsedcoreout                      !
  integer::ncsedcore_ntrec                                     ! count for netcdf datasets
  integer::ncsedcore_iou                                       ! io for netcdf restart

  ! I/O - netCDF parameters
  integer::ntrec_sout                                          ! count for netcdf datasets
  integer::ntrec_siou                                          ! io for netcdf datasets
  ! flag for Archer sediment iteration (singular matrix) failure error
  logical::error_Archer = .FALSE.

  ! *** Array definitions ***
  ! bioturbation mixing rate array
  REAL,ALLOCATABLE,DIMENSION(:)::par_sed_mix_k                 ! bioturbation mixing rate profile array
  ! look-up tables
  REAL,ALLOCATABLE,DIMENSION(:,:,:,:)   :: lookup_sed_dis_cal  ! CaCO3 diagensis look-up table [Ridgwell, 2001]
  REAL,ALLOCATABLE,DIMENSION(:,:,:,:,:) :: lookup_sed_dis_opal ! opal diagenesis look-up table [Ridgwell, 2001]
  REAL,ALLOCATABLE,DIMENSION(:)::lookup_vec_D                  ! lookup table dimension vector
  REAL,ALLOCATABLE,DIMENSION(:)::lookup_vec_dco3               ! lookup table dimension vector
  REAL,ALLOCATABLE,DIMENSION(:)::lookup_vec_frac               ! lookup table dimension vector
  REAL,ALLOCATABLE,DIMENSION(:)::lookup_vec_fCorg              ! lookup table dimension vector
  ! allocatable 2-D sediment arrays
  real,ALLOCATABLE,DIMENSION(:,:,:)::phys_sed                  ! sediment 'physics' (mainly grid details)
  LOGICAL,ALLOCATABLE,DIMENSION(:,:)::sed_mask                 ! sediment mask (.TRUE. == sediment grid point exists)
  LOGICAL,ALLOCATABLE,DIMENSION(:,:)::sed_mask_reef            ! shallow water sediment mask - coral reefs
  LOGICAL,ALLOCATABLE,DIMENSION(:,:)::sed_mask_muds            ! shallow water sediment mask - muds
  REAL,ALLOCATABLE,DIMENSION(:,:,:,:)::sed                     ! the sediment layer stack
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_top                   ! top sedimentary layer
  REAL,ALLOCATABLE,DIMENSION(:,:)::sed_top_h                   ! top height of sedimentary column (cm)
  REAL,ALLOCATABLE,DIMENSION(:,:)::sed_top_INTdth              ! integrated sediment height change (cm)
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_fsed                  ! rain flux to sediments (mol cm-2 yr-1)
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_fdis                  ! sediment dissolution flux - solids tracers (mol cm-2 yr-1)
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sedocn_fnet               ! net sediment->ocean flux - ocean tracers (mol cm-2 yr-1)
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_carb                  ! carbonate chemistry overlying sediment surface
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_carbconst             ! carbonate chemistry constants
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_carbalk               ! carbonate chemistry alkalinity
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_carbisor              ! carbonate (carbon) isotopic properties array
  LOGICAL,ALLOCATABLE,DIMENSION(:,:)::sed_save_mask            ! sediment data save mask (.TRUE. == save sediment grid point)
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_fsed_OLD              !
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sed_fdis_OLD              !
  ! allocatable sedcoe arrays
  REAL,ALLOCATABLE,DIMENSION(:,:,:)::sedcore                   ! sedcore sediment layer stack (num sedcores x layers x variables)
  ! sediments - conversion
  real,DIMENSION(n_sed)::conv_sed_cm3_mol                      ! convert solids volume to number of moles
  real,DIMENSION(n_sed)::conv_sed_mol_cm3                      ! convert number of moles to solids volume
  real,DIMENSION(n_sed)::conv_sed_cm3_g                        ! convert solids volume to mass
  real,DIMENSION(n_sed)::conv_sed_g_cm3                        ! convert mass to solids volume
  real,DIMENSION(n_sed)::conv_sed_mask                         ! mask for which sediment tracers contribute to total solids volume
  ! misc
  LOGICAL,ALLOCATABLE,DIMENSION(:) :: opt_sed                  ! options arrays

  ! *** sedcore array definition ***
  !
  integer::nv_sedcore                                          ! number of sedcores
  integer::n_sedcore_tracer                                    !
  integer::n_sedcore_tot                                       ! number of sedcore (store) layers
  !
  type fieldsedcore
     integer::i
     integer::j
     real::ht                                                  ! sedcore height
     logical::save                                             ! save?
     real,allocatable,DIMENSION(:)::top                        ! sedcore top layer
     real,allocatable,DIMENSION(:,:)::lay                      ! sedcore layer (#layer, #tracer)
  end type fieldsedcore
  !
  type(fieldsedcore),DIMENSION(:),ALLOCATABLE::vsedcore_store  !

  ! *** MISC sediment parameters ***
  ! sediment mixing and layer configuration
  REAL::par_sed_interf_th                ! sediment interface dissolution layer thickness (cm)
  REAL::par_sed_dporos_top               ! compaction factor (surface sediment layer porosity compared to underlying sediment stack)
  REAL::par_sed_mix_zmix                 ! depth scale for bioturbation (cm) [Archer et al., 2002]
  REAL::par_sed_mix_c0_O2                ! half saturation constant of O2 for bioturbation (mol kg-1) [Archer et al., 2002]
  ! CaCO3 and Corg diagenesis
  REAL::par_caldis_k                     ! carbonate dissolution "rate constant"
  REAL::par_caldis_exp                   ! carbonate dissolution exponent
  REAL::par_sed_presfrac_FeO             ! fraction of sedimentating scavenged Fe and POM Fe flux preserved in the sediments
  REAL::par_sed_FeO_fdis                 ! forced dussolution flux of Fe from the sediments
  REAL::par_sed_diagen_fPOCmax           ! max POC rain flux in carbonate dissolution (umol cm-2 yr-1)
  ! opal diagenesis
  REAL::par_sed_opal_KSi0                ! base opal dissolution rate constant (intercept at zero opal rain rate)
  REAL::par_sed_opal_Sitoopalmax         ! asymptotic [Si] %refrac/%opal ratio max limite sediments
  ! ash tracing
  logical::par_sed_ashevent              ! ash event?
  real::par_sed_ashevent_fash            ! ash flux (g cm-2 kyr-1)
  !GHC 20/05/09 time-series saving parameters
  INTEGER                                        :: tstep_count
  REAL                                           :: tsteps_per_year
  INTEGER,PARAMETER                              :: n_output_years_max = 10000
  REAL, DIMENSION(:), ALLOCATABLE                :: output_years_0d
  REAL, DIMENSION(:), ALLOCATABLE                :: output_years_2d
  INTEGER , DIMENSION(:), ALLOCATABLE            :: output_tsteps_0d
  INTEGER , DIMENSION(:), ALLOCATABLE            :: output_tsteps_2d
  INTEGER                                        :: output_counter_0d
  INTEGER                                        :: output_counter_2d
  REAL                                           :: year
  INTEGER                                        :: year_int, year_remainder
  CHARACTER(LEN=11)                              :: year_text
  real::sed_age                                           ! sediment age (years)
  real::sed_time                                          ! sediment time counter (years)
  real::sed_time_save                                     ! sediment time save counter (years)

CONTAINS

  ! ****************************************************************************************************************************** !
  ! CALCULATE TOTAL SEDIMENT VOLUME
  ! NOTE: this function calculate the total volume of sediment tracers (i.e., not taking into account sediment porosity)
  !       -> a mask array <conv_sed_mask> is applied so that only those sediment tracers that actually have a solid volume
  !          contribute to the total silids volume
  !         (isotopic properties, and carbonate 'age', for instance, do not have any real volume and are not counted)
  FUNCTION fun_calc_sed_vol(dum_sed)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_sed_vol
    ! dummy arguments
    REAL,INTENT(in),DIMENSION(n_sed)::dum_sed
    ! return value
    fun_calc_sed_vol = sum(conv_sed_mask(:)*dum_sed(:))
  END FUNCTION fun_calc_sed_vol
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE TOTAL SEDIMENT MASS
  ! NOTE: similar to 'fun_calc_sed_vol' above, but converted to total sediment tracer mass, rather than left as volume
  FUNCTION fun_calc_sed_mass(dum_sed)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_sed_mass
    ! dummy arguments
    REAL,INTENT(in),DIMENSION(n_sed)::dum_sed
    ! return value
    fun_calc_sed_mass = sum(conv_sed_mask(:)*conv_sed_cm3_g(:)*dum_sed(:))
  END FUNCTION fun_calc_sed_mass
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE SEDIMENT POROSITY (SUB-SURFACE)
  ! NOTE: as per Zeebe and Zachos [2007]
  FUNCTION fun_calc_sed_poros(dum_frac_CaCO3)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_sed_poros
    ! dummy arguments
    REAL,INTENT(in)::dum_frac_CaCO3
    ! local variables
    real::loc_F
    ! calculate local constants
    loc_F = (par_sed_poros_CaCO3 - par_sed_poros_det)/(1.0 - par_sed_poros_CaCO3)
    ! return value
    fun_calc_sed_poros = (par_sed_poros_det + dum_frac_CaCO3*loc_F)/(1.0 + dum_frac_CaCO3*loc_F)
  END FUNCTION fun_calc_sed_poros
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE SEDIMENT POROSITY (NEAR-SURFACE)
  ! NOTE: Archer [1996]
  FUNCTION fun_calc_sed_poros_nsur(dum_frac_CaCO3,dum_sed_th)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_sed_poros_nsur
    ! dummy arguments
    REAL,INTENT(in)::dum_frac_CaCO3 !
    REAL,INTENT(in)::dum_sed_th     ! sediment thickness (cm)
    ! local variables
    real::loc_sed_poros
    real::loc_sed_poros_alpha
    ! calculate local constants
    loc_sed_poros = fun_calc_sed_poros(dum_frac_CaCO3)
    loc_sed_poros_alpha = fun_calc_sed_poros_alpha(dum_frac_CaCO3)
    ! return value
    fun_calc_sed_poros_nsur = loc_sed_poros - &
         & (loc_sed_poros_alpha*(1.0 - loc_sed_poros)/dum_sed_th)*(exp(-dum_sed_th/loc_sed_poros_alpha) - 1.0)
  END FUNCTION fun_calc_sed_poros_nsur
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE SEDIMENT POROSITY DEPTH SCALE (ALPHA)
  ! NOTE: Archer [1996]
  FUNCTION fun_calc_sed_poros_alpha(dum_frac_CaCO3)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_sed_poros_alpha
    ! dummy arguments
    REAL,INTENT(in)::dum_frac_CaCO3
    ! return value
    fun_calc_sed_poros_alpha = 0.25*dum_frac_CaCO3 + 3.0*(1.0 - dum_frac_CaCO3)
  END FUNCTION fun_calc_sed_poros_alpha
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE NEAR-SURFACE-TO-SUB-SURFACE COMPACTION RATIO
  FUNCTION fun_calc_r_sed_por(dum_frac_CaCO3,dum_sed_th)
    IMPLICIT NONE
    ! result variable
    REAL::fun_calc_r_sed_por
    ! dummy arguments
    REAL,INTENT(in)::dum_frac_CaCO3 !
    REAL,INTENT(in)::dum_sed_th     ! sediment thickness (cm)
    ! return value
    fun_calc_r_sed_por = (1.0 - fun_calc_sed_poros(dum_frac_CaCO3))/(1.0 - fun_calc_sed_poros_nsur(dum_frac_CaCO3,dum_sed_th))
  END FUNCTION fun_calc_r_sed_por
  ! ****************************************************************************************************************************** !

  ! ****************************************************************************************************************************** !
  ! CALCULATE SEDIMENT CORE-TOP DATA
  function fun_sed_coretop()
    IMPLICIT NONE
    ! result variable
    REAL,DIMENSION(n_sed,n_i,n_j)::fun_sed_coretop
    ! local variables
    INTEGER::i,j,l,is
    REAL::loc_sed_tot_wt
    REAL::loc_sed_tot_vol
    REAL,DIMENSION(n_sed,n_i,n_j)::loc_sed
    real::loc_tot,loc_frac,loc_standard

    ! *** calculate core-top composition ***

    ! initialize local array
    loc_sed(:,:,:) = 0.0
    ! loop through all sediment grid points and convert data
    ! NOTE: normailze solid sediment components to a mass fraction and multiply by 100.0
    ! NOTE: convert isotopic composition to delta notation in units of (o/oo)
    ! NOTE: screen-out non-wet (i,j) grid points
    ! NOTE: force zero wt% bulk locations to have null associated age and isotope tracer values
    DO i=1,n_i
       DO j=1,n_j
          IF (sed_mask(i,j)) THEN
             loc_sed_tot_wt  = fun_calc_sed_mass(sed_top(:,i,j))
             loc_sed_tot_vol = fun_calc_sed_vol(sed_top(:,i,j))
             DO l=1,n_l_sed
                is = conv_iselected_is(l)
                SELECT CASE (sed_type(is))
                case (par_sed_type_bio,par_sed_type_abio)
                   ! solid components
                   IF (ctrl_data_save_wtfrac) THEN
                      loc_sed(is,i,j) = 100.0*conv_sed_cm3_g(is)*sed_top(is,i,j)/loc_sed_tot_wt
                   ELSE
                      loc_sed(is,i,j) = sed_top(is,i,j)/loc_sed_tot_vol
                   ENDIF
                case (par_sed_type_POM)
                   ! particulate organic matter components
                   ! NOTE: mass (or volume) fraction has little meaning for the P,N,Fe,O2 components of POM,
                   !       so just calculate the ratio of these components with POC
                   if (loc_sed(is_POC,i,j) > const_real_nullsmall) loc_sed(is,i,j) = sed_top(is,i,j)/sed_top(is_POC,i,j)
                case (par_sed_type_age)
                   ! age
                   ! NOTE: extract core-top age
                   IF (sed_top(sed_dep(is),i,j) > const_real_nullsmall) THEN
                      loc_sed(is,i,j) = sed_top(is,i,j)/sed_top(sed_dep(is),i,j)
                   ELSE
                      loc_sed(is,i,j) = const_real_null
                   ENDIF
                case (n_itype_min:n_itype_max)
                   ! isotopes
                   IF (sed_top(sed_dep(is),i,j) > const_real_nullsmall) THEN
                      loc_tot  = sed_top(sed_dep(is),i,j)
                      loc_frac = sed_top(is,i,j)
                      loc_standard = const_standards(sed_type(is))
                      loc_sed(is,i,j) = fun_calc_isotope_delta(loc_tot,loc_frac,loc_standard,.FALSE.,const_nulliso)
                   else
                      loc_sed(is,i,j) = const_real_null
                   endif
                case default
                   ! everything else
                   ! NOTE: assume no normalization
                   loc_sed(is,i,j) = sed_top(is,i,j)
                end select
             END DO
          ENDIF
       END DO
    END DO
    ! set coretop sediment composition
    fun_sed_coretop(:,:,:) = loc_sed(:,:,:)

  END function fun_sed_coretop
  ! ****************************************************************************************************************************** !

END MODULE sedgem_lib

max_stars_repo_name: derpycode/cgenie

max_issues_repo_name: derpycode/cgenie

max_forks_repo_name: derpycode/cgenie


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE EL5(JJA,JJB,JA,JB,JC,JD,ICOLBREI)
!                                                                  *
!   --------------  SECTION METWO    SUBPROGRAM 11  -------------  *
!                                                                  *
!     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
!     OF TWO PARTICLE OPERATOR IN CASE :    N'1 = N1 (+-) 1        *
!                                           N'2 = N2 (+-) 1        *
!                                           N'3 = N3 (+-) 1        *
!                                           N'4 = N4 (+-) 1        *
!                                                                  *
!      SUBROUTINE CALLED: EL51,EL52,EL53                           *
!                                                                  *
!   Written by  G. Gaigalas                                        *
!   Transform to fortran 90/95 by G. Gaigalas       December 2012  *
!   The last modification made by G. Gaigalas       October  2017  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE m_C,          ONLY: NPEEL
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE el51_I
      USE el52_I
      USE el53_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: JJA,JJB,JA,JB,JC,JD,ICOLBREI
!-----------------------------------------------
      IF(NPEEL <= 3)RETURN
      IF(JB < JC) THEN
        CALL EL51(JJA,JJB,JA,JB,JC,JD,1,ICOLBREI)
      ELSE IF(JA > JD.AND.JB > JD) THEN
        CALL EL51(JJA,JJB,JC,JD,JA,JB,2,ICOLBREI)
      ELSE IF(JB > JC.AND.JB < JD.AND.JA < JC) THEN
        CALL EL52(JJA,JJB,JA,JC,JB,JD,1,ICOLBREI)
      ELSE IF(JB > JC.AND.JB > JD.AND.JA > JC) THEN
        CALL EL52(JJA,JJB,JC,JA,JD,JB,2,ICOLBREI)
      ELSE IF(JB > JC.AND.JB > JD.AND.JA < JC) THEN
        CALL EL53(JJA,JJB,JA,JC,JD,JB,1,ICOLBREI)
      ELSE IF(JB > JC.AND.JB < JD.AND.JA > JC) THEN
        CALL EL53(JJA,JJB,JC,JA,JB,JD,2,ICOLBREI)
      ELSE
        WRITE(99,100)
        STOP
      END IF
      RETURN
  100 FORMAT(5X,'ERRO IN EL5 ')
      END SUBROUTINE EL5

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
MODULE m_lorentzian_smooth

   USE m_constants

   IMPLICIT NONE

   !PARAMETER FOR LORENTZIAN SMOOTHING
   REAL,    PARAMETER :: cut              = 1e-8

   CONTAINS

   !This is essentially smooth out of m_smooth but with a lorentzian distribution
   SUBROUTINE lorentzian_smooth(e,f,sigma,n)

      INTEGER, INTENT(IN)    :: n
      REAL,    INTENT(INOUT) :: f(:)
      REAL,    INTENT(IN)    :: sigma
      REAL,    INTENT(IN)    :: e(:)

      REAL :: dx
      REAL :: f0(n), ee(n)
      INTEGER :: ie , je , j1 , j2 , numPoints

      f0 = f
      f = 0.0
      ee = 0.0
      dx = e(2)-e(1)

      DO ie =1,  n

         ee(ie) = 1/pi_const * sigma/dx * 1.0/((ie-1)**2+(sigma/dx)**2)

         IF ( ee(ie).LT.cut ) EXIT
      ENDDO
      numPoints = ie - 1

      DO ie = 1, n

         j1 = ie - numPoints + 1
         j1 = MERGE(1,j1,j1.LT.1)
         j2 = ie + numPoints - 1
         j2 = MERGE(n,j2,j2.GT.n)
         DO je = j1 , j2
            f(ie) = f(ie) + ee(IABS(je-ie)+1)*f0(je)
         ENDDO

      ENDDO

   END SUBROUTINE lorentzian_smooth

END MODULE m_lorentzian_smooth

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR