text:
        INTEGER FUNCTION GETMENU( ITEMCNT, DEFAULT, PROMPT, CHOICES )

C...............................................................
C Version "$Id: getmenu.f 1 2017-06-10 18:05:20Z coats $"
C EDSS/Models-3 I/O API.
C Copyright (C) 1992-2002 MCNC and Carlie J. Coats, Jr.,
C (C) 2003-2013 Baron Advanced Meteorological Systems,
C (C) 2007-2013 Carlie J. Coats, Jr., and
C (C) 2014 UNC Institute for the Environment.
C Distributed under the GNU LESSER GENERAL PUBLIC LICENSE version 2.1
C See file "LGPL.txt" for conditions of use.
C.........................................................................
C       function body starts at line  87
C
C       CALLS:
C       RETURNS:    user-selected item number from menu choices.
C
C  FUNCTION:
C
C       Display a menu to the screen; then display a prompt, get the
C       user's response and check that it is within range.
C       NOTE:  prompt and menu choice strings  should have
C       length <= 72 characters.  Ideally, number of items should be <= 18;
C       should be <= 999 and >= 1 in any case.
C
C  ARGUMENT LIST DESCRIPTION:
C
C    Input arguments:
C
C       ITEMCNT     number of items in menu. Ideally, < 20
C       DEFAULT     default menu choice
C       PROMPT      menu prompt
C       CHOICES     array of menu choices
C
C    Output arguments:  none
C
C  REVISION HISTORY:
C
C       Created  9/1989 by CJC.
C       Modified 6/1990 by CJC for ROM 2.2 -- error exit via EXWST
C       Modified 8/1990 by CJC for ROM 2.2 -- TEMP.LET; exits on CONTROL-Z.
C       Modified 2/1993 by CJC for CRAY
C       Modified 8/1996 by CJC to treat "!" as a terminator
C       Modified 1/1997 by CJC:  logs result
C       Revised  5/2003 by CJC:  factor through M3MSG2 to ensure flush()
C       of log-messages.
C       Revised 6/2003 by CJC:  factor through M3PROMPT to ensure flush()
C       of PROMPT for IRIX F90v7.4
C       Modified 03/2010 by CJC: F9x changes for I/O API v3.1
C       Modified 02/2015 by CJC for I/O API 3.2: USE M3UTILIO.
C       Fix MH violation of coding-standards:  check status IOS from  ENVYN!!
C***********************************************************************

        IMPLICIT NONE

C.......   ARGUMENTS:

        INTEGER      , INTENT(IN   ) :: ITEMCNT         !  number of choices
        INTEGER      , INTENT(IN   ) :: DEFAULT         !  default response
        CHARACTER*(*), INTENT(IN   ) :: PROMPT          !  prompt string
        CHARACTER*(*), INTENT(IN   ) :: CHOICES ( * )   !  array of choice strings

C.......   EXTERNAL FUNCTION:

        LOGICAL, EXTERNAL :: ENVYN

C.......   LOCAL VARIABLES:

        INTEGER         LENGTH          !  length of default as ASCII string
        INTEGER         I , M
        INTEGER         ANSWER
        INTEGER         ERRCNT
        INTEGER         IOS
        CHARACTER*32    BUFFER
        CHARACTER*8     FMTSTR
        CHARACTER*256   MESG

        LOGICAL, SAVE :: PROMPTON
        LOGICAL, SAVE :: FIRSTIME = .TRUE.

        CHARACTER*16, PARAMETER :: PNAME = 'GETMENU'

C...............................................................
C   begin body of GETMENU

        IF( FIRSTIME ) THEN

            PROMPTON = ENVYN( 'PROMPTFLAG', 'Prompt for input flag',
     &                      .TRUE., IOS )
            IF ( IOS .GT. 0 ) THEN
                CALL M3EXIT( PNAME,0,0,'Bad env vble "PROMPTFLAG"', 2 )
            END IF
            FIRSTIME = .FALSE.

        END IF

        IF( .NOT. PROMPTON ) THEN
            GETMENU = DEFAULT
            MESG = 'Using default response "' //
     &          TRIM( CHOICES( DEFAULT ) ) // '" for query:'
            CALL M3MSG2( MESG )
            MESG = '"' // TRIM( PROMPT ) // '"'
            CALL M3MSG2( MESG )
            RETURN
        END IF

        IF ( ITEMCNT .GT. 9999 ) THEN
            WRITE ( 6,91010 )
     &              'Maximum number of items:  ', 9999,
     &              'Number of requested items:', ITEMCNT
            CALL M3EXIT( 'GETMENU', 0, 0,
     &                   'Maximum item count exceeded.', 1 )
        ELSE IF ( ITEMCNT .LT. 1 ) THEN
            WRITE ( 6,91010 )
     &              'Minimum number of items:  ', 1,
     &              'Number of requested items:', ITEMCNT
            CALL M3EXIT( 'GETMENU', 0, 0,
     &                   'Minimum item count exceeded.', 1 )
        END IF

        I       =  1
        LENGTH  =  0
10      CONTINUE

            I       =  I * 10
            LENGTH  =  LENGTH  +  1
            IF ( I .LE. DEFAULT )  GO TO  10

        ERRCNT =  0

100     CONTINUE

        WRITE ( 6,* )
        WRITE ( 6,92010 )
        IF ( ITEMCNT .LE. 18 )  WRITE ( 6,92020 )

        DO  200  I = 1 , ITEMCNT

            WRITE ( 6,92030 ) I , TRIM( CHOICES ( I ) )

200     CONTINUE

        IF ( ITEMCNT .LE. 18 )  WRITE ( 6,92020 )
        WRITE ( 6,92010 )

        IF ( ITEMCNT .LT. 10 )THEN      !  select prompt-format:
            BUFFER = '( 2A, I1, A )'
        ELSE IF ( ITEMCNT .LT. 100 )THEN
            BUFFER = '( 2A, I2, A )'
        ELSE
            BUFFER = '( 2A, I4, A )'
        END IF

        WRITE( MESG, BUFFER ) TRIM( PROMPT ), ' [', DEFAULT, '] >> '
        CALL M3PROMPT( MESG, BUFFER, IOS )

        IF ( IOS .NE. 0 )  THEN
            GO TO 900
        ELSE IF ( BUFFER ( 1:1 )  .EQ. ' ' )  THEN
            GETMENU  =  DEFAULT
            MESG = 'Using default "' // TRIM( CHOICES( DEFAULT ) )//'"'
        ELSE
            WRITE( FMTSTR, 94010 ) LEN_TRIM( BUFFER )

            READ( BUFFER, FMTSTR, IOSTAT=IOS, ERR=400 ) ANSWER

            IF ( ANSWER .LT. 1  .OR.  ANSWER .GT. ITEMCNT )  THEN

                ERRCNT  =  ERRCNT + 1
                WRITE ( 6,91020 ) ANSWER , ITEMCNT , ERRCNT
                CALL M3FLUSH( 6 )
                IF ( ERRCNT .LT. 5 )  GO TO  100

                GO TO  500      !  max error count exceeded

            END IF

            GETMENU  =  ANSWER
            MESG = 'Using response "'// TRIM( CHOICES( ANSWER ) )// '"'
        END IF
        CALL M3MSG2( MESG )

        RETURN

400     CONTINUE        !  error in read from BUFFER

        ERRCNT  =  ERRCNT + 1
        WRITE ( 6,91030 )  ERRCNT
        CALL M3FLUSH( 6 )
        IF ( ERRCNT .LT. 5 )  GO TO  100

500     CONTINUE        !  more than 5 entry errors

        CALL M3EXIT( 'GETMENU', 0, 0,
     &               'Maximum error count exceeded', 2 )

900     CONTINUE        !  error in read from terminal

        ERRCNT  =  ERRCNT + 1
        WRITE ( 6,91030 )  ERRCNT
        CALL M3FLUSH( 6 )
        IF ( ERRCNT .LT. 5 )  GO TO  100
        WRITE ( 6,91040 ) IOS
        CALL M3EXIT( 'GETMENU', 0, 0, 'MAX ERRCNT exceeded', 1 )

C.................   end body of GETMENU  ..............................

C...........   Error and warning message formats..... 91xxx

91010   FORMAT ( 2 ( /10X , A, I12 ) )

91020   FORMAT ( /5X , '>>>  RESPONSE ERROR  <<<'      ,
     &         //10X , 'Your response'          , I6 , 2X,
     &                 'not in the range 1 ...' , I6 ,
     &         /10X  , 'This is error ' , I1 , ' of 5 errors allowed ' ,
     &         /10X  , 'Please try again.'            , / )

91030   FORMAT ( /5X , '>>>  RESPONSE ERROR  <<<'         ,
     &         //10X , 'Did not understand your response' ,
     &         /10X  , 'This is error ' , I1 , ' of 5 errors allowed ' ,
     &         /10X  , 'Please try again.'                , / )

91040   FORMAT ( //10X  , 'Error number = ' , I3 )

C...........   Informational (LOG) message formats... 92xxx

92010   FORMAT ( 1X , 79('*') )

92020   FORMAT ( 1X , '*' , T79 , '*' )

92030   FORMAT ( 1X , '*  ' , I2 , ':  ' , A , T79 , '*' )

C...........   Internal buffering formats............ 94xxx

94010   FORMAT( '( I', I3, ')' )

        END FUNCTION GETMENU

ext: f

max_stars_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_issues_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork

max_forks_repo_name: mdhillmancmcl/TheWorldAvatar-CMCL-Fork


text:
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine quadraticsens(ipkon,lakon,kon,nobject,dgdxglob,
     &   xinterpol,nnodes,ne,nk,nodedesiinv,objectset)
!
!     interpolation of the sensitivitites of the midnodes to the
!     corner nodes - only valid for quadratic elements
!
      implicit none
!
      character*81 objectset(4,*)
      character*8 lakon(*)
!
      integer i,ii,j,l,ielem,nodecor,nk,ne,
     &   nope,indexe,ipkon(*),konl(26),ifaceq(2,20),
     &   ifacet(2,10),ifacew(2,15),kon(*),nnodes(nk),nobject,
     &   nodedesiinv(nk),start
!
      real*8 dgdxglob(2,nk,nobject),xinterpol(nk,nobject)
!
      intent(in) ipkon,lakon,kon,nobject,ne,nk,objectset
!
      intent(inout) nodedesiinv,dgdxglob,xinterpol,nnodes
!
!     cornernodes next to the midnode for quadratic hex element
!
      data ifaceq /0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             1,2,
     &             2,3,
     &             3,4,
     &             1,4,
     &             5,6,
     &             6,7,
     &             7,8,
     &             5,8,
     &             1,5,
     &             2,6,
     &             3,7,
     &             4,8/
!
!     cornernodes next to the midnode for quadratic tet elements
!
      data ifacet /0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             1,2,
     &             2,3,
     &             1,3,
     &             1,4,
     &             2,4,
     &             3,4/
!
!     cornernodes next to the midnode for quadratic wedge elements
!
      data ifacew /0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             0,0,
     &             1,2,
     &             2,3,
     &             1,3,
     &             4,5,
     &             5,6,
     &             4,6,
     &             1,4,
     &             2,5,
     &             3,6/
!
!     Loop over all elements
!
      do ielem=1,ne
!
         if(ipkon(ielem).lt.0) cycle
!
!        Check if element is quadratic
!
         if(lakon(ielem)(4:5).eq.'10') then
            nope=10
            start=5
         elseif(lakon(ielem)(4:5).eq.'20') then
            nope=20
            start=9
         elseif (lakon(ielem)(4:5).eq.'15') then
            nope=15
            start=7
         else
            cycle
         endif

         indexe=ipkon(ielem)
!
         do l=1,nope
            konl(l)=kon(indexe+l)
         enddo
!
!        Loop over all midnodes of the element
!
         do i=start,nope
!
            if(nodedesiinv(konl(i)).le.0) cycle
!
!           Loop over the 2 neighbors of the midnode
!
            do j=1,2
               if(lakon(ielem)(4:5).eq.'10') then
                  nodecor=konl(ifacet(j,i))
               elseif (lakon(ielem)(4:5).eq.'20') then
                  nodecor=konl(ifaceq(j,i))
               elseif (lakon(ielem)(4:5).eq.'15') then
                  nodecor=konl(ifacew(j,i))
               endif
               if(nodedesiinv(nodecor).eq.1) then
                  nnodes(nodecor)=nnodes(nodecor)+1
                  do ii=1,nobject
                     xinterpol(nodecor,ii)=xinterpol(nodecor,ii)+
     &                    dgdxglob(1,konl(i),ii)
                  enddo
               endif
            enddo
            nodedesiinv(konl(i))=-1
         enddo
      enddo
!
      do i=1,nk
         if(nnodes(i).gt.0) then
            do j=1,nobject
               if(objectset(1,j)(1:9).eq.'THICKNESS') cycle
               if(objectset(1,j)(1:9).eq.'FIXGROWTH') cycle
               if(objectset(1,j)(1:12).eq.'FIXSHRINKAGE') cycle
               dgdxglob(1,i,j)=xinterpol(i,j)/nnodes(i)
            enddo
         endif
      enddo
!
!     Correction of nodedesiinv
!
      do i=1,nk
         if(nodedesiinv(i).eq.-1) then
            nodedesiinv(i)=1
         endif
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
*DECK DPCOEF
      SUBROUTINE DPCOEF (L, C, TC, A)
C***BEGIN PROLOGUE  DPCOEF
C***PURPOSE  Convert the DPOLFT coefficients to Taylor series form.
C***LIBRARY   SLATEC
C***CATEGORY  K1A1A2
C***TYPE      DOUBLE PRECISION (PCOEF-S, DPCOEF-D)
C***KEYWORDS  CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
C***AUTHOR  Shampine, L. F., (SNLA)
C           Davenport, S. M., (SNLA)
C***DESCRIPTION
C
C     Abstract
C
C     DPOLFT  computes the least squares polynomial fit of degree  L  as
C     a sum of orthogonal polynomials.  DPCOEF  changes this fit to its
C     Taylor expansion about any point  C , i.e. writes the polynomial
C     as a sum of powers of (X-C).  Taking  C=0.  gives the polynomial
C     in powers of X, but a suitable non-zero  C  often leads to
C     polynomials which are better scaled and more accurately evaluated.
C
C     The parameters for  DPCOEF  are
C
C     INPUT -- All TYPE REAL variables are DOUBLE PRECISION
C         L -      Indicates the degree of polynomial to be changed to
C                  its Taylor expansion.  To obtain the Taylor
C                  coefficients in reverse order, input  L  as the
C                  negative of the degree desired.  The absolute value
C                  of L  must be less than or equal to NDEG, the highest
C                  degree polynomial fitted by  DPOLFT .
C         C -      The point about which the Taylor expansion is to be
C                  made.
C         A -      Work and output array containing values from last
C                  call to  DPOLFT .
C
C     OUTPUT -- All TYPE REAL variables are DOUBLE PRECISION
C         TC -     Vector containing the first LL+1 Taylor coefficients
C                  where LL=ABS(L).  If  L.GT.0 , the coefficients are
C                  in the usual Taylor series order, i.e.
C                    P(X) = TC(1) + TC(2)*(X-C) + ... + TC(N+1)*(X-C)**N
C                  If L .LT. 0, the coefficients are in reverse order,
C                  i.e.
C                    P(X) = TC(1)*(X-C)**N + ... + TC(N)*(X-C) + TC(N+1)
C
C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C                 Curve fitting by polynomials in one variable, Report
C                 SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED  DP1VLU
C***REVISION HISTORY  (YYMMDD)
C   740601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DPCOEF
C
      INTEGER I,L,LL,LLP1,LLP2,NEW,NR
      DOUBLE PRECISION A(*),C,FAC,SAVE,TC(*)
C***FIRST EXECUTABLE STATEMENT  DPCOEF
      LL = ABS(L)
      LLP1 = LL + 1
      CALL DP1VLU (LL,LL,C,TC(1),TC(2),A)
      IF (LL .LT. 2) GO TO 2
      FAC = 1.0D0
      DO 1 I = 3,LLP1
        FAC = FAC*(I-1)
 1      TC(I) = TC(I)/FAC
 2    IF (L .GE. 0) GO TO 4
      NR = LLP1/2
      LLP2 = LL + 2
      DO 3 I = 1,NR
        SAVE = TC(I)
        NEW = LLP2 - I
        TC(I) = TC(NEW)
 3      TC(NEW) = SAVE
 4    RETURN
      END


text:
      SUBROUTINE G1CLS3(IFUNCO,LENGTH,IDATA,RDATA)
C
C          ------------------------------------------------
C          ROUTINE NO. (1014)   VERSION (A8.7)    07:NOV:96
C          ------------------------------------------------
C
C          THIS PROVIDES PARAMETER-SETTING OPERATIONS.
C
C
C          THE ARGUMENTS ARE AS FOLLOWS:
C
C          <IFUNCO> GIVES THE FUNCTION CODE:
C                   =  1, OPEN DEVICE CHANNEL             (N),
C                   =  2, CLOSE DEVICE CHANNEL            (N),
C                   =  3, ERASE PICTURE WITH PROMPT/WAIT,
C                   =  4, SET LINE THICKNESS              (N),
C                   =  5, -NO LONGER USED-
C                   =  6, SET LINE PATTERN                (N1...N4),
C                   =  7, SET LINE HARDWARE MODE          (N),
C                   =  8, SET LINE BLINK                  (N),
C                   =  9, SET LINE COLOUR NUMBER          (N),
C                   = 10, EMPTY DEVICE BUFFERS,
C                   = 11, ERASE PICTURE WITHOUT PROMPT/WAIT,
C                   = 12, SET DEVICE RESOLUTION           (N).
C          <LENGTH> GIVES THE DATA LENGTH (IN BYTES),
C          [IDATA]  IS THE INTEGER DATA ARRAY, AND
C          [RDATA]  IS THE ALTERNATIVE REAL DATA ARRAY.
C
C
      REAL    RDATA(*)
      INTEGER IDATA(*)
      LOGICAL OPEN,LHRDW1
C
      COMMON /T1CPOS/ XCHAR,YCHAR
      COMMON /T1DENA/ OPEN
      COMMON /T1DEVS/ KCHAN1(5),IRESL1
      COMMON /T1DRES/ DRESX,DRESY
      COMMON /T1HRDL/ LHRDW1
      COMMON /T1KDEF/ KOLLID,KOLBAD,MAXCLS
      COMMON /T1LATT/ KOLIN1,ITHIK1
      COMMON /T1LINT/ ELWID,NFILLI,THIKST
      COMMON /T1LPAT/ MARKA1,MISSA1,MARKB1,MISSB1
      COMMON /T1MAPP/ X1MAP1,X2MAP1,Y1MAP1,Y2MAP1
      COMMON /T1PNUM/ NDPICS
      COMMON /T1PPBT/ XPLOTB,YPLOTB
      COMMON /T1PPOS/ XPLOT1,YPLOT1
C
C
      IF (IFUNCO.LT.1.OR.IFUNCO.GT.14) RETURN
      GO TO (1,2,3,4,5,6,7,8,9,10,11,13), IFUNCO
C
C          THIS OPENS THE DEVICE CHANNEL. (WITHOUT THIS
C          BEING DONE, NO GRAPHICAL OUTPUT IS POSSIBLE).
C
    1 IVALUE= IDATA(1)
      IF (IVALUE.LE.0.OR.IVALUE.GT.5) RETURN
C
      KCHAN1(IVALUE)= 1
      CALL G1FILB(0,0,-86,IVALUE)
      OPEN= .TRUE.
      RETURN
C
C          AND THIS CLOSES THE DEVICE CHANNEL.
C
    2 IVALUE= IDATA(1)
      IF (IVALUE.LE.0.OR.IVALUE.GT.5) RETURN
C
      KCHAN1(IVALUE)= 0
      CALL G1FILB(0,0,-87,IVALUE)
      OPEN= .FALSE.
C
      DO 100 ITEST= 1,5
        IF (KCHAN1(ITEST).EQ.1) OPEN= .TRUE.
  100 CONTINUE
      RETURN
C
C          THIS ERASES THE PICTURE WITH A PROMPT/WAIT.
C
    3 IF (.NOT.OPEN) RETURN
C
      KODERA= -81
      GO TO 12
C
C          THIS SETS THE LINE THICKNESS.
C
    4 ITHIK1= IDATA(1)
      NFILLI= (ITHIK1/ELWID)+0.5
      IF (ITHIK1.EQ.1.OR.NFILLI.LT.1) NFILLI= 1
C
      THIKST= (0.0005*ITHIK1)/NFILLI
      RETURN
C
C          THIS SETS THE LINE INTENSITY.
C          (IT IS NO LONGER IMPLEMENTED).
C
    5 RETURN
C
C          THIS SETS THE CURRENT LINE PATTERN, AND IN
C          THE CASE OF HARDWARE BROKEN LINES SETS THE
C          APPROPRIATE PATTERN VIA SUBROUTINE <G1HRDW>.
C
    6 MARKA1= IDATA(1)
      MISSA1= IDATA(2)
      MARKB1= IDATA(3)
      MISSB1= IDATA(4)
      RESMIN= 1000.0/AMIN1(DRESX,DRESY)
      IF (FLOAT(MARKA1).LT.RESMIN) MARKA1= 0
      IF (FLOAT(MISSA1).LT.RESMIN) MISSA1= 0
      IF (FLOAT(MARKB1).LT.RESMIN) MARKB1= 0
      IF (FLOAT(MISSB1).LT.RESMIN) MISSB1= 0
      IF (LHRDW1) CALL G1HRDW(1)
      RETURN
C
C          THIS SETS LINE HARDWARE MODE.
C
    7 LHRDW1= .FALSE.
      IF (IDATA(1).EQ.1) LHRDW1= .TRUE.
C
      CALL G1HRDW(1)
      RETURN
C
C          THIS SETS DISPLAY BLINKING.
C          (IT IS NOT IMPLEMENTED HERE).
C
    8 RETURN
C
C          THIS SETS THE LINE COLOUR NUMBER.
C
    9 KOLIN1= IDATA(1)
      IF (KOLIN1.EQ.0) KOLIN1= KOLLID
C
      CALL G1HRDW(4)
      RETURN
C
C          THIS CAUSES ALL THE INTERNAL BUFFERS TO
C          BE EMPTIED IMMEDIATELY RATHER THAN AFTER
C          THEY BECOME FULL (AS IS THE NORMAL CASE).
C
   10 CALL G1FILB(0,0,-84,0)
      RETURN
C
C          THIS ERASES THE DEVICE WITHOUT A PROMPT/WAIT.
C
   11 IF (.NOT.OPEN) RETURN
C
      KODERA= -82
   12 XPLOT1= X1MAP1
      YPLOT1= Y1MAP1
      CALL G1LINE(0.0,0.0,.FALSE.)
      XCHAR= XPLOTB
      YCHAR= YPLOTB
      NDPICS= NDPICS+1
      CALL G1FILB(0,0,KODERA,0)
      RETURN
C
C          THIS SETS THE DEVICE RESOLUTION.
C
   13 IRESL1= IDATA(1)
      CALL G1HRDW(9)
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
SUBROUTINE CALBED

  !**********************************************************************!
  ! ***  SUBROUTINE CALSND CALCULATES NONCOHESIVE SEDIMENT SETTLING,
  ! ***  DEPOSITION AND RESUSPENSION AND IS CALLED FOR SSEDTOX
  ! ***  NOT USED FOR SEDZLJ
  !**********************************************************************!
  !
  !----------------------------------------------------------------------!
  ! CHANGE RECORD
  ! DATE MODIFIED     BY               DESCRIPTION
  !----------------------------------------------------------------------!
  ! 2011-03           Paul M. Craig      Rewritten to F90 and added OMP
  !----------------------------------------------------------------------!

  USE GLOBAL

  IMPLICIT NONE

  INTEGER :: K,L,IFLAG,NS,KK,NX,KBTM1,ND,LF,LL
  REAL :: TMPVAL,WDENKGM3,WDENGMM3,TMPVALK
  REAL :: TMPVALKP,BETTMP,VOIDCON1,HBEDTMP,TMPVALO
  REAL :: TMPVALN,TMPEXP,TMPTOP,TMPBOT,FSTRSE,FDSTRSE
  REAL :: FHYDCN,DSTRESET,POROTMP

  REAL,SAVE,ALLOCATABLE,DIMENSION(:,:) :: SNDHYDCN

  IF(  .NOT. ALLOCATED(SNDHYDCN) )THEN
    ALLOCATE(SNDHYDCN(LCM,KBM))
    SNDHYDCN = 0.0
  ENDIF

  IF( ISTRAN(6) >= 1 .OR. ISTRAN(7) >= 1 )THEN
    HBEDMIN = 1.E-4
    IF( ISTRAN(7) >= 1 )THEN
      HBEDMIN = MAX(HBEDMIN,SNDDMX)
    END IF

    ! *********************************************************************
    ! *** CONSTANT POROSITY BED
    IF( IBMECH == 0 )THEN

      ! *** UPDATE TOP LAYER THICKNESS TO MAINTAIN CONSTANT POROSITY
      VOIDCON1 = BEDPORC/(1.-BEDPORC)

      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ND,K,L,LF,LL)              &
      !$OMP             PRIVATE(HBEDTMP,TMPVALO,TMPVALN)
      DO ND=1,NDM
        LF = 2 + (ND-1)*LDM
        LL = MIN(LF+LDM-1,LA)

        DO L=LF,LL
          K = KBT(L)
          HBEDTMP = (1. + VOIDCON1)*HBED(L,K)/(1. + VDRBED(L,K))
          TMPVALO = VDRBED(L,K)*HBED(L,K)/(1. + VDRBED(L,K))
          TMPVALN = VOIDCON1*HBEDTMP/(1. + VOIDCON1)
          QWBDTOP(L) = DELTI*(TMPVALO-TMPVALN)
          HBED(L,K) = HBEDTMP
          QWTRBED(L,K) = QWBDTOP(L) + QGW(L)/DXYP(L)
          VDRBED(L,K) = VOIDCON1     ! *** PMC 2010-12-06
        ENDDO
        DO K=0,KBT(L)-1
          DO L=LF,LL
            QWTRBED(L,K) = QGW(L)/DXYP(L)
          ENDDO
        END DO

      ENDDO   ! *** END OF DOMAIN LOOP
      !$OMP END PARALLEL DO

    ENDIF

    ! *********************************************************************
    ! *** SIMPLE CONSOLIDATING BED
    IF( IBMECH == 1 )THEN

      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ND,K,L,LF,LL)              &
      !$OMP             PRIVATE(TMPEXP,TMPTOP,TMPBOT,TMPVAL)
      DO ND=1,NDM
        LF = 2+(ND-1)*LDM
        LL = MIN(LF+LDM-1,LA)

        ! ***  IF SEDVRDT>0.0001 CONSOLIDATE TO SEDVRM (THE MINIMUM VOID RATIO)
        IF( SEDVRDT > 0.0001 )THEN
          TMPEXP = EXP(-DTSED/SEDVRDT)
          DO K=1,KB
            DO L=LF,LL
              IF( K <= KBT(L) )THEN
                VDRBED1(L,K) = VDRBED(L,K)
                HBED1(L,K) = HBED(L,K)
                VDRBED(L,K) = SEDVDRM + (VDRBED1(L,K)-SEDVDRM)*TMPEXP
                TMPTOP = 1. + VDRBED(L,K)
                TMPBOT = 1. + VDRBED1(L,K)
                HBED(L,K) = TMPTOP*HBED1(L,K)/TMPBOT
              ENDIF
            ENDDO
          ENDDO
        ENDIF

        ! ***  IF SEDVRDT>0.0001 CONSOLIDATE TO SEDVRM INSTANTANEOUSLY
        IF( SEDVRDT >= 0.0 .AND. SEDVRDT <= 0.0001 )THEN
          TMPEXP = 0.0
          DO K=1,KB
            DO L=LF,LL
              IF( K <= KBT(L) )THEN
                VDRBED1(L,K) = VDRBED(L,K)
                HBED1(L,K) = HBED(L,K)
                VDRBED(L,K) = SEDVDRM
                TMPTOP = 1. + VDRBED(L,K)
                TMPBOT = 1. + VDRBED1(L,K)
                HBED(L,K) = TMPTOP*HBED1(L,K)/TMPBOT
              ENDIF
            ENDDO
          ENDDO
        ENDIF

        ! ***  IF SEDVRDT<0.0 MAINTAIN INITIAL VOID RATIO (SAVED IN VDRBED2)
        IF( SEDVRDT < 0.0 )THEN
          TMPEXP = 1.0
          DO L=LF,LL
            K = KBT(L)
            VDRBED1(L,K) = VDRBED(L,K)
            HBED1(L,K) = HBED(L,K)
            VDRBED(L,K) = VDRBED2(L,K)
            TMPTOP = 1. + VDRBED(L,K)
            TMPBOT = 1. + VDRBED1(L,K)
            HBED(L,K) = TMPTOP*HBED1(L,K)/TMPBOT
          ENDDO
          DO L=LF,LL
            K = KBT(L)-1
            IF( K > 0 )THEN
              VDRBED1(L,K) = VDRBED(L,K)
              HBED1(L,K) = HBED(L,K)
              VDRBED(L,K) = VDRBED2(L,K)
              TMPTOP = 1. + VDRBED(L,K)
              TMPBOT = 1. + VDRBED1(L,K)
              HBED(L,K) = TMPTOP*HBED1(L,K)/TMPBOT
            ENDIF
          ENDDO
        ENDIF

        ! *** UPDATE POROSITY
        DO K=1,KB
          DO L=LF,LL
            IF( K <= KBT(L) )THEN
              PORBED(L,K)  = VDRBED(L,K) /(1. + VDRBED(L,K))
              PORBED1(L,K) = VDRBED1(L,K)/(1. + VDRBED1(L,K))
            ENDIF
          ENDDO
        ENDDO

        ! *** UPDATE PORE WATER FLOWS
        DO L=LF,LL
          QWTRBED(L,0) = QGW(L)/DXYP(L)
        ENDDO
        DO K=1,KB
          DO L=LF,LL
            IF( K <= KBT(L) )THEN
              TMPVAL = HBED(L,K)/(1. + VDRBED(L,K))
              QWTRBED(L,K) = QWTRBED(L,K-1)  -DELTI*TMPVAL*(VDRBED(L,K)-VDRBED1(L,K))
            ENDIF
          ENDDO
        ENDDO

      ENDDO   ! *** END OF DOMAIN LOOP
      !$OMP END PARALLEL DO

    ENDIF

    ! *********************************************************************
    ! *** CONSOLIDATE TO VOID RATIO/POROSITY DEFINED AS FUNCTION OF
    ! *** COHESIVE SEDIMENT FRACTION, FRACCOH(L,K)
    IF( IBMECH == 99 )THEN

      !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ND,K,L,LF,LL)              &
      !$OMP             PRIVATE(TMPEXP,TMPTOP,TMPBOT,TMPVAL,POROTMP)
      DO ND=1,NDM
        LF = 2+(ND-1)*LDM
        LL = MIN(LF+LDM-1,LA)

        DO K=1,KB
          DO L=LF,LL
            IF( K <= KBT(L) )THEN
              VDRBED1(L,K) = VDRBED(L,K)
              HBED1(L,K) = HBED(L,K)
              POROTMP = BMECH1*(FRACCOH(L,K)**BMECH2) + BMECH3  ! *** PMC - THIS APPROACH DISCONNECTS THE SEDIMENT MASS FROM THE VOID RATIO
              VDRBED(L,K) = POROTMP/(1.-POROTMP)
              TMPTOP = 1. + VDRBED(L,K)
              TMPBOT = 1. + VDRBED1(L,K)
              HBED(L,K) = TMPTOP*HBED1(L,K)/TMPBOT
            ENDIF
          ENDDO
        ENDDO

        ! *** UPDATE POROSITY
        DO K=1,KB
          DO L=LF,LL
            IF( K <= KBT(L) )THEN
              PORBED(L,K) = VDRBED(L,K)/(1. + VDRBED(L,K))
              PORBED1(L,K) = VDRBED1(L,K)/(1. + VDRBED1(L,K))
            ENDIF
          ENDDO
        ENDDO

        ! *** UPDATE PORE WATER FLOWS
        DO L=LF,LL
          QWTRBED(L,0) = QGW(L)/DXYP(L)
        ENDDO

        DO K=1,KB
          DO L=LF,LL
            IF( K <= KBT(L) )THEN
              TMPVAL = HBED(L,K)/(1. + VDRBED(L,K))
              QWTRBED(L,K) = QWTRBED(L,K-1)-DELTI*TMPVAL*(VDRBED(L,K)-VDRBED1(L,K))
            ENDIF
          ENDDO
        ENDDO

      ENDDO   ! *** END OF DOMAIN LOOP
      !$OMP END PARALLEL DO

    ENDIF     ! *** END OF IBMECH = 99

    ! *********************************************************************
    ! *** SET FLAG FOR FINITE STRAIN CONSOLIDATION
    IFLAG = 0
    IF( ISTRAN(6) >= 1 )IFLAG = IFLAG + 1
    IF( ISTRAN(7) >= 1 )IFLAG = IFLAG + 1

    ! *** FINITE STRAIN CONSOLIDATING HOMOGENEOUS BED
    IF( IBMECH >= 2 .AND. IBMECH < 98 .AND. IFLAG == 1 )THEN
      WDENKGM3 = 1.E3
      WDENGMM3 = 1.E6

      ! ***  SET PHYSICAL VERTICAL COORDINATES OF THE BED

      ! ***  ZBEDC = VERTICAL COORDINATE OF THE CENTER OF BED LAYER
      ! ***  ZBEDG = VERTICAL COORDINATE AT TOP OF BED LAYER
      DO L=2,LA
        ZBEDG(L,0) = ZELBEDA(L)
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDG(L,K) = ZBEDG(L,K-1) + HBED(L,K)
          ELSE
            ZBEDG(L,K) = HBED(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        ZBEDGT(L) = ZBEDG(L,KBT(L))
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDC(L,K) = 0.5*(ZBEDG(L,K) + ZBEDG(L,K-1))
          ELSE
            ZBEDC(L,K) = 0.5*ZBEDG(L,K)
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE TRANSFORMED THICKNESS OF BED LAYERS
      ! ***  DZBTR = TRANSFORMED THICKNESS OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            DZBTR(L,K) = HBED(L,K)/(1. + VDRBED(L,K))
            DZBTR1(L,K) = HBED1(L,K)/(1. + VDRBED1(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE WATER SPECIFIC WEIGHT NORMALIZED
      ! ***    EFFECTIVE STRESS USING FSTRSE
      ! ***  CALCULATE DERIVATIVE OF EFFECTIVE STRESS
      ! ***    WITH RESPECT TO VOID RATIO, DSTRSE USING
      ! ***    FUNCTION FDSTRSE
      ! ***  CALCULATE HYDRAULIC CONDUCTIVITY DIVIDED BY (1+VOID),
      ! ***    HYDCN =USING FUNCTION FHYDCN
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            STRSE(L,K)=  FSTRSE(VDRBED(L,K),BMECH1,BMECH2,BMECH3)
            DSTRSE(L,K)=  FDSTRSE(VDRBED(L,K),BMECH1,BMECH2,BMECH3)
            HYDCN(L,K)=  FHYDCN(VDRBED(L,K),BMECH4,BMECH5,BMECH6,IBMECHK)
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          KBTM1 = KBT(L)-1
          IF( K <= KBTM1 )THEN
            TMPVAL = ( DZBTR(L,K)/HYDCN(L,K) )  +( DZBTR(L,K+1)/HYDCN(L,K+1) )
            COEFK(L,K) = (DZBTR(L,K) + DZBTR(L,K+1))/TMPVAL
            DSTRESET = (DZBTR(L,K)*DSTRSE(L,K+1)  +DZBTR(L,K+1)*DSTRSE(L,K))  /(DZBTR(L,K) + DZBTR(L,K+1))
            COEFSK(L,K) = DSTRESET*COEFK(L,K)
          ENDIF
          IF( K == KBT(L) )THEN
            COEFK(L,K) = HYDCN(L,K)
            COEFSK(L,K) = DSTRSE(L,K)*HYDCN(L,K)
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE PRESSURE COMPONENTS
      DO K=1,KB
        DO L=2,LA
          SGSM1(L,K) = 0.
        ENDDO
      ENDDO
      IF( ISTRAN(6) > 0 )THEN
        DO NS=1,NSED
          DO K=1,KB
            DO L=2,LA
              IF( K <= KBT(L) )THEN
                SGSM1(L,K) = SGSM1(L,K) + SSG(NS)*VFRBED(L,K,NS)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(7) > 0 )THEN
        DO NX=1,NSND
          NS = NSED + NX
          DO K=1,KB
            DO L=2,LA
              IF( K <= KBT(L) )THEN
                SGSM1(L,K) = SGSM1(L,K) + SSG(NS)*VFRBED(L,K,NS)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDIF
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            SGSM1(L,K) = SGSM1(L,K)-1.
          ENDIF
        ENDDO
      ENDDO

      ! ***  NEW IMPLICIT CONSOLIDATION SOLUTION BEGINS HERE
      DO K=1,KB
        DO L=2,LA
          IF( K<KBT(L) )THEN
            ACOEF(L,K) = 2.*COEFSK(L,K)/(DZBTR(L,K) + DZBTR(L,K+1))
            TMPVALK = DZBTR(L,K)*SGSM1(L,K)
            TMPVALKP = DZBTR(L,K+1)*SGSM1(L,K+1)
            QCOEF(L,K) = (TMPVALK + TMPVALKP)*COEFK(L,K)/  (DZBTR(L,K) + DZBTR(L,K+1))
          ELSE
            ACOEF(L,K) = 0.0
            QCOEF(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        K = KBT(L)
        ACOEF(L,K) = 2.*COEFSK(L,K)/DZBTR(L,K)
        QCOEF(L,K) = SGSM1(L,K)*COEFK(L,K)
        QCOEF(L,0) = QGW(L)*DXYIP(L)
        QWTRBED(L,0) = QGW(L)*DXYIP(L)
      ENDDO
      DO L=2,LA
        ALOW(L,1) = 0.
        CUPP(L,KBT(L)) = 0.
        DO K=1,KBT(L)-1
          CUPP(L,K) = -DTSED*ACOEF(L,K)/DZBTR(L,K)
        ENDDO
        DO K=2,KBT(L)
          ALOW(L,K) = -DTSED*ACOEF(L,K-1)/DZBTR(L,K)
        ENDDO
        DO K=1,KBT(L)
          BMNN(L,K) = 1.0-ALOW(L,K)-CUPP(L,K)
        ENDDO
        K = KBT(L)
        BMNN(L,K) = BMNN(L,K) + DTSED*ACOEF(L,K)/DZBTR(L,K)
        DO K=1,KBT(L)
          RRHS(L,K) = VDRBED(L,K)  +DTSED*(QCOEF(L,K-1)-QCOEF(L,K))/DZBTR(L,K)
          VDRBED1(L,K) = VDRBED(L,K)
          HBED1(L,K) = HBED(L,K)
        ENDDO
      ENDDO

      IF( IBMECH == 2 )THEN
        DO L=2,LA
          K = KBT(L)
          RRHS(L,K) = RRHS(L,K) + DTSED*ACOEF(L,K)*VDRDEPO(1)/DZBTR(L,K)
        ENDDO
      ENDIF
      IF( IBMECH == 3 )THEN
        DO L=2,LA
          K = KBT(L)
          RRHS(L,K) = RRHS(L,K) + DTSED*ACOEF(L,K)*(VDRBED(L,K)  +(STRSE(L,K)/DSTRSE(L,K)))/DZBTR(L,K)
        ENDDO
      ENDIF
      DO L=2,LA
        BETTMP = BMNN(L,1)
        TOXTMP(L,1) = RRHS(L,1)/BETTMP
        DO KK=2,KBT(L)
          GAMTMP(L,KK) = CUPP(L,KK-1)/BETTMP
          BETTMP = BMNN(L,KK)-ALOW(L,KK)*GAMTMP(L,KK)
          TOXTMP(L,KK) = (RRHS(L,KK)-ALOW(L,KK)*TOXTMP(L,KK-1))/  BETTMP
        ENDDO
        DO KK=KBT(L)-1,1,-1
          TOXTMP(L,KK) = TOXTMP(L,KK)-GAMTMP(L,KK+1)*TOXTMP(L,KK+1)
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K<KBT(L) )THEN
            QWTRBED(L,K) = -ACOEF(L,K)*(TOXTMP(L,K+1)-TOXTMP(L,K))  +QCOEF(L,K)
          ELSE
            QWTRBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      IF( IBMECH == 2 )THEN
        DO L=2,LA
          K = KBT(L)
          QWTRBED(L,K) = -ACOEF(L,K)*(SEDVDRD-TOXTMP(L,K)) + QCOEF(L,K)
        ENDDO
      ELSE
        DO L=2,LA
          K = KBT(L)
          QWTRBED(L,K) = -ACOEF(L,K)*(VDRBED(L,K)  +(STRSE(L,K)/DSTRSE(L,K))-TOXTMP(L,K)) + QCOEF(L,K)
        ENDDO
      END IF

      ! ***  CALCULATE VOID RATIOS
      ! ***  VDRBED =  VOID RATIO OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            VDRBED(L,K) = VDRBED1(L,K)  -DTSED*(QWTRBED(L,K)-QWTRBED(L,K-1))/DZBTR1(L,K)
          ELSE
            VDRBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO

      ! ***  UPDATE LAYER THICKNESS
      ! ***  HBED = BED LAYER THICKNESS
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            HBED(L,K) = HBED1(L,K)*(1. + VDRBED(L,K))/(1. + VDRBED1(L,K))
          ELSE
            HBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        QWTRBED(L,0) = QGW(L)/DXYP(L)
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            QWTRBED(L,K) = QWTRBED(L,K-1)  -DELTI*(HBED(L,K)-HBED1(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! ***  ZBEDC = VERTICAL COORDINATE OF THE CENTER OF BED LAYER
      ! ***  ZBEDG = VERTICAL COORDINATE AT TOP OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDG(L,K) = ZBEDG(L,K-1) + HBED(L,K)
          ELSE
            ZBEDG(L,K) = HBED(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        ZBEDGT(L) = ZBEDG(L,KBT(L))
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDC(L,K) = 0.5*(ZBEDG(L,K) + ZBEDG(L,K-1))
          ELSE
            ZBEDC(L,K) = 0.5*ZBEDG(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            DZBTR(L,K) = HBED(L,K)/(1. + VDRBED(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! *** UPDATE POROSITY
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            PORBED(L,K) = VDRBED(L,K)/(1. + VDRBED(L,K))
            PORBED1(L,K) = VDRBED1(L,K)/(1. + VDRBED1(L,K))
          ELSE
            PORBED(L,K) = 0.0
            PORBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO

      ! *** ADD OR REMOVE LAYERS
      !CALL CALBLAY

    ENDIF

    ! *********************************************************************
    ! *** FINITE STRAIN CONSOLIDATING NON-HOMOGENEOUS BED
    IF( IBMECH >= 2 .AND. IBMECH<98 .AND. IFLAG == 2 )THEN
      WDENKGM3 = 1.E3
      WDENGMM3 = 1.E6

      ! ***  SET PHYSICAL VERTICAL COORDINATES OF THE BED

      ! ***  ZBEDC = VERTICAL COORDINATE OF THE CENTER OF BED LAYER
      ! ***  ZBEDG = VERTICAL COORDINATE AT TOP OF BED LAYER
      DO L=2,LA
        ZBEDG(L,0) = ZELBEDA(L)
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDG(L,K) = ZBEDG(L,K-1) + HBED(L,K)
          ELSE
            ZBEDG(L,K) = HBED(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        ZBEDGT(L) = ZBEDG(L,KBT(L))
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDC(L,K) = 0.5*(ZBEDG(L,K) + ZBEDG(L,K-1))
          ELSE
            ZBEDC(L,K) = 0.5*ZBEDG(L,K)
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE TRANSFORMED THICKNESS OF BED LAYERS
      ! ***  DZBTR = TRANSFORMED THICKNESS OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            DZBTR(L,K) = HBED(L,K)/(1. + VDRBED(L,K))
            DZBTR1(L,K) = HBED1(L,K)/(1. + VDRBED1(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE WATER SPECIFIC WEIGHT NORMALIZED
      ! ***    EFFECTIVE STRESS USING FSTRSE
      ! ***  CALCULATE DERIVATIVE OF EFFECTIVE STRESS
      ! ***    WITH RESPECT TO VOID RATIO, DSTRSE USING
      ! ***    FUNCTION FDSTRSE
      ! ***  CALCULATE HYDRAULIC CONDUCTIVITY DIVIED BY (1+VOID),
      ! ***    HYDCN =USING FUNCTION FHYDCN
      ! ***  NONCOHESIVE HYDRAULIC CONDUCTIVITY BASED ON R. R. RUMER, CHAP 3,
      ! ***  EQ 13 AND 14, IN 'FLOW THROUGH POROUS MEDIA' ED. R. J. M. DE WIEST
      ! ***  ACADEMIC PRESS, 1969
      ! ***  KH = 2854*(n**2)*(d**2)/(1-n) = 2854*(e**2)*(d**2)/(1+e)
      ! ***  FOR KH IN M/S AND D IN METERS
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            STRSE(L,K)=  FSTRSE(VDRBEDSED(L,K),BMECH1,BMECH2,BMECH3)
            DSTRSE(L,K)=  FDSTRSE(VDRBEDSED(L,K),BMECH1,BMECH2,BMECH3)
            HYDCN(L,K)=  FHYDCN(VDRBEDSED(L,K),BMECH4,BMECH5,BMECH6,IBMECHK)
            TMPVAL = VDRBEDSND(L,K)/(1.+VDRBEDSND(L,K))
            SNDHYDCN(L,K) = 2854.0*TMPVAL*TMPVAL*(SEDDIA50(L,K)**2)
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            COEFK(L,K) = (FRACCOH(L,K)/HYDCN(L,K))  +(FRACNON(L,K)/SNDHYDCN(L,K))
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            HYDCN(L,K) = 1./COEFK(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          COEFK(L,K) = 0.
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          KBTM1 = KBT(L)-1
          IF( K <= KBTM1 )THEN
            TMPVAL = ( DZBTR(L,K)/HYDCN(L,K) )  +( DZBTR(L,K+1)/HYDCN(L,K+1) )
            COEFK(L,K) = (DZBTR(L,K) + DZBTR(L,K+1))/TMPVAL
            DSTRESET = (DZBTR(L,K)*DSTRSE(L,K+1)  +DZBTR(L,K+1)*DSTRSE(L,K))  /(DZBTR(L,K) + DZBTR(L,K+1))
            COEFSK(L,K) = DSTRESET*COEFK(L,K)
          ENDIF
          IF( K == KBT(L) )THEN
            COEFK(L,K) = HYDCN(L,K)
            COEFSK(L,K) = DSTRSE(L,K)*HYDCN(L,K)
          ENDIF
        ENDDO
      ENDDO

      ! ***  CALCULATE PRESSURE COMPONENTS
      DO K=1,KB
        DO L=2,LA
          SGSM1(L,K) = 0.
        ENDDO
      ENDDO
      IF( ISTRAN(6) > 0 )THEN
        DO NS=1,NSED
          DO K=1,KB
            DO L=2,LA
              IF( K <= KBT(L) )THEN
                SGSM1(L,K) = SGSM1(L,K) + SSG(NS)*VFRBED(L,K,NS)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDIF
      IF( ISTRAN(7) > 0 )THEN
        DO NX=1,NSND
          NS = NSED + NX
          DO K=1,KB
            DO L=2,LA
              IF( K <= KBT(L) )THEN
                SGSM1(L,K) = SGSM1(L,K) + SSG(NS)*VFRBED(L,K,NS)
              ENDIF
            ENDDO
          ENDDO
        ENDDO
      ENDIF
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            SGSM1(L,K) = SGSM1(L,K)-1.
          ENDIF
        ENDDO
      ENDDO

      ! ***  NEW IMPLICIT CONSOLIDATION SOLUTION BEGINS HERE
      DO K=1,KB
        DO L=2,LA
          IF( K<KBT(L) )THEN
            ACOEF(L,K) = 2.*COEFSK(L,K)/(DZBTR(L,K) + DZBTR(L,K+1))
            TMPVALK = DZBTR(L,K)*SGSM1(L,K)
            TMPVALKP = DZBTR(L,K+1)*SGSM1(L,K+1)
            QCOEF(L,K) = (TMPVALK + TMPVALKP)*COEFK(L,K)/  (DZBTR(L,K) + DZBTR(L,K+1))
          ELSE
            ACOEF(L,K) = 0.0
            QCOEF(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        K = KBT(L)
        ACOEF(L,K) = 2.*COEFSK(L,K)/DZBTR(L,K)
        QCOEF(L,K) = SGSM1(L,K)*COEFK(L,K)
        QCOEF(L,0) = QGW(L)*DXYIP(L)
        QWTRBED(L,0) = QGW(L)*DXYIP(L)
      ENDDO
      DO L=2,LA
        ALOW(L,1) = 0.
        CUPP(L,KBT(L)) = 0.
        DO K=1,KBT(L)-1
          CUPP(L,K) = -DTSED*ACOEF(L,K)/DZBTR(L,K)
        ENDDO
        DO K=2,KBT(L)
          ALOW(L,K) = -DTSED*ACOEF(L,K-1)/DZBTR(L,K)
        ENDDO
        DO K=1,KBT(L)
          BMNN(L,K) = FRACCOH(L,K)-ALOW(L,K)-CUPP(L,K)
        ENDDO
        K = KBT(L)
        BMNN(L,K) = BMNN(L,K) + DTSED*ACOEF(L,K)/DZBTR(L,K)
        DO K=1,KBT(L)
          RRHS(L,K) = FRACCOH(L,K)*VDRBEDSED(L,K)  +DTSED*(QCOEF(L,K-1)-QCOEF(L,K))/DZBTR(L,K)
          VDRBED1(L,K) = VDRBED(L,K)
          HBED1(L,K) = HBED(L,K)
        ENDDO
      ENDDO
      IF( IBMECH == 2 )THEN
        DO L=2,LA
          K = KBT(L)
          RRHS(L,K) = RRHS(L,K) + DTSED*ACOEF(L,K)*VDRDEPO(1)/DZBTR(L,K)
        ENDDO
      ENDIF
      IF( IBMECH == 3 )THEN
        DO L=2,LA
          K = KBT(L)
          RRHS(L,K) = RRHS(L,K)  +DTSED*ACOEF(L,K)*(VDRBEDSED(L,K)  +(STRSE(L,K)/DSTRSE(L,K)))/DZBTR(L,K)
        ENDDO
      ENDIF
      DO L=2,LA
        BETTMP = BMNN(L,1)
        TOXTMP(L,1) = RRHS(L,1)/BETTMP
        DO KK=2,KBT(L)
          GAMTMP(L,KK) = CUPP(L,KK-1)/BETTMP
          BETTMP = BMNN(L,KK)-ALOW(L,KK)*GAMTMP(L,KK)
          TOXTMP(L,KK) = (RRHS(L,KK)-ALOW(L,KK)*TOXTMP(L,KK-1))/  BETTMP
        ENDDO
        DO KK=KBT(L)-1,1,-1
          TOXTMP(L,KK) = TOXTMP(L,KK)-GAMTMP(L,KK+1)*TOXTMP(L,KK+1)
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K<KBT(L) )THEN
            QWTRBED(L,K) = -ACOEF(L,K)*(TOXTMP(L,K+1)-TOXTMP(L,K))  +QCOEF(L,K)
          ELSE
            QWTRBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      IF( IBMECH == 2 )THEN
        DO L=2,LA
          K = KBT(L)
          QWTRBED(L,K) = -ACOEF(L,K)*(VDRDEPO(1) - TOXTMP(L,K)) + QCOEF(L,K)
        ENDDO
      ELSE
        DO L=2,LA
          K = KBT(L)
          QWTRBED(L,K) = -ACOEF(L,K)*(VDRBEDSED(L,K) + (STRSE(L,K)/DSTRSE(L,K))-TOXTMP(L,K)) + QCOEF(L,K)
        ENDDO
      END IF

      ! ***  CALCULATE VOID RATIOS
      ! ***  VDRBED =  VOID RATIO OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            IF( FRACCOH(L,K)>0.0 )THEN
            VDRBEDSED(L,K) = VDRBEDSED(L,K) - DTSED*(QWTRBED(L,K)-QWTRBED(L,K-1))/(FRACCOH(L,K)*DZBTR1(L,K))
              VDRBEDSED(L,K) = MAX(VDRBEDSED(L,K),VDRBEDSND(L,K))
            ELSE
              VDRBEDSED(L,K) = 0.0
            ENDIF
          ELSE
            VDRBEDSED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO

      ! ***  UPDATE VOID RATIO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            VDRBED(L,K) = FRACCOH(L,K)*VDRBEDSED(L,K) + FRACNON(L,K)*VDRBEDSND(L,K)
          ELSE
            VDRBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO

      ! ***  UPDATE LAYER THICKNESS
      ! ***  HBED = BED LAYER THICKNESS
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            HBED(L,K) = HBED1(L,K)*(1. + VDRBED(L,K))/(1. + VDRBED1(L,K))
          ELSE
            HBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        QWTRBED(L,0) = QGW(L)/DXYP(L)
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            QWTRBED(L,K) = QWTRBED(L,K-1)  -DELTI*(HBED(L,K)-HBED1(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! ***  ZBEDC = VERTICAL COORDINATE OF THE CENTER OF BED LAYER
      ! ***  ZBEDG = VERTICAL COORDINATE AT TOP OF BED LAYER
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDG(L,K) = ZBEDG(L,K-1) + HBED(L,K)
          ELSE
            ZBEDG(L,K) = HBED(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO L=2,LA
        ZBEDGT(L) = ZBEDG(L,KBT(L))
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            ZBEDC(L,K) = 0.5*(ZBEDG(L,K) + ZBEDG(L,K-1))
          ELSE
            ZBEDC(L,K) = 0.5*ZBEDG(L,K)
          ENDIF
        ENDDO
      ENDDO
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            DZBTR(L,K) = HBED(L,K)/(1. + VDRBED(L,K))
          ENDIF
        ENDDO
      ENDDO

      ! *** UPDATE POROSITY
      DO K=1,KB
        DO L=2,LA
          IF( K <= KBT(L) )THEN
            PORBED(L,K) = VDRBED(L,K)/(1. + VDRBED(L,K))
            PORBED1(L,K) = VDRBED1(L,K)/(1. + VDRBED1(L,K))
          ELSE
            PORBED(L,K) = 0.0
            PORBED(L,K) = 0.0
          ENDIF
        ENDDO
      ENDDO

      ! *** ADD OR REMOVE LAYERS
      !CALL CALBLAY
    ENDIF   ! *** END OF FINITE STRAIN OPION

    ! *********************************************************************
    ! *** UPDATE BULK DENSITY AND TOTAL SEDIMENT MASS
    !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ND,LF,LL,K,L,NS)
    DO ND=1,NDM
      LF = 2+(ND-1)*LDM
      LL = MIN(LF+LDM-1,LA)

      DO K=1,KB
        DO L=LF,LL
          SEDBT(L,K) = 0.0
          SNDBT(L,K) = 0.0
        ENDDO
      ENDDO

      IF( ISTRAN(6) >= 1 )THEN
        DO NS=1,NSED
          DO K=1,KB
            DO L=LF,LL
              SEDBT(L,K) = SEDBT(L,K) + SEDB(L,K,NS)
            ENDDO
          ENDDO
        ENDDO
      ENDIF

      IF( ISTRAN(7) >= 1 )THEN
        DO NS=1,NSND
          DO K=1,KB
            DO L=LF,LL
              SNDBT(L,K) = SNDBT(L,K) + SNDB(L,K,NS)
            ENDDO
          ENDDO
        ENDDO
      ENDIF

      DO K=1,KB
        DO L=LF,LL
          IF( HBED(L,K) > 0. )THEN
            ! *** UPDATE TOTAL/WET DENSITY
            BDENBED(L,K) = 1000.*PORBED(L,K) + 0.001*(SEDBT(L,K) + SNDBT(L,K))/HBED(L,K)
          ELSE
            BDENBED(L,K) = 0.
          ENDIF
        ENDDO
      ENDDO
    ENDDO  ! END OF DOMAIN
    !$OMP END PARALLEL DO

  ENDIF

  8011 FORMAT('CBED',I12,6E13.5)

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

  RETURN
  END

max_stars_repo_name: dsi-llc/EFDCPlus

max_forks_repo_name: dsi-llc/EFDCPlus


text:
	SUBROUTINE GCFILLF ( kx, ky, grid, ioffx, ioffy, iskip, nlvl,
     +			     clvl, icolr, linlbl, lintyp, zgrid, iret )
C************************************************************************
C* GCFILLF								*
C* 									*
C* This subroutine draws contours through a grid of data.  The 		*
C* algorithm used is based on a two-dimensional Lagrangian fitting of	*
C* the grid points.  It is the original GEMPAK contouring program.	*
C* The original subroutine, GCONTR, had a different calling sequence.	*
C* 									*
C* GCFILLF ( KX, KY, GRID, IOFFX, IOFFY, ISKIP, NLVL, CLVL, ICOLR,	*
C*           LINLBL, LINTYP, ZGRID, IRET )				*
C*									*
C* Input parameters:							*
C*									*
C*	KX		INTEGER		Number of x grid points 	*
C*	KY		INTEGER		Number of y grid points		*
C*	GRID (KX,KY)	REAL		Grid data array			*
C*	IOFFX		INTEGER		X offset to first point		*
C*	IOFFY		INTEGER		Y offset to first point		*
C*	ISKIP		INTEGER		Skip factor in original grid	*
C*	NLVL		INTEGER		Number of contour levels	*
C*	CLVL   (NLVL)	REAL		Contour level values		*
C*	ICOLR  (NLVL+1)	INTEGER		Contour color numbers		*
C*	LINLBL (NLVL+1)	INTEGER		Contour label types		*
C*	LINTYP (NLVL+1)	INTEGER		Contour fill types		*
C*                                                                      *
C* Work parameters:                                                     *
C*      ZGRID (KX*KY)   REAL            Array for grid points           *
C*									*
C* Output parameters:							*
C*									*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* M. desJardins/NMC	11/91	From GCONTR				*
C* S. Jacobs/NCEP	 1/96	Added CADJST to adjust grid values	*
C* T. Lee/GSC		 9/97	Fixed typo, NOMOMO -> NOMONO		*
C* S. Jacobs/NCEP	 2/00	Eliminate subboxing when CONTUR=1	*
C* T. Lee/SAIC		10/01	Added contour fill types to calling seq.*
C* S. Gilbert/NCEP      03/07   Added work arrays to call sequence      *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE 	'ERROR.PRM'
	INCLUDE		'CONTUR.CMN'
C*
	REAL		grid (*), clvl (*), zgrid (*)
	INTEGER		icolr (*), linlbl (*), lintyp (*)
C*
	LOGICAL		linflg, bad
	REAL		clev (LLCLEV)
	INTEGER		linwid (LLCLEV)
C------------------------------------------------------------------------
	iret = NORMAL
C
C*	Move the contour levels into a local array in order to add a
C*	maximum level.
C
	IF  ( nlvl .lt. LLCLEV )  THEN
	    nlev = nlvl
	  ELSE
	    nlev = LLCLEV - 1
	END IF
	bad = .false.
	DO  i = 1, nlev
	    clev (i) = clvl (i)
	    IF  ( i .gt. 1 )  THEN
		IF  ( clev (i) .le. clev (i-1) )  bad = .true.
	    END IF
	END DO
C
C*	If levels are not in the correct order, set error and return.
C
	IF  ( bad )  THEN
	    iret = NOMONO
	    RETURN
	END IF
C
C*	Build dummy maximum level.
C
	CALL GR_STAT  ( grid, kx, ky, 1, 1, kx, ky, gmin, gmax, gavg,
     +			gdev, ier )
	high = 2 * ABS ( gmax ) + ABS ( gmin )
	IF  ( high .le. clev (nlev) ) high = clev (nlev) + 1
	nlev = nlev + 1
	clev (nlev) = high
C
C*	Read in the box subset factor.  Check that it is valid and set
C*	it in common for the rest of the original contouring code.
C
	IF  ( jbxsub .lt. 2 )  THEN
	    numsub = 2
	  ELSE IF  ( jbxsub .gt. 48 )  THEN
	    numsub = 50
	  ELSE
	    numsub = jbxsub + 2
	END IF
C
C*	Move offsets and skip factor into common.
C
	offx = FLOAT ( ioffx )
	offy = FLOAT ( ioffy )
	skip = FLOAT ( iskip )
C
C*	This program expects the data in the Z array in common.
C
	DO  i = 1, kx * ky
	    zgrid (i) = grid (i)
	END DO
C
C*	Adjust the grid values.
C
	CALL CADJST ( kx, ky, nlev, clev, zgrid, ier )
C
C*	Save grid size as ISIZE, JSIZE.
C
	isize = kx
	jsize = ky
C
C*	Compute distance along each subbox.
C
	fincxy = 1. / FLOAT ( numsub - 1 )
C
C*	Call the driver for the fill contouring.
C
	linflg = .false.
	CALL FFDRIV  ( nlev, clev, icolr, lintyp, linwid, linlbl,
     +		       linflg, zgrid, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
! =============================================================================
!                      Write parcel diagnostics to NetCDF
! =============================================================================
module parcel_diagnostics_netcdf
    use constants, only : one
    use netcdf_utils
    use netcdf_writer
    use netcdf_reader
    use parcel_container, only : parcels, n_parcels
    use parcel_diagnostics
    use parameters, only : lower, extent, nx, ny, nz
    use config, only : package_version, cf_version
    use timer, only : start_timer, stop_timer
    use options, only : write_netcdf_options
    use physics, only : write_physical_quantities
    implicit none

    private

    character(len=512) :: ncfname
    integer            :: ncid
    integer            :: t_axis_id, t_dim_id, n_writes,            &
                          pe_id, ke_id, te_id, npar_id, nspar_id,   &
                          rms_x_vor_id, rms_y_vor_id, rms_z_vor_id, &
                          avg_lam_id, std_lam_id,                   &
                          avg_vol_id, std_vol_id, sum_vol_id,       &
                          psi_id

    double precision   :: restart_time

    integer :: parcel_stats_io_timer

    public :: create_netcdf_parcel_stats_file,  &
              write_netcdf_parcel_stats,        &
              parcel_stats_io_timer

    contains

        ! Create the parcel diagnostic file.
        ! @param[in] basename of the file
        ! @param[in] overwrite the file
        subroutine create_netcdf_parcel_stats_file(basename, overwrite, l_restart)
            character(*), intent(in)  :: basename
            logical,      intent(in)  :: overwrite
            logical,      intent(in)  :: l_restart
            logical                   :: l_exist

            ncfname =  basename // '_parcel_stats.nc'

            restart_time = -one
            n_writes = 1

            call exist_netcdf_file(ncfname, l_exist)

            if (l_restart .and. l_exist) then
                call open_netcdf_file(ncfname, NF90_NOWRITE, ncid)
                call get_num_steps(ncid, n_writes)
                call get_time(ncid, restart_time)
                call read_netcdf_parcel_stats_content
                call close_netcdf_file(ncid)
                n_writes = n_writes + 1
                return
            endif

            call create_netcdf_file(ncfname, overwrite, ncid)

            ! define global attributes
            call write_netcdf_info(ncid=ncid,                    &
                                   epic_version=package_version, &
                                   file_type='parcel_stats',     &
                                   cf_version=cf_version)

            call write_netcdf_box(ncid, lower, extent, (/nx, ny, nz/))

            call write_physical_quantities(ncid)

            call write_netcdf_options(ncid)

            call define_netcdf_temporal_dimension(ncid, t_dim_id, t_axis_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='pe',                                                  &
                long_name='potential energy',                               &
                std_name='',                                                &
                unit='m^5/s^2',                                             &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=pe_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='ke',                                                  &
                long_name='kinetic energy',                                 &
                std_name='',                                                &
                unit='m^5/s^2',                                             &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=ke_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='te',                                                  &
                long_name='total energy',                                   &
                std_name='',                                                &
                unit='m^5/s^2',                                             &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=te_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='psi',                                                 &
                long_name='enstrophy',                                      &
                std_name='',                                                &
                unit='m^3/s^2',                                             &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=psi_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='n_parcels',                                           &
                long_name='number of parcels',                              &
                std_name='',                                                &
                unit='1',                                                   &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=npar_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='n_small_parcel',                                      &
                long_name='number of small parcels',                        &
                std_name='',                                                &
                unit='1',                                                   &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=nspar_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='avg_lam',                                             &
                long_name='average aspect ratio',                           &
                std_name='',                                                &
                unit='1',                                                   &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=avg_lam_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='std_lam',                                             &
                long_name='standard deviation aspect ratio',                &
                std_name='',                                                &
                unit='1',                                                   &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=std_lam_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='avg_vol',                                             &
                long_name='average volume',                                 &
                std_name='',                                                &
                unit='m^3',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=avg_vol_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='std_vol',                                             &
                long_name='standard deviation volume',                      &
                std_name='',                                                &
                unit='m^3',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=std_vol_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='sum_vol',                                             &
                long_name='total volume',                                   &
                std_name='',                                                &
                unit='m^3',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=sum_vol_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='x_rms_vorticity',                                     &
                long_name='root mean square of x vorticity component',      &
                std_name='',                                                &
                unit='1/s',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=rms_x_vor_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='y_rms_vorticity',                                     &
                long_name='root mean square of y vorticity component',      &
                std_name='',                                                &
                unit='1/s',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=rms_y_vor_id)

            call define_netcdf_dataset(                                     &
                ncid=ncid,                                                  &
                name='z_rms_vorticity',                                     &
                long_name='root mean square of z vorticity component',      &
                std_name='',                                                &
                unit='1/s',                                                 &
                dtype=NF90_DOUBLE,                                          &
                dimids=(/t_dim_id/),                                        &
                varid=rms_z_vor_id)

            call close_definition(ncid)

        end subroutine create_netcdf_parcel_stats_file

        ! Pre-condition: Assumes an open file
        subroutine read_netcdf_parcel_stats_content

            call get_dim_id(ncid, 't', t_dim_id)

            call get_var_id(ncid, 't', t_axis_id)

            call get_var_id(ncid, 'pe', pe_id)

            call get_var_id(ncid, 'ke', ke_id)

            call get_var_id(ncid, 'te', te_id)

            call get_var_id(ncid, 'psi', psi_id)

            call get_var_id(ncid, 'n_parcels', npar_id)

            call get_var_id(ncid, 'n_small_parcel', nspar_id)

            call get_var_id(ncid, 'avg_lam', avg_lam_id)

            call get_var_id(ncid, 'std_lam', std_lam_id)

            call get_var_id(ncid, 'avg_vol', avg_vol_id)

            call get_var_id(ncid, 'std_vol', std_vol_id)

            call get_var_id(ncid, 'sum_vol', sum_vol_id)

            call get_var_id(ncid, 'x_rms_vorticity', rms_x_vor_id)

            call get_var_id(ncid, 'y_rms_vorticity', rms_y_vor_id)

            call get_var_id(ncid, 'z_rms_vorticity', rms_z_vor_id)

        end subroutine read_netcdf_parcel_stats_content

        ! Write a step in the parcel diagnostic file.
        ! @param[in] t is the time
        subroutine write_netcdf_parcel_stats(t)
            double precision, intent(in)    :: t

            call start_timer(parcel_stats_io_timer)

            if (t <= restart_time) then
                call stop_timer(parcel_stats_io_timer)
                return
            endif

            call open_netcdf_file(ncfname, NF90_WRITE, ncid)

            ! write time
            call write_netcdf_scalar(ncid, t_axis_id, t, n_writes)

            !
            ! write diagnostics
            !
            call write_netcdf_scalar(ncid, pe_id, pe, n_writes)
            call write_netcdf_scalar(ncid, ke_id, ke, n_writes)
            call write_netcdf_scalar(ncid, te_id, ke + pe, n_writes)
            call write_netcdf_scalar(ncid, psi_id, psi, n_writes)
            call write_netcdf_scalar(ncid, npar_id, n_parcels, n_writes)
            call write_netcdf_scalar(ncid, nspar_id, n_small, n_writes)
            call write_netcdf_scalar(ncid, avg_lam_id, avg_lam, n_writes)
            call write_netcdf_scalar(ncid, std_lam_id, std_lam, n_writes)
            call write_netcdf_scalar(ncid, avg_vol_id, avg_vol, n_writes)
            call write_netcdf_scalar(ncid, std_vol_id, std_vol, n_writes)
            call write_netcdf_scalar(ncid, sum_vol_id, sum_vol, n_writes)
            call write_netcdf_scalar(ncid, rms_x_vor_id, rms_zeta(1), n_writes)
            call write_netcdf_scalar(ncid, rms_y_vor_id, rms_zeta(2), n_writes)
            call write_netcdf_scalar(ncid, rms_z_vor_id, rms_zeta(3), n_writes)

            ! increment counter
            n_writes = n_writes + 1

            call close_netcdf_file(ncid)

            call stop_timer(parcel_stats_io_timer)

        end subroutine write_netcdf_parcel_stats
end module parcel_diagnostics_netcdf

max_stars_repo_name: matt-frey/epic

max_issues_repo_name: matt-frey/epic

max_forks_repo_name: matt-frey/epic


text:
      subroutine hbbdecay(p,ib,ibb,msq)
************************************************************************
*     Author: J.M. Campbell, June 2012                                 *
*                                                                      *
*     matrix element squared for the process of                        *
*     Higgs decay  H --> b(ib)+b~(ibb)                                 *
*     with bottom mass included                                        *
************************************************************************
      implicit none
      include 'constants.f'
      include 'masses.f'
      integer ib,ibb
      double precision p(mxpart,4),s56,msq,msqhbb

      s56=2d0*(p(ib,4)*p(ibb,4)-p(ib,1)*p(ibb,1)
     &        -p(ib,2)*p(ibb,2)-p(ib,3)*p(ibb,3))+2d0*mb**2

      msq=msqhbb(s56)

      return
      end

      double precision function msqhbb(s)
      implicit none
      double precision s
      include 'constants.f'
      include 'masses.f'
      include 'couple.f'
      include 'part.f'
      include 'msbarmasses.f'
      include 'ewcouple.f'
      include 'first.f'
      double precision mb_eff,massfrun
      save mb_eff
!$omp threadprivate(mb_eff)

      if (first) then
c--- run mb to appropriate scale
        if (part .eq. 'lord') then
          mb_eff=massfrun(mb_msbar,hmass,amz,1)
        else
          mb_eff=massfrun(mb_msbar,hmass,amz,2)
        endif
        first=.false.
      endif

      msqhbb=xn*gwsq*mb_eff**2/(4d0*wmass**2)*2d0*(s-4d0*mb**2)

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20100628 10:38:32 nbsp Howdy Elizabeth, and Welcome to the Wiki.  My names Evan, pleased to meet you!  I saw you mentioned you had a standard poodle; they are wonderful dogs.  My wife is a staunch Labrador retriever fan though, so theres probably not one in the near future for us.   Once again, welcome to the wiki!  Users/JabberWokky Evan JabberWokky Edwards

20100805 23:13:44 nbsp welcome to the wiki, have you heard about 3 degrees of Davis? Users/StevenDaubert

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module GwfNpfOptionsModule
  use KindModule, only: I4B, DP
  use ConstantsModule, only: DONE

  implicit none
  private

  !> Data structure and helper methods for passing NPF options
  !! into npf_df, as an alternative to reading those from file
  !<
  type, public :: GwfNpfOptionsType
    integer(I4B) :: icellavg    !< same as npf variable
    integer(I4B) :: ithickstrt  !< same as npf variable
    integer(I4B) :: iperched    !< same as npf variable
    integer(I4B) :: ivarcv      !< same as npf variable
    integer(I4B) :: idewatcv    !< same as npf variable
    integer(I4B) :: irewet      !< same as npf variable
    real(DP)     :: wetfct      !< same as npf variable
    integer(I4B) :: iwetit      !< same as npf variable
    integer(I4B) :: ihdwet      !< same as npf variable
  contains
    procedure, pass(this) :: construct
    procedure, pass(this) :: destroy
  end type GwfNpfOptionsType

contains

  !> @brief construct the input options, set variables to their defaults
  !<
  subroutine construct(this)
    class(GwfNpfOptionsType), intent(inout) :: this  !< the NPF options, as in the input OPTIONS block

    this%icellavg = 0
    this%ithickstrt = 0
    this%iperched = 0
    this%ivarcv = 0
    this%idewatcv = 0
    this%irewet = 0
    this%wetfct = DONE
    this%iwetit = 1
    this%ihdwet = 0

  end subroutine construct

  !> @brief cleans up
  !<
  subroutine destroy(this)
    class(GwfNpfOptionsType), intent(inout) :: this  !< the NPF options, as in the input OPTIONS block

    ! nothing to be done here for now...

  end subroutine destroy

end module GwfNpfOptionsModule


! { dg-do compile }
!
! PR 79311: [OOP] ICE in generate_finalization_wrapper, at fortran/class.c:1992
!
! Contributed by DIL <liakhdi@ornl.gov>

module tensor_recursive
  implicit none

  type :: tens_signature_t
  contains
    final :: tens_signature_dtor
  end type

  type :: tens_header_t
    type(tens_signature_t) :: signature
  contains
    final :: tens_header_dtor
  end type

contains

  subroutine tens_signature_dtor(this)
    type(tens_signature_t) :: this
  end subroutine

  subroutine tens_header_dtor(this)
    type(tens_header_t) :: this
  end subroutine

end


text:
C***********************************************************************
C    Module:  symgen.f  (part of font generation tools in Xplot/sym)
C
C    Copyright (C) 1996 Harold Youngren, Mark Drela
C
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C
C    Report problems to:    guppy@maine.com
C                        or drela@mit.edu
C***********************************************************************

      PROGRAM SYMGEN
C------------------------------------------------------------------------
C
C     Interactive program for creating and modifying vector fonts.
C
C
C     A font set is kept stored in a database file XXXX.FNT, which
C     is read in when the program is invoked with
C
C      % symgen XXXX
C
C     The first line of XXXX.FNT is a list of the defined characters.
C     Each subsequent line then defines one character in the list:
C
C  ABCD...
C  21616 12168 13280 ...      <- definition for "A"
C  21616 12280 14880 ...      <- definition for "B"
C   .
C   .
C
C     Each 5-digit number defines a vector node x,y location on a 96x96
C     grid, with 16,16 being the lower-left "origin" of the character.
C     See SUBR. SYMPLT header below for the encoding format.
C     If the numbers are zero or missing, then the font character
C     is not plotted (blank).
C
C     Once the database has been edited, command "W" will generate
C     a Fortran INCLUDE file XXXX.INC, with the vector font information
C     placed in DATA statements.  This is for use in a SYMBOL-type
C     routine to allow plotting of the font.  Four sample databases,
C     include files, and routines already exist:
C
C       CHAR.FNT    CHAR.INC    SUBR. PLCHAR    (Upper,lower case letters)
C       SLAN.FNT    SLAN.INC    SUBR. PLSLAN    (Slanted version of CHAR )
C       MATH.FNT    MATH.INC    SUBR. PLMATH    (Latex-like Greek symbols)
C       SYMB.FNT    SYMB.INC    SUBR. PLSYMB    (Versatec plotting symbols)
C
C             also available:   SUBR. PLNUMB    (Whole floating-point numbers)
C
C     SUBR. PLCHAR, PLSLAN, and PLMATH (in ../plt_font.f) are identical
C     except for their INCLUDE statements.  One can generate analogous
C     files and plot routines for any custom font, e.g.
C
C       USER.FNT    USER.INC    SUBR. PLUSER
C
C
C     SUBR. PLSYMB is nearly the same, except that it takes an integer
C     argument ISYM = 0,1... to select the symbol defined on line ISYM+1:
C
C  ABCD...
C  21616 12168 13280 ...      <- definition for ISYM = 0
C  21616 12280 14880 ...      <- definition for ISYM = 1
C   .
C   .
C     The ASCII characters on the first line are ignored in PLSYMB.
C
C     PLSYMB also shifts the symbol down and left by 1/2 the symbol
C     size, so the symbol is centered on 0,0.  This duplicates the
C     Versatec SYMBOL routine.  Note that PLSYMB can only plot one
C     symbol at a time.
C
C
C     WARNING:  The CHAR font implemented in PLCHAR is assumed to be
C               a "WYSIWIG" font, in which each font character closely
C               represents its ASCII index character.  It is used by
C               other libPlt routines (and also this program!), and
C               hence should not be significantly modified.
C
C
C------------------------------------------------------------------------
C
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
      INTEGER NODET(20)
C
      LOGICAL SAVED
C
      CHARACTER*80 ARGP1, STRING, RANGE
C
      CHARACTER*1 COPT, ANS, ALPH, ALPHO, ALPHT, KCHAR
      LOGICAL LABORT, LDONE
C
      SIZE = 6.0
      IDEV = 1
      IPSLU = 0
C
      CH = 0.7
      X0 = 0.0
      Y0 = 0.0
C
      NCHARS = 0
C
C---- get the XXXX working set name
      CALL GETARG(1,ARGP1)
      CALL LOAD(ARGP1)
C
C---- initial overlay character, plot-test string
      ALPHO = ' '
      STRING = ' '
C
C---- initial target character is first one in database
      ALPH = CHARS(1:1)
      SAVED = .TRUE.
C
      CALL PLINITIALIZE
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
C
      CALL GETCOLOR(ICOL0)
C
 1000 FORMAT(A)
C
 900  CONTINUE
      WRITE(*,1005)
 1005 FORMAT(/'   C hange target character'
     &       /'   I nput  target character vectors'
     &       /'   M odify target character vectors'
     &       /'   T ranslate target character'
     &       /'   O verlay other character'
     &       /'   L ean character(s)'
     &       /'   P lot specified character string'
     &       /'   A dd a new character to database'
     &       /'   D elete a character from database'
     &       /'   S ave database'
     &       /'   W rite include file from current database'
     &       /'   Z oom'
     &       /'   U nzoom'
     &       /'   Q uit')
C
 910  WRITE(*,1015) ALPH
 1015 FORMAT(/' Select option ( target character = ',A1,' ):  ', $)
      READ(*,1000) COPT
C
      IF(INDEX('Qq',COPT).NE.0) THEN
        IF(.NOT.SAVED) THEN
          WRITE(*,*)
          WRITE(*,*) 'Database not saved.  Really quit ?  N'
          READ(*,1000) ANS
          IF(INDEX('Yy',ANS) .EQ. 0) GO TO 900
        ENDIF
        CALL PLCLOSE
        STOP
      ENDIF
C
      NOPT = INDEX('CIMTOLPADSWZU',COPT) + INDEX('cimtolpadswzu',COPT)
C
      GOTO (10,20,30,40,50,60,70,80,90,100,110,120,130) NOPT
      GOTO 900
C
C===========================================
C---- select character
 10   CONTINUE
      WRITE(*,1010) CHARS(1:NCHAR)
 1010 FORMAT(/1X,A,
     &      //1X,'Select new target character: ',$)
      READ(*,1000) ALPH
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- input character
 20   CONTINUE
C
      KC = INDEX(CHARS,ALPH)
      IF(KC.EQ.0) THEN
        WRITE(*,*) 'Character is not in database.'
        GO TO 900
      ENDIF
C
C---- clear node-input buffer array
      DO 201 K=1, 20
        NODET(K) = 0
 201  CONTINUE
C
C
C---- small node-select symbol size
      CHS = 0.010*CH
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PABORT
      CALL PDONE
C
ccc      CALL PLCHAR(X0,Y0,CH,ALPH,0.0,1)
C
      CALL PLFLUSH
C
      WRITE(*,*)
      WRITE(*,*) 'Input stroke points (double-click to end a stroke)...'
C
      ILAST = -999
      JLAST = -999
      ISTAT = 2
C
      DO 205 K=1, 20
 204    CALL GETCURSORXY(XC,YC,KCHAR)
C
C------ exit if clicked inside "ABORT" box
        IF(LABORT(XC,YC)) THEN
         WRITE(*,*) 'Database unchanged.'
         GO TO 209
        ENDIF
C
C------ finish up if clicked inside "DONE" box
        IF(LDONE(XC,YC)) GO TO 206
C
C------ set integer coordinates and clip to within field
        I = INT(64.0*(XC-X0)/CH + 16.5)
        J = INT(64.0*(YC-Y0)/CH + 16.5)
        I = MAX( 0 , MIN(I,96) )
        J = MAX( 0 , MIN(J,96) )
C
        IF(ILAST.EQ.I .AND. JLAST.EQ.J) THEN
          ISTAT = 2
          ILAST = -999
          JLAST = -999
          WRITE(*,*) 'Will start new stroke...'
          GO TO 204
        ENDIF
C
        WRITE(*,1026) K, I, J
 1026   FORMAT(1X,'   node',I3,':   ', 2I5)
C
C------ set new node location in buffer array
        NODET(K) = 10000*ISTAT + 100*I + J
C
C------ put small symbol over actual node location
        XNODE = CH*FLOAT(I-16)/64.0
        YNODE = CH*FLOAT(J-16)/64.0
        IF(ILAST.EQ.-999 .OR. JLAST.EQ.-999) THEN
          CALL PLOT(XNODE,YNODE,3)
        ELSE
          CALL PLOT(XNODE,YNODE,2)
        ENDIF
        CALL NEWCOLORNAME('orange')
        CALL PLSYMB(XNODE,YNODE,CHS,1,0.0,0)
        CALL NEWCOLOR(ICOL0)
        CALL PLFLUSH
C
C------ set up for next input point
        ILAST = I
        JLAST = J
        ISTAT = 1
C
 205  CONTINUE
      WRITE(*,*) '20-node array limit reached'
 206  CONTINUE
C
C---- put buffer array into database
      KC = INDEX(CHARS,ALPH)
      DO 207 K=1, 20
        NODE(K,KC) = NODET(K)
 207  CONTINUE
      WRITE(*,*)
      WRITE(*,*) 'Database updated.'
      SAVED = .FALSE.
C
C---- replot whole character
 209  CALL PLEND
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
C
      CALL PLEND
      GO TO 900
C
C===========================================
C---- modify character
 30   CONTINUE
      WRITE(*,1030) ALPH
 1030 FORMAT(/1X, 'Modifying character: ', A)
C
      KC = INDEX(CHARS,ALPH)
      IF(KC.EQ.0) THEN
        WRITE(*,*) 'Character is not in database.'
        GO TO 900
      ENDIF
C
C---- save database node array for restoration after abort
      DO 301 K=1, 20
        NODET(K) = NODE(K,KC)
 301  CONTINUE
C
      CALL PLTINI
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
C
ccc      CALL PLSYMB(X0,Y0,CH,ALPH,0.0,1)
c
C
C---- small node-select symbol size
      CHS = 0.010*CH
C
      CALL PABORT
      CALL PDONE
      CALL PLFLUSH
C
      WRITE(*,*)
      WRITE(*,*) 'Click on old/new point pairs...'
      DO 305 IPASS=1, 12345
C
        CALL GETCURSORXY(XC,YC,KCHAR)
C
C------ restore and exit if clicked inside "ABORT" box
        IF(LABORT(XC,YC)) THEN
         DO 3052 K=1, 20
           NODE(K,KC) = NODET(K)
 3052    CONTINUE
         WRITE(*,*) 'Database unchanged.'
         GO TO 309
        ENDIF
C
C------ finish up if clicked inside "DONE" box
        IF(LDONE(XC,YC)) GO TO 306
C
C------ set integer coordinates and clip to within field
        I = INT(64.0*(XC-X0)/CH + 16.5)
        J = INT(64.0*(YC-Y0)/CH + 16.5)
        I = MAX( 0 , MIN(I,96) )
        J = MAX( 0 , MIN(J,96) )
C
C------ find nearest vector node
        IDMIN = 1000000
        KMIN = 0
        DO 3054 K=1, 20
          ISTAT = NODE(K,KC) / 10000
          IF(ISTAT.EQ.0) GO TO 3055
C
          NODEB = NODE(K,KC) - ISTAT*10000
          IT = NODEB / 100
          JT = NODEB - 100*IT
C
          IDIST = (I - IT)**2 + (J - JT)**2
          IF(IDIST.LT.IDMIN) THEN
            IDMIN = IDIST
            KMIN = K
          ENDIF
 3054   CONTINUE
 3055   CONTINUE
C
        KT = KMIN
        IF(KT.EQ.0) THEN
          WRITE(*,*) 'Nearest point not found.'
          GO TO 900
        ENDIF
C
C------ plot small symbol on nearest vector node to identify it
        ISTAT = NODE(KT,KC) / 10000
        NODEB = NODE(KT,KC) - ISTAT*10000
        IT = NODEB / 100
        JT = NODEB - 100*IT
C
        CALL NEWCOLORNAME('red')
        XNODE = CH*FLOAT(IT-16)/64.0
        YNODE = CH*FLOAT(JT-16)/64.0
        CALL PLSYMB(XNODE,YNODE,CHS,1,0.0,0)
        FK = FLOAT(KT)
        CALL PLNUMB(XNODE+1.5*CHS,YNODE+1.5*CHS,3.0*CHS,FK,0.0,-1)
        CALL NEWCOLOR(ICOL0)
C
C------ now get new location for vector node KT
        CALL GETCURSORXY(XC,YC,KCHAR)
C
        IF(LABORT(XC,YC)) THEN
         DO 3056 K=1, 20
           NODE(K,KC) = NODET(K)
 3056    CONTINUE
         WRITE(*,*) 'No changes made.'
         GO TO 309
        ENDIF
C
        IF(LDONE(XC,YC)) GO TO 306
C
C------ integer coordinates for new location
        I = INT(64.0*(XC-X0)/CH + 16.5)
        J = INT(64.0*(YC-Y0)/CH + 16.5)
        I = MAX( 0 , MIN(I,96) )
        J = MAX( 0 , MIN(J,96) )
C
C------ encode coordinates into database
        NODE(KT,KC) = 10000*ISTAT + 100*I + J
C
C------ replot everything
        CALL PLEND
        CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
        CALL PABORT
        CALL PDONE
        CALL PLFLUSH
C
 305  CONTINUE
 306  CONTINUE
C
      WRITE(*,*)
      WRITE(*,*) 'Database updated.'
      SAVED = .FALSE.
C
C---- replot final state
 309  CALL PLEND
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
C
      CALL PLEND
      GO TO 900
C
C===========================================
C---- translate current character
 40   CONTINUE
C
      DO K=1, 20
        NODET(K) = 0
      ENDDO
C
 41   WRITE(*,1040)
 1040 FORMAT(/' Enter dX, dY (in points):  ', $)
      READ(*,*,ERR=41) ID, JD
C
      IC = INDEX(CHARS,ALPH)
      DO 42 K=1, 20
C
C------ unpack node coordinates from database
        ISTAT = NODE(K,IC) / 10000
        IF(ISTAT.EQ.0) GO TO 43
C
        NODEB = NODE(K,IC) - ISTAT*10000
        I = NODEB / 100
        J = NODEB - 100*I
C
C------ shift coordinates
        I = I + ID
        J = J + JD
C
        IF(I.LT.0 .OR. I.GT.96) THEN
          WRITE(*,*) 'X movement puts character outside field 0..96'
          GO TO 40
        ENDIF
C
        IF(J.LT.0 .OR. J.GT.96) THEN
          WRITE(*,*) 'Y movement puts character outside field 0..96'
          GO TO 40
        ENDIF
C
C------ encode new coordinates in buffer array
        NODET(K) = 10000*ISTAT + 100*I + J
C
 42   CONTINUE
 43   CONTINUE
C
C---- store shifted nodes in database from buffer array
      DO 46 K=1, 20
        NODE(K,IC) = NODET(K)
 46   CONTINUE
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- select overlay character
 50   CONTINUE
      ALPHO = ' '
      WRITE(*,1050) CHARS(1:NCHAR)
 1050 FORMAT(/1X,A,
     &      //1X,'Select character to overlay (<return> if none: ',$)
      READ(*,1000) ALPHO
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- lean character(s)
 60   CONTINUE
      WRITE(*,1060)
 1060 FORMAT(/' Enter character range to lean (e.g. AZ): ',$)
      READ(*,1000) RANGE
      IF(INDEX(RANGE,' ').EQ.1) GO TO 900
C
 62   WRITE(*,1062)
 1062 FORMAT(/1X,'Enter tan(lean_angle), + to right:  ',$)
      READ (*,*,ERR=62) TANA
C
      CALL LEAN(RANGE,TANA)
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- enter plot test string
 70   CONTINUE
      WRITE(*,1070)
 1070 FORMAT(/' Enter string: ',$)
      READ(*,1000) STRING
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- add character to database
 80   CONTINUE
      IF(NCHAR+1.GE.100) THEN
        WRITE(*,*)
        WRITE(*,*) 'Array limit reached.  Cannot add character.'
        GO TO 900
      ENDIF
C
      WRITE(*,1080)
 1080 FORMAT(/1X,'Enter character to be added: ',$)
      READ(*,1000) ALPH
      IC = INDEX(CHARS,ALPH)
      IF(IC.NE.0) THEN
        WRITE(*,*) 'That is already in database.'
        GO TO 900
      ENDIF
C
 81   WRITE(*,1081) CHARS(1:NCHAR), ALPH
 1081 FORMAT(/1X,A,
     &      //1X,' ... insert ',A1,
     &       ' before which character ? (<return> to append) : ',$)
      READ(*,1000) ALPHT
      IC = INDEX(CHARS,ALPHT)
      IF(IC.EQ.0) IC = NCHAR+1
C
C---- move up all characters above IC to make room for new character
      IF(IC.LE.NCHAR) THEN
        CHARS(IC+1:NCHAR+1) = CHARS(IC:NCHAR)
        DO 802 KC=NCHAR, IC, -1
          DO 8024 K=1, 20
            NODE(K,KC+1) = NODE(K,KC)
 8024     CONTINUE
 802    CONTINUE
      ENDIF
C
      NCHAR = NCHAR+1
C
C---- set new character and clear its vectors
      CHARS(IC:IC) = ALPH
      DO 804 K=1, 20
        NODE(K,IC) = 0
 804  CONTINUE
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- delete character from database
 90   CONTINUE
C
      WRITE(*,1090)
 1090 FORMAT(/1X,'Select character to be deleted: ',$)
      READ(*,1000) ALPH
C
      IF(INDEX(' ',ALPH).EQ.1) THEN
       WRITE(*,*) 'No action taken'
       GO TO 900
      ENDIF
C
      IC = INDEX(CHARS,ALPH)
C
      IF(IC.EQ.0) THEN
       WRITE(*,*) 'Character not in database.'
       GO TO 900
      ENDIF
C
C---- pull down all characters above the one to be deleted
      CHARS(IC:NCHAR) = CHARS(IC+1:NCHAR+1)
      DO 902 KC=IC, NCHAR-1
        DO 9024 K=1, 20
          NODE(K,KC) = NODE(K,KC+1)
 9024   CONTINUE
 902  CONTINUE
      NCHAR = NCHAR-1
      IC = MIN(IC,NCHAR)
C
C---- set new current character
      ALPH = CHARS(IC:IC)
C
      CALL PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      CALL PLFLUSH
      CALL PLEND
      GO TO 900
C
C===========================================
C---- write out database
 100  CONTINUE
      CALL SAVE(ARGP1,SAVED)
      GO TO 900
C
C===========================================
C---- write out include file
 110  CONTINUE
      CALL WRIT(ARGP1)
      GO TO 910
C
C===========================================
C---- set zoom
 120  CONTINUE
      CALL USETZOOM(.FALSE.,.TRUE.)
      CALL REPLOT(1)
      GO TO 910
C
C===========================================
C---- clear zoom
 130  CONTINUE
      CALL CLRZOOM
      CALL REPLOT(1)
      GO TO 910
C
      END

      SUBROUTINE LOAD(ARGP1)
      CHARACTER*80 ARGP1,FNAME
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
 1000 FORMAT(A)
C
      K = INDEX(ARGP1,' ') - 1
      IF(K.LT.1) THEN
       WRITE(*,*) 'SYMGEN argument error: ', ARGP1
       STOP
      ENDIF
C
      FNAME = ARGP1(1:K) // '.FNT'
C
      OPEN(3,FILE=FNAME,STATUS='OLD',ERR=90)
      READ(3,1000) CHARS
      NCHAR = INDEX(CHARS,' ') - 1
C
      DO IC=1, NCHAR
        DO K=1, 20
          NODE(K,IC) = 0
        ENDDO
      ENDDO
C
      DO IC=1, NCHAR
        READ(3,*,END=80) (NODE(K,IC),K=1,20)
C
C------ clear any invalid points (should be zero anyway)
        DO K=1, 20
          ISTAT = NODE(K,IC) / 10000
          IF(ISTAT.EQ.0) NODE(K,IC) = 0
        ENDDO
      ENDDO
C
 80   CONTINUE
      CLOSE(3)
C
      WRITE(*,1200) CHARS(1:NCHAR)
 1200 FORMAT(/1X,'Database read for the following character set...'
     &      //1X, A)
C
      RETURN
C
 90   CONTINUE
      WRITE(*,*) 'Database file not found:  ',FNAME(1:60)
      RETURN
C
      END

      SUBROUTINE SAVE(ARGP1,SAVED)
      CHARACTER*80 ARGP1,FNAME
      CHARACTER*1 ANS
      CHARACTER*4 CNUM
      LOGICAL SAVED
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
C
 1000 FORMAT(A)
 1500 FORMAT(1X,20I6)
 2100 FORMAT(/' File ',A,' exits.  Overwrite ?  Y')
C
C
      KARG = INDEX(ARGP1,' ') - 1
C
C---- write new database file for next SYMGEN call
      FNAME = ARGP1(1:KARG) // '.FNT'
      OPEN(3,FILE=FNAME,STATUS='OLD',ERR=10)
C
      WRITE(*,2100) FNAME(1:KARG+4)
      READ(*,1000) ANS
      IF(INDEX('Nn',ANS) .NE. 0) THEN
        CLOSE(3)
        RETURN
      ENDIF
C
 10   OPEN(3,FILE=FNAME,STATUS='UNKNOWN')
C
      WRITE(3,1000) CHARS(1:NCHAR)
      DO IC=1, NCHAR
        WRITE(3,1500) (NODE(K,IC),K=1,20)
      ENDDO
      CLOSE(3)
      WRITE(*,*)
      WRITE(*,*) 'New database file written: ', FNAME(1:KARG+4)
      SAVED = .TRUE.
C
      RETURN
      END

      SUBROUTINE WRIT(ARGP1)
      CHARACTER*80 ARGP1,FNAME,LINE
      CHARACTER*1 ANS
      CHARACTER*4 CNUM
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
C---- write new include file for SYMBOL-type routine
C
 1000 FORMAT(A)
 1500 FORMAT(1X,20I6)
 2100 FORMAT(/' File ',A,' exits.  Overwrite ?  Y')
C
      KARG = INDEX(ARGP1,' ') - 1
C
      FNAME = ARGP1(1:KARG) // '.INC'
      OPEN(4,FILE=FNAME,STATUS='OLD',ERR=20)
C
      WRITE(*,2100) FNAME(1:KARG+4)
      READ(*,1000) ANS
      IF(INDEX('Nn',ANS) .NE. 0) THEN
        CLOSE(4)
        GO TO 28
      ENDIF
C
 20   OPEN(4,FILE=FNAME,STATUS='UNKNOWN')
      WRITE(4,*) ' '
      WRITE(4,3010) NCHAR, NCHAR
C
      WRITE(4,*) ' '
      NARR = (NCHAR-1)/10
      DO N=0, NARR
        NDIM = MIN( 10 , NCHAR-N*10)
        WRITE(4,3201) N, NDIM
      ENDDO
      DO N=0, NARR
        WRITE(4,3202) 10*N+1, N
      ENDDO
C
C---Write character translation data statements
      WRITE(4,3020) NCHAR
      DO L=1, 8
        K1 = 26*L - 25
        K2 = 26*L
        IF(K1 .GT. NCHAR) GO TO 21
        K2 = MIN(K2,NCHAR)
        WRITE(4,3021) K1,K2, CHARS(K1:K2)
      ENDDO
 21   CONTINUE
C
C---Write out character node data
      DO N=0, NARR
        IC1 = 10*N + 1
        IC2 = 10*N + 10
        IC2 = MIN(IC2,NCHAR)
C
        IC = IC1
        WRITE(4,3210) N
        DO IC=IC1, IC2-1
          WRITE(4,3220) (NODE(K,IC),K= 1,10)
          WRITE(4,3220) (NODE(K,IC),K=11,20)
        ENDDO
        IC = IC2
        WRITE(4,3220) (NODE(K,IC),K= 1,10)
        WRITE(4,3230) (NODE(K,IC),K=11,20)
      ENDDO
C
      CLOSE(4)
C
 3010 FORMAT('      CHARACTER*',I3,' CHARS'
     &      /'      INTEGER NODE(20,',I3,')')
C
 3201 FORMAT('      DIMENSION NODE',I1,'(20,',I2,')')
 3202 FORMAT('      EQUIVALENCE ( NODE(1,',I3,') , NODE',I1,'(1,1) )')
C
 3020 FORMAT('      DATA NCHARS / ',I3,' /')
 3021 FORMAT('      DATA CHARS(',I3,':',I3,') / ''',A,''' /')
C
 3210 FORMAT( '      DATA NODE',I1,' /')
 3220 FORMAT( '     & ',10(I5,',') )
 3230 FORMAT( '     & ', 9(I5,','),I5,' /' )
C
      CLOSE(4)
      WRITE(*,*)
      WRITE(*,*) 'New include file written: ', FNAME(1:KARG+4)
      WRITE(*,*) 'Put include file in Xplot directory,',
     &          '  make libPlt.a  to implement new font.'
C
 28   CONTINUE
      RETURN
C
      END

      SUBROUTINE LEAN(RANGE,TANA)
      CHARACTER*(*) RANGE
C
      CHARACTER*1 ALPH
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
      KC1 = INDEX(CHARS,RANGE(1:1))
      KC2 = INDEX(CHARS,RANGE(2:2))
C
      IF(KC1.EQ.0 .OR. KC2.EQ.0) THEN
       WRITE(*,*) 'Specified range not in current character set.'
       WRITE(*,*) 'No action taken.'
       RETURN
      ENDIF
C
C---- go over each character...
      DO 12 KC=KC1, KC2
C
C------ go over each node
        DO K=1, 20
C
C-------- strip off leading point-status digit
          ISTAT = NODE(K,KC) / 10000
          NODEB = NODE(K,KC) - ISTAT*10000
C
          IF(ISTAT.GT.0) THEN
C
C---------- decode x and y location coordinates
            I = NODEB / 100
            J = NODEB - 100*I
C
C---------- perform tilt
            I = I + INT( FLOAT(J-16)*TANA )
C
C---------- encode coordinates back into database
            NODE(K,KC) = 10000*ISTAT + 100*I + J
C
          ENDIF
C
        ENDDO
C
 12   CONTINUE
C
      RETURN
      END

      SUBROUTINE PLTINI
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
C
      RELWSIZ = 0.8
      CALL PLOPEN(RELWSIZ,IPSLU,IDEV)
      CALL PLOTABS(0.125,0.125,-3)
      CALL NEWFACTOR(SIZE)
      CALL PLOT(0.25*CH,0.25*CH,-3)
C
      RETURN
      END ! PLTINI

      SUBROUTINE PLTALL(X0,Y0,ALPH,ALPHO,STRING)
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
      CHARACTER*1 ALPH, ALPHO
      CHARACTER*80 STRING
C
      CHSEQ = 0.10*CH
      CHSTR = 0.05*CH
C
      XSEQ = CH + CH/4.0 + 0.03*CH
      YSEQ = CH + CH/4.0 - 1.10*CHSEQ
C
      XCHR = CH + CH/4.0 + 0.03*CH
      YCHR = CH          - 0.07*CH
C
      XSTR = XCHR
      YSTR = YCHR - 0.65*CH
C
      CALL PLTINI
      CALL PLTGRD(X0,Y0,CH)
      CALL SYMPLT(X0,Y0,CH,ALPH,0.0,-1)
      CALL SYMPLT(X0,Y0,CH,ALPHO,0.0,1)
      CALL SYMPLT(XSTR,YSTR,CHSTR,STRING,0.0,80)
      CALL PLTCHR(XCHR,YCHR)
C
      XQ  = XSEQ
      CHQ = CHSEQ
c      CALL PLOT(XQ,YSEQ,3)
      DO ISCAL=0, 10
        CALL SYMPLT(XQ,YSEQ,CHQ,ALPH,0.0,1)
c        CALL PLCHAR(999.,999.,CHQ,ALPH,0.0,1)
        XQ  = XQ + CHQ
        CHQ = CHQ * 0.75
      ENDDO
C
      RETURN
      END

      SUBROUTINE PLTCHR(XCHR,YCHR)
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
C
      CHARACTER*1 ALPH
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
      CHP = 0.040*CH
C
      DO 10 ILIN=1, 5
        YPLT = YCHR - CHP - 2.5*CHP*(ILIN-1)
        DO 102 KC=1, 26
          XPLT = XCHR + 0.5*CHP + 1.0*CHP*(KC-1)
          XPLT = XCHR + 1.0*CHP*(KC-1)
C
          IC = KC + 26*(ILIN-1)
          IF(IC.GT.NCHAR) GO TO 15
C
          ALPH = CHARS(IC:IC)
          CALL SYMPLT(XPLT,YPLT,CHP,ALPH,0.0,1)
 102    CONTINUE
 10   CONTINUE
C
 15   CONTINUE
      RETURN
      END

      SUBROUTINE SYMPLT(XC,YC,CH,STRING,ANGLE,NC1)
C----------------------------------------------------------------
C     Plots character string using vector font database NODE(..)
C
C     Each NODE(..) value has the form  sxxyy
C
C        xx   = location of polyline point   0 ... 96
C        yy
C        s    = 2  if point is the start of a new polyline stroke
C             = 1  if point is inside or at the end of stroke
C             = 0  if point is not valid
C----------------------------------------------------------------
      CHARACTER*(*) STRING
C
      CHARACTER*1 ALPH
C
      CHARACTER*100 CHARS
      COMMON /CDATA/ CHARS
      COMMON /NDATA/ NCHAR, NODE(20,100)
C
      NC = ABS(NC1)
      CHS = 0.010*CH
C
      CALL GETCOLOR(ICOL0)
C
C---- go over each character...
      DO 12 IC=1, NC
C
C------ set plot location
        X0 = XC + CH*FLOAT(IC-1)
        Y0 = YC
C
        ALPH = STRING(IC:IC)
        KC = INDEX(CHARS,ALPH)
C
        IF(KC.NE.0) THEN
C
C-------- decode and plot each node
          DO K=1, 20
C
C---------- strip off leading point-status digit
            ISTAT = NODE(K,KC) / 10000
            NODEB = NODE(K,KC) - ISTAT*10000
C
C---------- decode x and y location coordinates
            I = NODEB / 100
            J = NODEB - 100*I
C
            X = X0 + CH*FLOAT(I-16)/64.0
            Y = Y0 + CH*FLOAT(J-16)/64.0
C
            IF     (ISTAT.EQ.0) THEN
             GOTO 12
            ELSEIF (ISTAT.EQ.1) THEN
             CALL PLOT(X,Y,2)
            ELSE
             CALL PLOT(X,Y,3)
            ENDIF
C
C---------- plot symbol at vector nodes
            IF(NC1.LT.0) THEN
              CALL NEWCOLORNAME('blue')
              CALL PLSYMB(X,Y,CHS,1,0.0,0)
              FK = FLOAT(K)
              CALL PLNUMB(X+1.5*CHS,Y+1.5*CHS,3.0*CHS,FK,0.0,-1)
              CALL NEWCOLOR(ICOL0)
              CALL PLOT(X,Y,3)
            ENDIF
C
          ENDDO
        ENDIF
C
 12   CONTINUE
C
      RETURN
      END

      SUBROUTINE PLTGRD(X0,Y0,CH)
      DATA LMASK1, LMASK2, LMASK3 / -32640, -30584, -21846 /
C
      CALL NEWPEN(1)
C
      NX = 24
      NY = 24
      DX = CH/16.0
      DY = CH/16.0
      CALL PLGRID(X0-0.25*CH,Y0-0.25*CH, NX,DX, NY,DY, LMASK2 )
C
      NX = 2
      NY = 2
      DX = CH/2.0
      DY = CH/2.0
      CALL PLGRID(X0,Y0, NX,DX, NY,DY, LMASK3 )
C
      DO I=0,2
        DO J=0,2
          XPLT = X0 + DX*FLOAT(I)
          YPLT = Y0 + DY*FLOAT(J)
          CALL PLSYMB(XPLT,YPLT,CH/96.0,5,0.0,0)
        ENDDO
      ENDDO
C
      CALL PLFLUSH
C
      RETURN
      END

      SUBROUTINE PDONE
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
      COMMON /DONE/ XDONE(2), YDONE(2)
C
C---- set DONE window
      XDONE(1) = CH + CH/4.0 + 0.125/SIZE
      XDONE(2) = CH + CH/4.0 + 1.025/SIZE
      YDONE(1) =    - CH/4.0 + 0.525/SIZE
      YDONE(2) =    - CH/4.0 + 1.425/SIZE
C
C---- plot DONE window
      CALL GETCOLOR(ICOL0)
      CALL NEWCOLORNAME('green')
      CALL PLOT(XDONE(1),YDONE(1),3)
      CALL PLOT(XDONE(2),YDONE(1),2)
      CALL PLOT(XDONE(2),YDONE(2),2)
      CALL PLOT(XDONE(1),YDONE(2),2)
      CALL PLOT(XDONE(1),YDONE(1),2)
C
      CHA = MIN( (XDONE(2)-XDONE(1))/8.0 , (YDONE(2)-YDONE(1))/1.5 )
      XCA = 0.5*(XDONE(2)+XDONE(1)) - 2.0*CHA
      YCA = 0.5*(YDONE(2)+YDONE(1)) - 0.5*CHA
      CALL PLCHAR(XCA,YCA,CHA,'DONE',0.0,4)
      CALL NEWCOLOR(ICOL0)
C
      RETURN
      END

      FUNCTION LDONE(XC,YC)
      COMMON /DONE/ XDONE(2), YDONE(2)
      LOGICAL LDONE
C
C---- return T if location XC,YC falls within DONE window
C
      LDONE = XC .GE. XDONE(1) .AND.
     &        XC .LE. XDONE(2) .AND.
     &        YC .GE. YDONE(1) .AND.
     &        YC .LE. YDONE(2)
C
      RETURN
      END

      SUBROUTINE PABORT
      COMMON /PLTC/ IDEV, IPSLU, SIZE, CH
      COMMON /ABRT/ XABORT(2), YABORT(2)
C
C---- set abort window
      XABORT(1) =  CH + CH/4.0 + 0.125/SIZE
      XABORT(2) =  CH + CH/4.0 + 1.025/SIZE
      YABORT(1) =     - CH/4.0
      YABORT(2) =     - CH/4.0 + 0.400/SIZE
C
C---- plot abort window
      CALL GETCOLOR(ICOL0)
      CALL NEWCOLORNAME('red')
      CALL PLOT(XABORT(1),YABORT(1),3)
      CALL PLOT(XABORT(2),YABORT(1),2)
      CALL PLOT(XABORT(2),YABORT(2),2)
      CALL PLOT(XABORT(1),YABORT(2),2)
      CALL PLOT(XABORT(1),YABORT(1),2)
C
      CHA = MIN( (XABORT(2)-XABORT(1))/8.0 , (YABORT(2)-YABORT(1))/1.5 )
      XCA = 0.5*(XABORT(2)+XABORT(1)) - 2.5*CHA
      YCA = 0.5*(YABORT(2)+YABORT(1)) - 0.5*CHA
      CALL PLCHAR(XCA,YCA,CHA,'ABORT',0.0,5)
      CALL NEWCOLOR(ICOL0)
C
      RETURN
      END

      FUNCTION LABORT(XC,YC)
      COMMON /ABRT/ XABORT(2), YABORT(2)
      LOGICAL LABORT
C
C---- return T if location XC,YC falls within abort window
C
      LABORT = XC .GE. XABORT(1) .AND.
     &         XC .LE. XABORT(2) .AND.
     &         YC .GE. YABORT(1) .AND.
     &         YC .LE. YABORT(2)
C
      RETURN
      END