text:
C$Procedure ZZDDHPPF ( Private --- DDH Prepare Preexisting File )

      SUBROUTINE ZZDDHPPF ( UNIT, ARCH, BFF )

C$ Abstract
C
C     SPICE Private routine intended solely for the support of SPICE
C     routines.  Users should not call this routine directly due
C     to the volatile nature of this routine.
C
C     Prepare preexisting binary file for entry into the handle
C     table.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     None.
C
C$ Keywords
C
C     PRIVATE
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE              'zzddhman.inc'
      INCLUDE              'zzftprms.inc'

      INTEGER               UNIT
      INTEGER               ARCH
      INTEGER               BFF

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     UNIT       I   Logical unit attached to the binary file.
C     ARCH       I   Integer code indicating the file architecture.
C     BFF        O   Integer code indicating the binary file format.
C
C$ Detailed_Input
C
C     UNIT       is a logical unit attached to the binary file to be
C                prepared for inclusion into the handle table.
C
C     ARCH       is an integer that indicates the architecture of
C                the file attached to UNIT.  Acceptable values are
C                the parameters:
C
C                   DAF
C                   DAS
C
C                defined in ZZDDHMAN.INC.
C
C$ Detailed_Output
C
C     BFF        is an integer that indicates the binary file format
C                of the DAF attached to UNIT.  Possible values are
C                the parameters:
C
C                   BIGI3E
C                   LTLI3E
C                   VAXGFL
C                   VAXDFL
C
C                defined in ZZDDHMAN.INC.
C
C$ Parameters
C
C     See the include file ZZDDHMAN.INC.
C
C$ Exceptions
C
C     1) SPICE(UNKNOWNFILARC) is signaled when ARCH is not in the
C        range of codes for known file architectures or the binary
C        file's ID word is unknown to IDW2AT.  BFF is set to 0 when
C        this error is signaled.  UNIT is not closed.
C
C     2) SPICE(FILEREADFAILED) is signaled when either of the two
C        READ statements in the module returns non-zero IOSTAT, thus
C        indicating read failure.  BFF is set to 0 in this case.  Unit
C        is not closed.
C
C     3) SPICE(FILARCMISMATCH) is signaled when the file attached to
C        UNIT is determined to utilize an architecture that is
C        different from the one to which the input argument ARCH
C        refers.  Unit is not closed.
C
C     4) SPICE(UNKNOWNBFF) is signaled whenever the binary file
C        format detection algorithm reaches a state of uncertainty
C        for DAFs.  This can be the result of several conditions,
C        an empty pre-N0052 DAF, reading a DAF with an unknown BFF
C        from a future toolkit, etc.  In all cases, BFF is set to 0.
C        Unit is not closed.
C
C     5) If a pre-FTP string binary is loaded, no FTP based
C        diagnostics are performed, and the file is assumed to be
C        in proper, working order.
C
C$ Files
C
C     This routine reads at least one, and potentially, several records
C     from the file attached to UNIT.
C
C$ Particulars
C
C     This routine exists to prepare a binary file for inclusion
C     in the handle table in ZZDDHMAN.  This includes verifying
C     that the file is suitable to load and determining the binary
C     file format where possible.
C
C     For DAF files:
C
C        The binary file format of old (pre-N0050) binaries is
C        detectable if the file is non-empty and undamaged.
C        New files contain the binary file format identification
C        string in the file record along with the FTP error
C        detection string.  They are correctly identified in most
C        cases, including damaged.
C
C     For DAS files:
C
C        The binary file format of old (pre-N0052) binaries is
C        not detectable.  This this module will assume that any
C        old DAS binaries are of the native format.  New binaries
C        include the binary file format identification string as
C        well as the FTP error detection string.  They are
C        correctly identified in most cases as well.
C
C     FTP Error Detection:
C
C        FTP error detection occurs when at least part of the
C        detection string is detected in the file record.  When
C        absent, no errors are signaled and the file is then
C        assumed to be an old binary.  In the event that the FTP
C        detection string is present, and additional unknown
C        sequences are present, diagnostics are only performed on
C        sequences known to this version of the toolkit.
C
C$ Examples
C
C     See ZZDDHMAN for sample usage.
C
C$ Restrictions
C
C     1) The file attached to UNIT was written on a platform whose
C        characters are of a single byte in length.
C
C     2) Numeric data when read as characters from the UNIT
C        preserves the bit patterns present in the file in
C        memory.
C
C     3) The intrinsic ICHAR preserves the bit pattern of the
C        character byte read from a file.  Namely if one examines
C        the integer created the 8 least significant bits will be
C        precisely those found in the character.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     F.S. Turner     (JPL)
C
C$ Version
C
C-    SPICELIB Version 2.25.0, 10-MAR-2014 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-INTEL.
C
C-    SPICELIB Version 2.24.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-LINUX-64BIT-IFORT.
C
C-    SPICELIB Version 2.23.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-GFORTRAN.
C
C-    SPICELIB Version 2.22.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-64BIT-GFORTRAN.
C
C-    SPICELIB Version 2.21.0, 10-MAR-2014 (BVS)
C
C        Updated for PC-CYGWIN-64BIT-GCC_C.
C
C-    SPICELIB Version 2.20.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL.
C
C-    SPICELIB Version 2.19.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL-CC_C.
C
C-    SPICELIB Version 2.18.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-INTEL-64BIT-CC_C.
C
C-    SPICELIB Version 2.17.0, 13-MAY-2010 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-NATIVE_C.
C
C-    SPICELIB Version 2.16.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-WINDOWS-64BIT-IFORT.
C
C-    SPICELIB Version 2.15.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-LINUX-64BIT-GFORTRAN.
C
C-    SPICELIB Version 2.14.0, 13-MAY-2010 (BVS)
C
C        Updated for PC-64BIT-MS_C.
C
C-    SPICELIB Version 2.13.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-INTEL_C.
C
C-    SPICELIB Version 2.12.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-IFORT.
C
C-    SPICELIB Version 2.11.0, 13-MAY-2010 (BVS)
C
C        Updated for MAC-OSX-64BIT-GFORTRAN.
C
C-    SPICELIB Version 2.10.0, 18-MAR-2009 (BVS)
C
C        Updated for PC-LINUX-GFORTRAN.
C
C-    SPICELIB Version 2.9.0, 18-MAR-2009 (BVS)
C
C        Updated for MAC-OSX-GFORTRAN.
C
C-    SPICELIB Version 2.8.0, 19-FEB-2008 (BVS)
C
C        Updated for PC-LINUX-IFORT.
C
C-    SPICELIB Version 2.7.0, 14-NOV-2006 (BVS)
C
C        Updated for PC-LINUX-64BIT-GCC_C.
C
C-    SPICELIB Version 2.6.0, 14-NOV-2006 (BVS)
C
C        Updated for MAC-OSX-INTEL_C.
C
C-    SPICELIB Version 2.5.0, 14-NOV-2006 (BVS)
C
C        Updated for MAC-OSX-IFORT.
C
C-    SPICELIB Version 2.4.0, 14-NOV-2006 (BVS)
C
C        Updated for PC-WINDOWS-IFORT.
C
C-    SPICELIB Version 2.3.0, 26-OCT-2005 (BVS)
C
C        Updated for SUN-SOLARIS-64BIT-GCC_C.
C
C-    SPICELIB Version 2.2.0, 03-JAN-2005 (BVS)
C
C        Updated for PC-CYGWIN_C.
C
C-    SPICELIB Version 2.1.0, 03-JAN-2005 (BVS)
C
C        Updated for PC-CYGWIN.
C
C-    SPICELIB Version 2.0.1, 17-JUL-2002 (BVS)
C
C        Added MAC-OSX environments.
C
C-    SPICELIB Version 2.0.0, 06-FEB-2002 (FST)
C
C        This routine was updated to load binaries created by
C        N0051 versions of Sun Solaris Native C Toolkits.  See
C        the Revisions section for details.
C
C-    SPICELIB Version 1.0.0, 04-OCT-2001 (FST)
C
C
C-&

C$ Revisions
C
C-    SPICELIB Version 2.0.0, 06-FEB-2002 (FST)
C
C        Shortly after releasing N0052, a few of our users
C        discovered that they were unable to load binary
C        DAFs created with the N0051 Sun Solaris Native C
C        (SUN-SOLARIS-NATIVE_C) Toolkits.  The reason for this
C        is the previous version of ZZDDHPPF released with N0052
C        assumed that if a DAF file record possessed a valid
C        FTP error detection string, then it must contain a
C        binary file format ID string as well.  Both were added
C        to the DAF file record in N0050.
C
C        However, a bug in the N0051 version of the ZZPLATFM
C        master file, the source of the binary file format ID
C        string for a given platform, neglected to assign a
C        value to the string.  Since it was a C environment
C        and the implementation of ZZPLATFM resulted in the
C        string being a static variable, it was initialized
C        to nulls and written into the file.
C
C        This version of ZZDDHPPF has been extended to recognize
C        a null binary file format ID, and apply the byte
C        examination algorithm used on pre-N0050 DAFs to determine
C        its format.
C
C-&

C
C     SPICELIB Functions
C
      INTEGER               ISRCHC
      INTEGER               POS

      LOGICAL               RETURN

C
C     Local Parameters
C
C     Number of characters to be read in from a record.
C
      INTEGER               NUMCHR
      PARAMETER           ( NUMCHR = 1000 )

C
C     Bounding indices for the window that brackets the FTP
C     error detection string in the file record.
C
      INTEGER               FTPSTR
      PARAMETER           ( FTPSTR = 500 )

      INTEGER               FTPEND
      PARAMETER           ( FTPEND = 1000 )

C
C     Index of the start of the binary file format identification
C     string in DAF binaries.
C
      INTEGER               DAFBFF
      PARAMETER           ( DAFBFF = 89 )

C
C     Index of the start of the binary file format identification
C     string in DAS binaries.
C
      INTEGER               DASBFF
      PARAMETER           ( DASBFF = 85 )

C
C     Size of the binary format identification string.
C
      INTEGER               BFFSIZ
      PARAMETER           ( BFFSIZ = 8 )

C
C     Index of the first byte of NI in the DAF file record.
C
      INTEGER               NIPOS
      PARAMETER           ( NIPOS = 13 )

C
C     Index of the first byte of NSUM in the DAF descriptor record.
C
      INTEGER               NSMPOS
      PARAMETER           ( NSMPOS = 17 )

C
C     Index of the first byte of FDREC in the DAF file record.
C
      INTEGER               FDRPOS
      PARAMETER           ( FDRPOS = 77 )

C
C     Integer code such that CHAR(INTNUL) produces the NULL character.
C
      INTEGER               INTNUL
      PARAMETER           ( INTNUL = 0 )

C
C     IDW2AT Output Argument Lengths.
C
      INTEGER               IDWLEN
      PARAMETER           ( IDWLEN = 4 )

C
C     NULLID is the index of the extended STRBFF "NULL" string ID.
C
      INTEGER               NULLID
      PARAMETER           ( NULLID = NUMBFF+1 )

C
C     Local Variables
C
      CHARACTER*(1)         CARG
      CHARACTER*(NUMCHR)    CHRREC
      CHARACTER*(BFFSIZ)    BFFIDW
      CHARACTER*(IDWLEN)    FILARC
      CHARACTER*(IDWLEN)    FILTYP
      CHARACTER*(SIZDLM)    FTPDLM
      CHARACTER*(SIZEND)    FTPLFT
      CHARACTER*(SIZSTR)    FTPMEM
      CHARACTER*(SIZEND)    FTPRGT
      CHARACTER*(1)         NULL
      CHARACTER*(STRSIZ)    STRARC ( NUMARC )
      CHARACTER*(STRSIZ)    STRBFF ( NULLID )

      INTEGER               FDREC
      INTEGER               FTPPOS
      INTEGER               I
      INTEGER               IOSTAT
      INTEGER               TSTARC

      LOGICAL               FIRST
      LOGICAL               FOUND
      LOGICAL               FTPERR

C
C     Statement Functions
C
      INTEGER               ZZICHR

C
C     Saved Variables
C
      SAVE                  FIRST
      SAVE                  FTPDLM
      SAVE                  FTPLFT
      SAVE                  FTPMEM
      SAVE                  FTPRGT
      SAVE                  NULL
      SAVE                  STRARC
      SAVE                  STRBFF

C
C     Data Statements
C
      DATA                  FIRST  / .TRUE. /

C
C     Statement Function Definitions
C
C     This function controls the conversion of characters to integers.
C     On some supported environments, ICHAR is not sufficient to
C     produce the desired results.  This however, is not the case
C     with this particular environment.
C
      ZZICHR(CARG) = ICHAR(CARG)

C
C     Standard SPICE error handling.
C
      IF ( RETURN() ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'ZZDDHPPF' )
      END IF

C
C     If this is the first time into the routine, populate local
C     copies of reference values.  This includes the names of the
C     BFF parameters, the names of the ARCH parameters, and the
C     local copy of the FTP string.
C
      IF ( FIRST ) THEN

C
C        Construct and store the NULL valued byte.
C
         NULL = CHAR(INTNUL)

C
C        Retrieve the BFF and ARCH names.
C
         DO I = 1, NUMBFF
            CALL ZZDDHGSD ( 'BFF', I, STRBFF(I) )
         END DO

         DO I = 1, NUMARC
            CALL ZZDDHGSD ( 'ARCH', I, STRARC(I) )
         END DO

C
C        Extend STRBFF to include the null BFFID.  This addresses
C        the N0051 Sun Solaris Native C toolkit binary files.
C
         DO I = 1, BFFSIZ
            STRBFF(NULLID) (I:I) = NULL
         END DO

C
C        Fetch the FTP string.
C
         CALL ZZFTPSTR ( FTPMEM, FTPLFT, FTPRGT, FTPDLM )

C
C        Set FIRST to FALSE so we will not reassign any of these values.
C
         FIRST = .FALSE.

      END IF

C
C     Get the simple consistency checks out of the way first.  Is
C     the input ARCH value valid?
C
      IF ( ( ARCH .LE. 0 ) .OR. ( ARCH .GT. NUMARC ) ) THEN

         BFF = 0

         CALL SETMSG ( 'The integer code, ''#'' indicating the file '
     .   //            'architecture to examine is out of range.'     )
         CALL ERRINT ( '#', ARCH                                      )
         CALL SIGERR ( 'SPICE(UNKNOWNFILARC)'                         )
         CALL CHKOUT ( 'ZZDDHPPF'                                     )
         RETURN

      END IF

C
C     Read the first record from the file as a string of NUMCHR
C     characters.
C
      READ ( UNIT, REC=1, IOSTAT=IOSTAT ) CHRREC

C
C     Check for read failure.
C
      IF ( IOSTAT .NE. 0 ) THEN

         BFF = 0

         CALL SETMSG ( 'Error reading the file record from the '
     .   //            'binary DAF file ''#''.  IOSTAT = #.'     )
         CALL ERRFNM ( '#', UNIT                                 )
         CALL ERRINT ( '#', IOSTAT                               )
         CALL SIGERR ( 'SPICE(FILEREADFAILED)'                   )
         CALL CHKOUT ( 'ZZDDHPPF'                                )
         RETURN

      END IF

C
C     First check the ID word from the input file.
C
      CALL IDW2AT ( CHRREC(1:8), FILARC, FILTYP )

C
C     Now locate FILARC in the STRARC array.
C
      TSTARC = ISRCHC ( FILARC, NUMARC, STRARC )

C
C     If FILARC was not found, signal an appropriate error.
C
      IF ( TSTARC .EQ. 0 ) THEN

         BFF = 0

         CALL SETMSG ( 'The file, #, has a unidentified file '
     .   //            'architecture.  Check that this file '
     .   //            'is a properly created binary SPICE '
     .   //            'kernel.'                               )
         CALL ERRFNM ( '#', UNIT                               )
         CALL SIGERR ( 'SPICE(UNKNOWNFILARC)'                  )
         CALL CHKOUT ( 'ZZDDHPPF'                              )
         RETURN

C
C     Otherwise we have an architecture mismatch error, if
C     FILARC does not agree with ARCH.
C
      ELSE IF ( TSTARC .NE. ARCH ) THEN

         BFF = 0

         CALL SETMSG ( 'A request to load the # file, $, has '
     .   //            'been made by the % system.  This operation '
     .   //            'is not permitted.'                           )
         CALL ERRCH  ( '#', STRARC(TSTARC)                           )
         CALL ERRFNM ( '$', UNIT                                     )
         CALL ERRCH  ( '%', STRARC(ARCH)                             )
         CALL SIGERR ( 'SPICE(FILARCHMISMATCH)'                      )
         CALL CHKOUT ( 'ZZDDHPPF'                                    )
         RETURN

      END IF

C
C     Now check for possible FTP transfer errors.
C
      CALL ZZFTPCHK ( CHRREC(FTPSTR:FTPEND), FTPERR )

      IF ( FTPERR ) THEN

         BFF = 0

         CALL SETMSG ( 'FTP transfer error detected.  This '
     .   //            'binary $, ''#'', has most likely been '
     .   //            'corrupted by an ASCII mode FTP transfer. '
     .   //            'Obtain the file using IMAGE or BINARY '
     .   //            'transfer mode from the source.'            )
         CALL ERRCH  ( '$', STRARC(TSTARC)                         )
         CALL ERRFNM ( '#', UNIT                                   )
         CALL SIGERR ( 'SPICE(FTPXFERERROR)'                       )
         CALL CHKOUT ( 'ZZDDHPPF'                                  )
         RETURN

      END IF

C
C     Now this search is redundant, but the presence of the
C     FTPLFT string in the latter half of the file record
C     is fairly conclusive evidence that this is a "new" binary,
C     and we can expect to locate the binary file format
C     identification string.
C
      FTPPOS = POS ( CHRREC(500:1000), FTPLFT, 1 )

C
C     Check to see if we found FTPLFT.  If so extract the binary
C     file format ID word from the file record.
C
      IF ( FTPPOS .NE. 0 ) THEN

C
C        Extract BFFIDW from CHRREC.
C
         IF ( ARCH .EQ. DAF ) THEN
            BFFIDW = CHRREC ( DAFBFF : (DAFBFF+BFFSIZ-1) )
         ELSE IF ( ARCH .EQ. DAS ) THEN
            BFFIDW = CHRREC ( DASBFF : (DASBFF+BFFSIZ-1) )
         END IF

C
C        See if we can find BFFIDW in the STRBFF list.
C
         BFF = ISRCHC ( BFFIDW, NUMBFF + 1, STRBFF )

C
C        Check to see if BFF is 0, if it is, signal an error since
C        this indicates an unrecognized BFF.
C
         IF ( BFF .EQ. 0 ) THEN

            CALL SETMSG ( 'The file ''#'' utilizes the binary file '
     .      //            'format ''#''.  This format is currently '
     .      //            'unknown to this toolkit.  A toolkit '
     .      //            'update may be in order.'                  )
            CALL ERRFNM ( '#', UNIT                                  )
            CALL ERRCH  ( '#', BFFIDW                                )
            CALL SIGERR ( 'SPICE(UNKNOWNBFF)'                        )
            CALL CHKOUT ( 'ZZDDHPPF'                                 )
            RETURN

         END IF

C
C        See if we have a NULLID situation, if not check out and
C        return as swe have identified the BFF.
C
         IF ( BFF .NE. NULLID ) THEN
            CALL CHKOUT ( 'ZZDDHPPF' )
            RETURN
         END IF

      END IF

C
C     There is no FTP string, if the file is a DAS, we have to
C     assume it is of the native architecture.
C
      IF ( ARCH .EQ. DAS ) THEN

         CALL ZZPLATFM ( 'FILE_FORMAT', BFFIDW )
         CALL UCASE    ( BFFIDW,        BFFIDW )

         BFF = ISRCHC ( BFFIDW, NUMBFF, STRBFF )

         IF ( BFF .EQ. 0 ) THEN

            CALL SETMSG ( 'The native architecture for this '
     .      //            'platform is unknown to this version of '
     .      //            'the toolkit. This is a severe problem '
     .      //            'that should never occur, please contact '
     .      //            'NAIF.'                                    )
            CALL SIGERR ( 'SPICE(BUG)'                               )
            CALL CHKOUT ( 'ZZDDHPPF'                                 )
            RETURN

         END IF

         CALL CHKOUT ( 'ZZDDHPPF' )
         RETURN

      END IF

C
C     If we reach this point, then we are either dealing with
C     an old DAF (created by a pre-N0050 toolkit) or one of the
C     DAFs created by the N0051 Sun Solaris Native C version of
C     the toolkit.  This requires an examination of the bits
C     and bytes in the file that works this way:
C
C        Since in a valid DAF, 2 <= NI <= 250, we can easily
C        determine whether the 4 bytes used to store NI in the
C        file record are little or big endian.  If we discover
C        that the integer is encoded as big-endian, then stop
C        as this file must be 'BIG-IEEE'.  If it is little
C        endian, then locate the first descriptor record
C        in the file.
C
C        Read the first descriptor record.  Extract NSUM, the
C        3rd DP from the record.  If it is 0.0D0, signal an error
C        as this is an empty DAF and we can not determine its
C        type.  If it's non-zero, then check to see if the first
C        4 bytes are "0s".  If they are it must be 'LTL-IEEE'.
C        Otherwise pass it off to ZZDDHIVF to discriminate between
C        'VAX-GFLT' and 'VAX-DFLT'.  We know the first 4 bytes must
C        be "0s" in the 'LTL-IEEE" case, since NSUM is subject to
C        the following inequality: 1 <= NSUM <= 125
C
C     Having laid out the scheme, let's get to it.  First take a
C     look at the four character bytes that hold NI.  These bytes
C     be one of the following:
C
C        Little Endian:  VAL, 0, 0, 0
C           Big Endian:    0, 0, 0, VAL
C
C     where VAL is some non-zero value.
C
      IF ( ( CHRREC(     NIPOS:NIPOS     ) .EQ. NULL ) .AND.
     .     ( CHRREC( (NIPOS+1):(NIPOS+1) ) .EQ. NULL ) .AND.
     .     ( CHRREC( (NIPOS+2):(NIPOS+2) ) .EQ. NULL ) .AND.
     .     ( CHRREC( (NIPOS+3):(NIPOS+3) ) .NE. NULL )       ) THEN

         BFF = BIGI3E

      ELSE IF ( ( CHRREC(     NIPOS:NIPOS     ) .NE. NULL ) .AND.
     .          ( CHRREC( (NIPOS+1):(NIPOS+1) ) .EQ. NULL ) .AND.
     .          ( CHRREC( (NIPOS+2):(NIPOS+2) ) .EQ. NULL ) .AND.
     .          ( CHRREC( (NIPOS+3):(NIPOS+3) ) .EQ. NULL )       ) THEN

C
C        At this point we know we are dealing with a little endian
C        file.  Locate the first descriptor record.
C
         FDREC = ZZICHR( CHRREC(     FDRPOS:FDRPOS     ) )
         FDREC = ZZICHR( CHRREC( (FDRPOS+1):(FDRPOS+1) ) )*16   + FDREC
         FDREC = ZZICHR( CHRREC( (FDRPOS+2):(FDRPOS+2) ) )*256  + FDREC
         FDREC = ZZICHR( CHRREC( (FDRPOS+3):(FDRPOS+3) ) )*4096 + FDREC

C
C        Read the record into CHRREC.
C
         READ ( UNIT, REC=FDREC, IOSTAT=IOSTAT ) CHRREC

C
C        Check for read failure.
C
         IF ( IOSTAT .NE. 0 ) THEN

            BFF = 0

            CALL SETMSG ( 'Error reading a descriptor record from '
     .      //            'the binary DAF file ''#''.  IOSTAT = #.' )
            CALL ERRFNM ( '#', UNIT                                 )
            CALL ERRINT ( '#', IOSTAT                               )
            CALL SIGERR ( 'SPICE(FILEREADFAILED)'                   )
            CALL CHKOUT ( 'ZZDDHPPF'                                )
            RETURN

         END IF

C
C        Now examine the NSUM DP in this record to determine the
C        architecture.
C
         IF ( ( CHRREC(     NSMPOS:NSMPOS )     .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+1):(NSMPOS+1) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+2):(NSMPOS+2) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+3):(NSMPOS+3) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+4):(NSMPOS+4) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+5):(NSMPOS+5) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+6):(NSMPOS+6) ) .EQ. NULL ) .AND.
     .        ( CHRREC( (NSMPOS+7):(NSMPOS+7) ) .EQ. NULL )       ) THEN

C
C           In this case we have an empty DAF, and can not distinguish
C           between little endian formats.  Signal an error and return.
C
            BFF = 0

            CALL SETMSG ( 'The DAF, ''#'', appears to contain no '
     .      //            'data.  As such, its binary file format '
     .      //            'can not be determined which prevents '
     .      //            'it from being loaded.'                   )
            CALL ERRFNM ( '#', UNIT                                 )
            CALL SIGERR ( 'SPICE(UNKNOWNBFF)'                       )
            CALL CHKOUT ( 'ZZDDHPPF'                                )
            RETURN

         ELSE IF ( ( CHRREC(     NSMPOS:NSMPOS     ) .EQ. NULL ) .AND.
     .             ( CHRREC( (NSMPOS+1):(NSMPOS+1) ) .EQ. NULL ) .AND.
     .             ( CHRREC( (NSMPOS+2):(NSMPOS+2) ) .EQ. NULL ) .AND.
     .             ( CHRREC( (NSMPOS+3):(NSMPOS+3) ) .EQ. NULL )  ) THEN

C
C           In this case the file is little endian IEEE.  Set BFF.
C
            BFF = LTLI3E

         ELSE

C
C           We are probably looking at a VAX file.  Find out which
C           format.
C
            CALL ZZDDHIVF ( CHRREC(NSMPOS:(NSMPOS+7)), BFF, FOUND )

            IF ( .NOT. FOUND ) THEN

               BFF = 0

               CALL SETMSG ( 'Unable to determine the binary file '
     .         //            'format of DAF ''#''.'                 )
               CALL ERRFNM ( '#', UNIT                              )
               CALL SIGERR ( 'SPICE(UNKNOWNBFF)'                    )
               CALL CHKOUT ( 'ZZDDHPPF'                             )
               RETURN

            END IF

         END IF

      ELSE

         BFF = 0

      END IF

      CALL CHKOUT ( 'ZZDDHPPF' )
      RETURN

      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
!> \file npc_sw_proc.f90
!> Contains module npc_surfacewater_processes.

!>Nitrogen, phosphorus and organic carbon processes in surface water in HYPE
MODULE NPC_SURFACEWATER_PROCESSES

  !Copyright 2012-2017 SMHI
  !
  !This file is part of HYPE.
  !HYPE is free software: you can redistribute it and/or modify it under the terms of the Lesser GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
  !HYPE 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 Lesser GNU General Public License for more details.
  !You should have received a copy of the Lesser GNU General Public License along with HYPE. If not, see <http://www.gnu.org/licenses/>.

  !Used modules
  USE STATETYPE_MODULE, ONLY : riverstatetype,lakestatetype
  USE GENERAL_WATER_CONCENTRATION
  USE GENERAL_FUNCTIONS
  !Subroutines also uses modvar and hypevariables
  IMPLICIT NONE
  PRIVATE
  !----------------------------------------------
  ! Private procedures
  !----------------------------------------------
  ! denitrification_water
  ! production_mineralisation
  ! lake_sedimentation
  ! sedimentation_resuspension
  ! calculate_lake_tpmean
  ! calculate_river_tpmean
  ! internal_lake_load
  ! oc_production_mineralisation
  ! calculate_wetland_np
  !----------------------------------------------
  PUBLIC :: initiate_river_npc_state, &
       initiate_lake_npc_state, &
       add_dry_deposition_to_lake, &
       add_dry_deposition_to_river, &
       np_processes_in_river, &
       np_processes_in_lake, &
       oc_processes_in_river, &
       oc_processes_in_lake, &
       add_diffuse_source_to_local_river, &
       add_point_sources_to_main_river, &
       calculate_river_wetland, &
       set_lake_slowwater_maxvolume
CONTAINS

  !>\brief Initiation river variables for nutrients and organic
  !>carbon simulations. Concentration (mg/L)
  !>
  !>\b Consequences Module hypevariables variable Qmax, Q2max,
  !>iQmax, iQ2max may be allocated and set.
  !-----------------------------------------------------------------
  SUBROUTINE initiate_river_npc_state(initN,initP,initC,initS,initT1,riverstate)

    USE HYPEVARIABLES, ONLY : Qmax,       &   !OUT
                              Q2max,      &   !OUT
                              iQmax,      &   !OUT
                              iQ2max,     &   !OUT
                              m_iniT2,    &
                              m_tpmean,   &
                              m_ldTPmean
    USE MODVAR,      ONLY : nsub, &
                            basin, &
                            regpar, &
                            regiondivision, &
                            lakedatapar, &
                            lakedataparindex

    !Argument declarations
    LOGICAL, INTENT(IN) :: initN  !<flag for initiation of nitrogen model
    LOGICAL, INTENT(IN) :: initP  !<flag for initiation of phosphorus model
    LOGICAL, INTENT(IN) :: initC  !<flag for initiation of organic carbon model
    LOGICAL, INTENT(IN) :: initS  !<flag for simulation of sediment model
    LOGICAL, INTENT(IN) :: initT1 !<flag for initiation of tracer model T1
    TYPE(riverstatetype),INTENT(INOUT)   :: riverstate  !<River states

    !Local variables
    INTEGER isb

    !>\b Algorithm \n
    !>Allocate and initialize river sediment variables
    IF(initN.OR.initP.OR.initT1)THEN !For calculation of bankful flow and Qmean
      IF(.NOT.ALLOCATED(Qmax))   ALLOCATE(Qmax(2,nsub))
      IF(.NOT.ALLOCATED(Q2max))  ALLOCATE(Q2max(2,nsub))
      IF(.NOT.ALLOCATED(iQmax))  ALLOCATE(iQmax(2,nsub))
      IF(.NOT.ALLOCATED(iQ2max)) ALLOCATE(iQ2max(2,nsub))
      riverstate%Q365 = 0.0001
      riverstate%Qdayacc = 0.000
      Qmax = 0.0001; Q2max = 0.0001
      iQmax = 365; iQ2max = 364
    ENDIF

    !>Set TPmean-variable if phosphorus is not calculated by HYPE
    IF((initN.OR.initC.OR.initS).AND.(.NOT.initP))THEN
      DO isb = 1,nsub
        IF(basin(isb)%parregion(regiondivision(m_tpmean))>0) riverstate%TPmean(:,isb) = regpar(m_tpmean,basin(isb)%parregion(regiondivision(m_tpmean)))
        IF(ALLOCATED(lakedatapar)) riverstate%TPmean(1,isb) = lakedatapar(lakedataparindex(isb,1),m_ldtpmean)
        IF(ALLOCATED(lakedatapar)) riverstate%TPmean(2,isb) = lakedatapar(lakedataparindex(isb,2),m_ldtpmean)
      ENDDO
    ENDIF

  END SUBROUTINE initiate_river_npc_state

  !>Initiation lake for nutrients and organic carbon. Concentration in
  !(mg/L)
  !!
  !>\b Consequences Module hypevariables variable slowlakeini
  !> may be allocated and set.
  !-----------------------------------------------------------------
  SUBROUTINE initiate_lake_npc_state(initN,initP,initC,initS,initT,initilake,initolake,lakestate)

    USE HYPEVARIABLES, ONLY : slowlakeini,  &   !OUT
                              m_tpmean,     &
                              m_tnmean,     &
                              m_tocmean,    &
                              m_ldtpmean,   &
                              m_ldtnmean,   &
                              m_ldtocmean,  &
                              m_lddeeplake, &
                              m_iniT2
    USE MODVAR, ONLY : basin, &
                       nsub, &
                       genpar, &
                       regpar, &
                       regiondivision, &
                       lakedatapar, &
                       lakedataparindex, &
                       i_in,i_on, &
                       i_pp,i_oc,i_t2

    !Argument declarations
    LOGICAL, INTENT(IN) :: initN     !<flag for initiation of nitrogen model
    LOGICAL, INTENT(IN) :: initP     !<flag for initiation of phosphorus model
    LOGICAL, INTENT(IN) :: initC     !<flag for initiation of organic carbon model
    LOGICAL, INTENT(IN) :: initS     !<flag for simulation of sediment model
    LOGICAL, INTENT(IN) :: initT     !<flag for initiation of tracer model (T1 or T2)
    LOGICAL, INTENT(IN) :: initilake !<flag for ilake class existance
    LOGICAL, INTENT(IN) :: initolake !<flag for olake class existance
    TYPE(lakestatetype),INTENT(INOUT)    :: lakestate   !<Lake states

    !Local variables
    INTEGER j,isb     !loop-variables (class,subbasin)
    REAL, ALLOCATABLE :: inidepth(:,:)

    !Allocate variables for lake partitioning in fastflow (lakewi) and slow lake parts.
    IF(initN.OR.initP.OR.initC.OR.initT)THEN
      CALL set_lake_slowwater_maxvolume(nsub,basin%lakedepth(1),basin(:)%lakedepth(2), &
                   lakedataparindex,lakedatapar(:,m_lddeeplake),initilake,initolake)
      ALLOCATE(inidepth(2,nsub))

      !Initiate default lake volumes (mm)
      inidepth = lakestate%water     ! water is earlier initiated as whole lake wolume
      !Redistribute initial lake volume by lake partitioning
      IF(initilake)THEN
        DO isb = 1,nsub
          IF(inidepth(1,isb)>slowlakeini(1,isb))THEN
            lakestate%slowwater(1,isb) = slowlakeini(1,isb)
            lakestate%water(1,isb) = inidepth(1,isb) - lakestate%slowwater(1,isb)
          ELSEIF(inidepth(1,isb)<=slowlakeini(1,isb))THEN
            lakestate%slowwater(1,isb) = inidepth(1,isb)
            lakestate%water(1,isb) = 0
          ENDIF
        ENDDO
      ENDIF
      IF(initolake)THEN
        DO isb = 1,nsub
          IF(inidepth(2,isb)>slowlakeini(2,isb))THEN
            lakestate%slowwater(2,isb) = slowlakeini(2,isb)
            lakestate%water(2,isb) = inidepth(2,isb) - lakestate%slowwater(2,isb)
          ELSEIF(inidepth(2,isb)<=slowlakeini(2,isb))THEN
            lakestate%slowwater(2,isb) = inidepth(2,isb)
            lakestate%water(2,isb) = 0
          ENDIF
        ENDDO
      ENDIF
      DEALLOCATE(inidepth)
    ENDIF

    !Initialize lake concentration (mg/L)
    IF(initN.OR.initP)THEN
      DO isb = 1,nsub
        IF(initP)THEN
          IF(ALLOCATED(lakedatapar))THEN
            lakestate%conc(i_pp,:,isb) = lakedatapar(lakedataparindex(isb,:),m_ldtpmean)
          ELSEIF(basin(isb)%parregion(regiondivision(m_tpmean))>0)THEN
            lakestate%conc(i_pp,:,isb) = regpar(m_tpmean,basin(isb)%parregion(regiondivision(m_tpmean)))
          ENDIF
        ENDIF
        IF(initN)THEN
          IF(ALLOCATED(lakedatapar))THEN
            lakestate%conc(i_on,:,isb) = lakedatapar(lakedataparindex(isb,:),m_ldtnmean)*0.5
            lakestate%conc(i_in,:,isb) = lakedatapar(lakedataparindex(isb,:),m_ldtnmean)*0.5
          ELSEIF(basin(isb)%parregion(regiondivision(m_tnmean))>0)THEN
            lakestate%conc(i_on,:,isb) = regpar(m_tnmean,basin(isb)%parregion(regiondivision(m_tnmean)))*0.5
            lakestate%conc(i_in,:,isb) = regpar(m_tnmean,basin(isb)%parregion(regiondivision(m_tnmean)))*0.5
          ENDIF
        ENDIF
      ENDDO
    ENDIF

    IF(initC)THEN
      DO isb = 1,nsub
        IF(ALLOCATED(lakedatapar))THEN
          lakestate%conc(i_oc,:,isb) = lakedatapar(lakedataparindex(isb,:),m_ldtocmean)
        ELSEIF(basin(isb)%parregion(regiondivision(m_tocmean))>0)THEN
          lakestate%conc(i_oc,:,isb) = regpar(m_tocmean,basin(isb)%parregion(regiondivision(m_tocmean)))
        ENDIF
      ENDDO
    ENDIF

    IF(i_t2>0)THEN
      DO isb = 1,nsub
        lakestate%conc(i_t2,:,isb) = genpar(m_iniT2)
      ENDDO
      lakestate%lowertemp(:,:) = genpar(m_iniT2)
      lakestate%uppertemp(:,:) = genpar(m_iniT2)
    ENDIF

    IF(initN.OR.initP.OR.initC.OR.initT)THEN
      DO isb = 1,nsub
        DO j = 1,2
          lakestate%concslow(:,j,isb) = lakestate%conc(:,j,isb)
        ENDDO
      ENDDO
    ENDIF

    !Set TPmean-variable if phosphorus is not calculated by HYPE
    IF((initN.OR.initC.OR.initS).AND..NOT.initP)THEN
      DO isb = 1,nsub
        IF(ALLOCATED(lakedatapar))THEN
          lakestate%TPmean(1,isb) = lakedatapar(lakedataparindex(isb,1),m_ldtpmean)
          lakestate%TPmean(2,isb) = lakedatapar(lakedataparindex(isb,2),m_ldtpmean)
        ELSEIF(basin(isb)%parregion(regiondivision(m_tpmean))>0)THEN
          lakestate%TPmean(:,isb) = regpar(m_tpmean,basin(isb)%parregion(regiondivision(m_tpmean)))
        ENDIF
      ENDDO
    ENDIF

  END SUBROUTINE initiate_lake_npc_state

  !>Calculate the maximum slowwater volume of lakes. Used for simulation of nutrients and organic carbon.
  !!The variable is used for lake partitioning in fastflow and slow lake parts.
  !>
  !>\b Consequences Module hypevariables variable slowlakeini may be set.
  !>
  !>\b Reference ModelDescription Chapter Rivers and lakes (Basic assumptions)
  !-----------------------------------------------------------------
  SUBROUTINE set_lake_slowwater_maxvolume(n,ildepth,oldepth,ldparindex,ldpar,initilake,initolake)

    USE HYPEVARIABLES, ONLY : slowlakeini  !OUT

    !Argument declarations
    INTEGER, INTENT(IN) :: n          !<number of subbasins
    REAL, INTENT(IN)    :: ildepth(n) !<ilake depth
    REAL, INTENT(IN)    :: oldepth(n) !<olake depth
    INTEGER, INTENT(IN) :: ldparindex(n,2)  !<lakedataparindex
    REAL, INTENT(IN)    :: ldpar(:)   !<deeplake parameter (lakedatapar(:,m_lddeeplake))
    LOGICAL, INTENT(IN) :: initilake  !<flag for ilake class existance
    LOGICAL, INTENT(IN) :: initolake  !<flag for olake class existance

    !Local variables
    INTEGER isb         !loop-variable
    REAL    deeplake    !model parameter deeplake
    REAL    inidepth(n)

    !Allocate variable for maximum slow lake part
    IF(.NOT. ALLOCATED(slowlakeini))  ALLOCATE(slowlakeini(2,n))
    slowlakeini = 0.    !mm
    !Set variable for maximum slow lake part for ilake to ilake depth or adjusted by parameter deeplake.
    IF(initilake)THEN
      inidepth = ildepth*1000.
      slowlakeini(1,:) = inidepth   !default
      DO isb = 1,n
        deeplake = ldpar(ldparindex(isb,1))
        IF(deeplake>0) slowlakeini(1,isb) = deeplake * inidepth(isb)           !ilake: slowlake target water stage (mm)
      ENDDO
    ENDIF
    !Set variable for maximum slow lake part for olake to olake depth or adjusted by parameter deeplake.
    IF(initolake)THEN
      inidepth = oldepth * 1000.
      slowlakeini(2,:) = inidepth   !default
      DO isb = 1,n
        deeplake = ldpar(ldparindex(isb,2))
        IF(deeplake>0) slowlakeini(2,isb) = deeplake * inidepth(isb)           !olake: slowlake target water stage (mm)
      ENDDO
    ENDIF

  END SUBROUTINE set_lake_slowwater_maxvolume

  !>\brief Calculate atmospheric dry deposition of N and P and add it
  !>to lakewater
  !>
  !>\b Reference ModelDescription Chapter Processes above ground (Atmospheric deposition of nitrogen and phosphorus)
  !----------------------------------------------------------------------------
  SUBROUTINE add_dry_deposition_to_lake(i,areaij,pooltype,source,veg,dryin,drypp,lakestate)

    USE MODVAR, ONLY : numsubstances,     &
                       i_in,i_pp

    !Argument declarations
    INTEGER, INTENT(IN) :: i                          !<index of subbasin
    REAL, INTENT(IN)    :: areaij                     !<classarea (km2)
    INTEGER, INTENT(IN) :: pooltype                   !<laketype: 1=ilake, 2=olake
    REAL, INTENT(OUT)   :: source(numsubstances)      !<dry deposition (kg/timestep)
    INTEGER, INTENT(IN) :: veg                        !<vegetation index
    REAL, INTENT(IN)    :: dryin(3)                   !<dry deposition IN (kg/km2/timestep)
    REAL, INTENT(IN)    :: drypp                      !<dry deposition PP (kg/km2/timestep)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate    !<Lake state

    !Local variables
    REAL :: sourcedd(numsubstances)

    source = 0.
    IF(i_in==0.AND.i_pp==0)RETURN

    !Prepare dry deposition
    sourcedd = 0.
    IF(i_in>0 .AND. veg>0) sourcedd(i_in) = dryin(veg)
    IF(i_pp>0) sourcedd(i_pp) = drypp

    !Add dry deposition of Inorg-N and PartP to lake
    IF(lakestate%water(pooltype,i)>0)THEN
      CALL add_source_to_water(lakestate%water(pooltype,i),numsubstances,lakestate%conc(:,pooltype,i),sourcedd)
    ELSE
      CALL add_source_to_water(lakestate%slowwater(pooltype,i),numsubstances,lakestate%concslow(:,pooltype,i),sourcedd)
    ENDIF

    !Calculate atmospheric dry deposition loads (kg/timestep)
    source = sourcedd*areaij

    RETURN

  END SUBROUTINE add_dry_deposition_to_lake

  !>\brief Calculate atmospheric dry deposition of N and P and add it
  !>to riverwater
  !>
  !\b Reference ModelDescription Chapter Processes above ground (Atmospheric
  !>deposition of nitrogen and phosphorus)
  !----------------------------------------------------------------------------
  SUBROUTINE add_dry_deposition_to_river(i,areaij,pooltype,source,veg,dryin,drypp,riverstate)

    USE MODVAR, ONLY : numsubstances,     &
                       i_in,i_pp
    USE HYPEVARIABLES, ONLY : ttpart,ttstep

    !Argument declarations
    INTEGER, INTENT(IN) :: i                          !<index of subbasin
    REAL, INTENT(IN)    :: areaij                     !<classarea (km2)
    INTEGER, INTENT(IN) :: pooltype                   !<rivertype: 1=lriver, 2=mriver
    REAL, INTENT(OUT)   :: source(numsubstances)      !< dry deposition (kg/timestep)
    INTEGER, INTENT(IN) :: veg                        !<vegetation index
    REAL, INTENT(IN)    :: dryin(3)                   !<dry deposition IN (kg/km2/timestep)
    REAL, INTENT(IN)    :: drypp                      !<dry deposition PP (kg/km2/timestep)
    TYPE(riverstatetype),INTENT(INOUT) :: riverstate  !<River state

    !Local variables
    INTEGER l
    REAL sourcedd(numsubstances)
    REAL totvol
    REAL waterfrac

    source = 0.
    IF(i_in==0.AND.i_pp==0)RETURN

    !Prepare dry deposition
    IF(i_in>0 .AND. veg>0) source(i_in) = dryin(veg)*areaij   !kg/timestep
    IF(i_pp>0) source(i_pp) = drypp*areaij
    sourcedd = source*1000.   !g=mg/L*m3

    IF(SUM(sourcedd)>0.)THEN
      !Fractions of river water i different components
      totvol = riverstate%water(pooltype,i) + (SUM(riverstate%qqueue(1:ttstep(pooltype,i),pooltype,i)) + riverstate%qqueue(ttstep(pooltype,i)+1,pooltype,i) * ttpart(pooltype,i))
      IF(totvol<=0) RETURN
      waterfrac = riverstate%water(pooltype,i)/totvol

      !Add dry deposition of Inorg-N and PartP to river watercourse
      IF(riverstate%water(pooltype,i)>0)THEN
        CALL add_source_to_water(riverstate%water(pooltype,i),numsubstances,riverstate%conc(:,pooltype,i),waterfrac*sourcedd)
      ENDIF
      DO l = 1,ttstep(pooltype,i)
        IF(riverstate%qqueue(l,pooltype,i)>0)THEN
          waterfrac = riverstate%qqueue(l,pooltype,i)/totvol
          CALL add_source_to_water(riverstate%qqueue(l,pooltype,i),numsubstances,riverstate%cqueue(:,l,pooltype,i),waterfrac*sourcedd)
        ENDIF
      ENDDO
      IF(ttpart(pooltype,i)>0)THEN
        l = ttstep(pooltype,i) + 1
        IF(riverstate%qqueue(l,pooltype,i)>0)THEN
          waterfrac = riverstate%qqueue(l,pooltype,i)/totvol    !Note whole volume so that pool get correct concentration change
          CALL add_source_to_water(riverstate%qqueue(l,pooltype,i),numsubstances,riverstate%cqueue(:,l,pooltype,i),waterfrac*sourcedd)
        ENDIF
      ENDIF
    ENDIF

  END SUBROUTINE add_dry_deposition_to_river

  !>\brief Calculate nutrient and sediment processes in river
  !> This include denitrification, mineralisation, primary production,
  !!sedimentation, resuspension, exchange with sediment
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Common things in lakes and river,
  !> Denitrification, Primary production and mineralization, and Sedimentation/Resuspension)
  !---------------------------------------------------------------------------
  SUBROUTINE np_processes_in_river(i,itype,area,depth,transq,Qbank,denpar, &
                                   denparl,hsatINpar,prodNpar,prodPpar,hsatTPpar, &
                                   sedexppar,limpppar,riverstate)

    USE MODVAR, ONLY : i_in,i_sp,conductP,conductS

    !Argument declarations
    INTEGER, INTENT(IN) :: i         !<index of current subbasin
    INTEGER, INTENT(IN) :: itype     !<river type (local or main)
    REAL, INTENT(IN)    :: area      !<river area (m2)
    REAL, INTENT(IN)    :: depth     !<river depth (m)
    REAL, INTENT(IN)    :: transq    !<flow out of translation box chain (m3/s)
    REAL, INTENT(IN)    :: qbank     !<bank full river flow
    REAL, INTENT(IN)    :: denpar    !<model parameter denitrification rate (kg/m2/day)
    REAL, INTENT(IN)    :: denparl   !<model parameter denitrification rate, local river (kg/m2/day)
    REAL, INTENT(IN)    :: hsatINpar !<model parameter half saturation IN (mg/L)
    REAL, INTENT(IN)    :: prodNpar  !<model parameter production ON
    REAL, INTENT(IN)    :: prodPpar  !<model parameter production PP
    REAL, INTENT(IN)    :: hsatTPpar !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)    :: sedexppar !<sedimentation/resuspension parameter (mg/L)
    REAL, INTENT(IN)    :: limpppar  !<limitation of sedimentation parameter (mg/L)
    TYPE(riverstatetype),INTENT(INOUT) :: riverstate  !<River states

    !Local parameters
    INTEGER, PARAMETER :: systemtype = 2    !river system

    !Calculate the nutrient processes
    IF(area>0)THEN
      IF(i_sp>0) CALL calculate_river_tpmean(i,itype,riverstate)
      IF(i_in>0) THEN
        IF(itype==1)THEN
          CALL denitrification_water(i,itype,systemtype,area,denparl,hsatINpar,RIVERSTATE=riverstate)  !denitrification in local rivers
        ELSE
          CALL denitrification_water(i,itype,systemtype,area,denpar,hsatINpar,RIVERSTATE=riverstate)  !denitrification in rivers
        ENDIF
      ENDIF
      CALL production_mineralisation(i,itype,systemtype,area,prodNpar,prodPpar,hsatTPpar,limpppar,RIVERSTATE=riverstate,DEPTH=depth) !mineraliation and primary production in rivers
      IF(conductP.OR.conductS) CALL sedimentation_resuspension(i,itype,area,sedexppar,transq,Qbank,depth,riverstate) !sedimentation and resuspension of part P in rivers
       !    No sediment SRP exchange. The concentration smoothed with deadvolume
    ENDIF

  END SUBROUTINE np_processes_in_river

  !>\brief Calculate nutrient and sediment processes in lake:
  !!denitrification, mineralisation, primary production,
  !!sedimentation, internal load
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Common things in lakes and river,
  !> Denitrification, Primary production and mineralization, Sedimentation/Resuspension and Internal load)
  !------------------------------------------------------------------
  SUBROUTINE np_processes_in_lake(i,itype,area,denpar,hsatINpar,prodNpar, &
                                  prodPpar,hsatTPpar,sedonpar,sedpppar,sedsspar,sedaepar,  &
                                  limonpar,limpppar,limsspar,lakestate)

    USE MODVAR, ONLY : conductN,conductP,conductS, &
                       i_on,i_pp,i_ss,i_ae

    !Argument declarations
    INTEGER, INTENT(IN) :: i          !<index of subbasin
    INTEGER, INTENT(IN) :: itype      !<lake type (ilake or olake)
    REAL, INTENT(IN)    :: area       !<lake area (m2)
    REAL, INTENT(IN)    :: denpar     !<model parameter denitrification rate (kg/m2/day)
    REAL, INTENT(IN)    :: hsatINpar  !<model parameter half saturation IN (mg/L)
    REAL, INTENT(IN)    :: prodNpar   !<model parameter production ON
    REAL, INTENT(IN)    :: prodPpar   !<model parameter production PP
    REAL, INTENT(IN)    :: sedonpar   !<ON sedimentation rate  (lakes)
    REAL, INTENT(IN)    :: sedpppar   !<PP sedimentation rate  (lakes)
    REAL, INTENT(IN)    :: sedsspar   !<SS sedimentation rate  (lakes)
    REAL, INTENT(IN)    :: sedaepar   !<AE sedimentation rate  (lakes)
    REAL, INTENT(IN)    :: hsatTPpar  !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)    :: limonpar   !<limitation of sedimentation ON parameter (mg/L)
    REAL, INTENT(IN)    :: limpppar   !<limitation of sedimentation PP parameter (mg/L)
    REAL, INTENT(IN)    :: limsspar   !<limitation of sedimentation SS parameter (mg/L)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate  !<Lake state

    !Local parameters
    INTEGER, PARAMETER :: systemtype = 1    !lake

    !Calculate the nutrient processes
    IF(conductP) CALL calculate_lake_tpmean(i,itype,lakestate)
    IF(conductN)THEN
      CALL denitrification_water(i,itype,systemtype,area,denpar,hsatINpar,LAKESTATE=lakestate) !denitrification in lakes
    ENDIF
    CALL production_mineralisation(i,itype,systemtype,area,prodNpar,prodPpar,hsatTPpar,limpppar,LAKESTATE=lakestate)  !primary production and mineralisation in lakes
    IF(conductN) CALL lake_sedimentation(i,i_on,itype,sedonpar,limonpar,lakestate)
    IF(conductP) CALL lake_sedimentation(i,i_pp,itype,sedpppar,limpppar,lakestate)
    IF(conductS)THEN
      CALL lake_sedimentation(i,i_ss,itype,sedsspar,limsspar,lakestate)
      CALL lake_sedimentation(i,i_ae,itype,sedaepar,0.,lakestate)
    ENDIF
    IF(conductP) CALL internal_lake_load(i,itype,systemtype,area,lakestate)  !internal load of phosphorus

  END SUBROUTINE np_processes_in_lake

  !>\brief Calculates the denitrification in river and lakes
  !!Lake processes in slow turnover lake part.
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Denitrification)
  !-----------------------------------------------------------------------
  SUBROUTINE denitrification_water(i,watertype,systemtype,area,denpar,halfsatINwater,riverstate,lakestate)

    USE MODVAR, ONLY : i_in
    USE HYPEVARIABLES, ONLY : maxdenitriwater

    !Argument declarations
    INTEGER, INTENT(IN) :: i          !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype  !<Lake or river type (1=local, 2=main/outlet)
    INTEGER, INTENT(IN) :: systemtype !<aquatic system type (1=lake, 2=river)
    REAL, INTENT(IN)    :: area       !<lake surface area/river bottom area (m2)
    REAL, INTENT(IN)    :: denpar     !<model parameter denitrification rate (kg/m2/day)
    REAL, INTENT(IN)    :: halfsatINwater  !<model parameter half saturation IN (mg/L)
    TYPE(riverstatetype),INTENT(INOUT),OPTIONAL :: riverstate  !<River states
    TYPE(lakestatetype),INTENT(INOUT),OPTIONAL  :: lakestate   !<Lake states

    !Local variables
    REAL, DIMENSION(1) :: denitri_water, inorganicNpool
    REAL tmpfcn, concfcn, watertemp, waterconc,vol

    !Local parameters
    INTEGER, PARAMETER :: pooldim = 1

    !Initial pools and values
    IF(systemtype==1) THEN   !lakes
      vol = lakestate%slowwater(watertype,i) * area * 1.0E-6    !0.001 m3
      waterconc = lakestate%concslow(i_in,watertype,i)          !mg/L
      inorganicNpool = vol * waterconc                          !kg
      watertemp = lakestate%temp(watertype,i)
    ELSE                     !rivers
      vol = riverstate%water(watertype,i) * 1.0E-3
      waterconc = riverstate%conc(i_in,watertype,i)
      inorganicNpool = vol * waterconc     !kg
      watertemp = riverstate%temp(watertype,i)
    ENDIF

    !Temperature and concentration dependence factor
    tmpfcn  = tempfactor(watertemp)
    concfcn = halfsatconcfactor(waterconc,halfsatINwater)

    !Denitrification
    denitri_water = denpar * area * concfcn * tmpfcn   !kg
    denitri_water = MIN(maxdenitriwater*inorganicNpool, denitri_water)    !max 50% kan be denitrified
    CALL retention_pool(pooldim, inorganicNPool, denitri_water)
    IF(systemtype==1) THEN   !lakes
      CALL new_concentration(inorganicNpool(1),vol,lakestate%concslow(i_in,watertype,i))
    ELSE                     !rivers
      IF(riverstate%water(watertype,i) > 0.) THEN
        CALL new_concentration(inorganicNpool(1),vol,riverstate%conc(i_in,watertype,i))
      ENDIF
    ENDIF

  END SUBROUTINE denitrification_water

  !>\brief Calculates transformation between IN/ON and SRP/PP in water.
  !!Also simulated algae production (=ON-production) for sediment simulation.
  !!Simulating the combined processes of primary production and
  !!mineralisation. Lake process in slow turnover lake part.
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes
  !>in rivers and lakes (Primary production and mineralization)
  !-----------------------------------------------------------------------
  SUBROUTINE production_mineralisation(i,watertype,systemtype,area, &
                                       prodNpar,prodPpar,halfsatTPwater,  &
                                       limpppar,riverstate,lakestate,depth)

    USE MODVAR, ONLY : conductN,conductP,conductS, &
                       i_in,i_on,i_sp,i_pp,i_ae
    USE HYPEVARIABLES, ONLY : maxprodwater,   &
                              maxdegradwater, &
                              NPratio

    !Argument declarations
    INTEGER, INTENT(IN) :: i          !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype  !<Lake or river type (1=local, 2=main/outlet)
    INTEGER, INTENT(IN) :: systemtype !<aquatic system type (1=lake, 2=river)
    REAL, INTENT(IN)    :: area       !<lake surface area/ river bottom area (m2)
    REAL, INTENT(IN)    :: prodNpar   !<model parameter production rate ON in water
    REAL, INTENT(IN)    :: prodPpar   !<model parameter production rate PP in water
    REAL, INTENT(IN)    :: halfsatTPwater !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)    :: limpppar   !<limitation of sedimentation parameter (mg/L)
    TYPE(riverstatetype),INTENT(INOUT),OPTIONAL :: riverstate !<River states
    TYPE(lakestatetype),INTENT(INOUT),OPTIONAL  :: lakestate  !<Lake states
    REAL, INTENT(IN), OPTIONAL :: depth      !<river depth (m)

    !Local variables
    REAL, DIMENSION(1) :: ONpool, INpool, SRPpool,PPpool,AEpool
    REAL, DIMENSION(1) :: minprodN, minprodP,minprodAE
    REAL watertemp,waterTPmean,temp10,temp20
    REAL tmpfcn, tmpfcn1, tmpfcn2, TPfcn
    REAL vol
    REAL waterdepth     !m

    !Local parameters
    INTEGER, PARAMETER :: pooldim = 1

    IF(.NOT.(conductN.OR.conductP.OR.conductS)) RETURN

    !Pools of nutrients in the water, water temperature and fraction of depth of water volume that is active
    IF (systemtype==1) THEN   !lakes
      vol = lakestate%slowwater(watertype,i) * area / 1.0E6
      IF(conductN) THEN
        INpool = vol * lakestate%concslow(i_in,watertype,i)   !kg
        ONpool = vol * lakestate%concslow(i_on,watertype,i)   !kg
      ENDIF
      IF(conductP) THEN
        SRPpool = vol * lakestate%concslow(i_sp,watertype,i) !kg
        PPpool  = vol * lakestate%concslow(i_pp,watertype,i)  !kg
      ENDIF
      IF(conductS) AEpool = vol * lakestate%concslow(i_ae,watertype,i) !kg N
    ELSE                     !rivers
      vol = riverstate%water(watertype,i) / 1.0E3
      IF(conductN) THEN
        INpool = vol * riverstate%conc(i_in,watertype,i) !kg
        ONpool = vol * riverstate%conc(i_on,watertype,i) !kg
      ENDIF
      IF(conductP) THEN
        SRPpool = vol * riverstate%conc(i_sp,watertype,i)    !kg
        PPpool  = vol * riverstate%conc(i_pp,watertype,i)    !kg
      ENDIF
      IF(conductS) AEpool = vol * riverstate%conc(i_ae,watertype,i)    !kg N
    ENDIF

    !Set help variables
    IF (systemtype==1) THEN   !lakes
      watertemp = lakestate%temp(watertype,i)
      waterdepth = lakestate%slowwater(watertype,i) * 1.E-3
      waterTPmean = lakestate%TPmean(watertype,i)
      temp10 = lakestate%temp10(watertype,i)
      temp20 = lakestate%temp20(watertype,i)
    ELSE                     !rivers
      watertemp = riverstate%temp(watertype,i)
      waterdepth = depth
      waterTPmean = riverstate%TPmean(watertype,i)
      temp10 = riverstate%temp10(watertype,i)
      temp20 = riverstate%temp20(watertype,i)
    ENDIF

    IF(watertemp > 0.)THEN
      !Total phosphorus concentration dependent factor
      TPfcn = halfsatconcfactor(MAX(waterTPmean-limpppar,0.),halfsatTPwater)

      !Temperature dependent factor
      tmpfcn1 = watertemp / 20.
      tmpfcn2 = (temp10 - temp20) / 5.
      tmpfcn = tmpfcn1*tmpfcn2

      !Production/mineralisation
      IF(conductN.OR.conductS)THEN
        minprodN = prodNpar * TPfcn * tmpfcn * waterdepth * area  !kg
        IF(conductN)THEN
          IF(minprodN(1) > 0.) THEN  !production
            minprodN = MIN(maxprodwater * INpool, minprodN)
          ELSE                       !mineralisation
            minprodN = MAX(-maxdegradwater * ONpool, minprodN)
          ENDIF
          CALL retention_pool(pooldim,INpool,minprodN)   !minprodN may be negative
          CALL production_pool(pooldim,ONpool,minprodN)
        ENDIF
        IF(conductS)THEN
          minprodAE = MAX(-AEpool, minprodN)
          CALL production_pool(pooldim,AEpool,minprodAE) !AEpool is a part of ON pool
        ENDIF
      ENDIF
      IF(conductP)THEN
        minprodP = prodPpar * NPratio * TPfcn * tmpfcn * waterdepth * area  !kg
        IF(minprodP(1) > 0.) THEN  !production
          minprodP = MIN(maxprodwater * SRPpool,minprodP)
        ELSE                       !mineralisation
          minprodP = MAX(-maxdegradwater * PPpool,minprodP)
        ENDIF
        CALL retention_pool(pooldim,SRPpool,minprodP)    !minprodP may be negative
        CALL production_pool(pooldim,PPpool,minprodP)
      ENDIF

      !New concentration due to changes in pools
      IF(systemtype==1) THEN            !lakes
        IF(conductN) CALL new_concentration(INpool(1),vol,lakestate%concslow(i_in,watertype,i))
        IF(conductN) CALL new_concentration(ONpool(1),vol,lakestate%concslow(i_on,watertype,i))
        IF(conductP) CALL new_concentration(SRPpool(1),vol,lakestate%concslow(i_sp,watertype,i))
        IF(conductP) CALL new_concentration(PPpool(1),vol,lakestate%concslow(i_pp,watertype,i))
        IF(conductS) CALL new_concentration(AEpool(1),vol,lakestate%concslow(i_ae,watertype,i))
      ELSE                              !rivers
        IF(riverstate%water(watertype,i) > 0.) THEN
          IF(conductN) CALL new_concentration(INpool(1),vol,riverstate%conc(i_in,watertype,i))
          IF(conductN) CALL new_concentration(ONpool(1),vol,riverstate%conc(i_on,watertype,i))
          IF(conductP) CALL new_concentration(SRPpool(1),vol,riverstate%conc(i_sp,watertype,i))
          IF(conductP) CALL new_concentration(PPpool(1),vol,riverstate%conc(i_pp,watertype,i))
          IF(conductS) CALL new_concentration(AEpool(1),vol,riverstate%conc(i_ae,watertype,i))
        ENDIF
      ENDIF
    ENDIF

  END SUBROUTINE production_mineralisation

  !>\brief Calculate sedimentation of substance in lake
  !!Only substance in in slow turnover lake part sediments.
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in
  !>rivers and lakes (Sedimentation/Resuspension) and Organic carbon (River and lakes - Sedimentation)
  !--------------------------------------------------------------------------
  SUBROUTINE lake_sedimentation(i,substance,watertype,sedrate,limsedpar,lakestate)

    !Argument declarations
    INTEGER, INTENT(IN) :: i          !<current index of subbasin
    INTEGER, INTENT(IN) :: substance  !<current index of substance (PP,ON,SS,AE)
    INTEGER, INTENT(IN) :: watertype  !<Lake type (1=local, 2=outlet)
    REAL, INTENT(IN)    :: sedrate    !<sedimentation rate  (lakes) (m/ts)
    REAL, INTENT(IN)    :: limsedpar  !<concentration limit for sedimentation (mg/L)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate  !<Lake state

    !Local variables
    REAL, DIMENSION(1) :: pool     !substance pool in water per area (g/m2) ((kg))
    REAL, DIMENSION(1) :: delta    !sedimentation (g/m2/ts)
    REAL depth                     !lake average depth (m)

    !>/b Algoritm
    !>Calculate water volume, amount of substance in water (pool) and the sedimentation
    depth = lakestate%slowwater(watertype,i)  * 1.0E-3        !m
    pool = depth * lakestate%concslow(substance,watertype,i) !mg/l*m
    delta = sedrate * MAX(lakestate%concslow(substance,watertype,i)-limsedpar,0.) !m/ts*mg/l

    !>Remove sedimentation from the pool
    CALL retention_pool(1,pool,delta)

    !>Calculate the new concentration in the water due to the change in the pool
    CALL new_concentration(pool(1),depth,lakestate%concslow(substance,watertype,i))

  END SUBROUTINE lake_sedimentation

  !>\brief Calculate sedimentation and resuspension of PP and SS in rivers.
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Sedimentation/Resuspension)
  !-----------------------------------------------------------------
  SUBROUTINE sedimentation_resuspension(i,watertype,area,sedexppar,riverq,qbank,depth,riverstate)

    USE MODVAR, ONLY : i_pp,i_ss

    !Argument declaration
    INTEGER, INTENT(IN) :: i          !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype  !<river type (1=local, 2=main)
    REAL, INTENT(IN)    :: area       !<river surface area (m2)
    REAL, INTENT(IN)    :: sedexppar  !<sedimentation/resuspension parameter
    REAL, INTENT(IN)    :: riverq     !<river discharge (m3/s)
    REAL, INTENT(IN)    :: qbank      !<bank full flow (m3/s)
    REAL, INTENT(IN)    :: depth      !<river depth (m)
    TYPE(riverstatetype),INTENT(INOUT) :: riverstate  !<River states

    !Local variables
    REAL, DIMENSION(1) :: PPpool,SSpool    !pools in water (kg)
    REAL, DIMENSION(1) :: transfer         !changes (kg/d)
    REAL, DIMENSION(1) :: tempsed,tempsed2 !temparary variables for river sediment pools (kg)
    REAL sedresp, help, qbankcorr

    !>/b Algoritm
    !>Initial check if calculation is to be made
    IF(sedexppar == 0) RETURN
    IF(area>0. .AND. qbank>0.) THEN

      !>Assess pool of PP and SS in water and sediment
      IF(i_pp>0)THEN
        tempsed(1) = riverstate%Psed(watertype,i)
        PPpool = (riverstate%water(watertype,i) * riverstate%conc(i_pp,watertype,i))* 1.0E-3 !kg
      ENDIF
      IF(i_ss>0)THEN
        tempsed2(1) = riverstate%Ssed(watertype,i)
        SSpool = (riverstate%water(watertype,i) * riverstate%conc(i_ss,watertype,i))* 1.0E-3 !kg
      ENDIF

      !>Calculate sedimentation and resuspension factor
      qbankcorr = 0.7*qbank
      help = 0.
      IF(qbankcorr-riverq>0.) help = help + ((qbankcorr-riverq)/qbankcorr)**sedexppar   !sedimentation at low flow
      IF(riverq>0) help = help - (riverq/qbankcorr)**sedexppar  !resuspension at all flows
      sedresp = MAX(-1., MIN(1.,help))

      !>Transfer PP and SS between sediment and water pools
      IF(sedresp>0.)THEN  !sedimentation
        IF(i_pp>0)THEN
          transfer = sedresp * (riverstate%conc(i_pp,watertype,i) * MIN(riverstate%water(watertype,i),area * depth)) / 1.0E3
          CALL retention_pool(1,PPpool,transfer)           !transfer may change
          CALL production_pool(1,tempsed,transfer)
        ENDIF
        IF(i_ss>0)THEN
          transfer = sedresp * (riverstate%conc(i_ss,watertype,i) * MIN(riverstate%water(watertype,i),area * depth)) / 1.0E3
          CALL retention_pool(1,SSpool,transfer)           !transfer may change
          CALL production_pool(1,tempsed2,transfer)
        ENDIF
      ELSE                !resuspension
        IF(i_pp>0)THEN
          transfer = - sedresp * tempsed
          CALL retention_pool(1,tempsed,transfer)          !transfer may change
          CALL production_pool(1,PPpool,transfer)
        ENDIF
        IF(i_ss>0)THEN
          transfer = - sedresp * tempsed2
          CALL retention_pool(1,tempsed2,transfer)         !transfer may change
          CALL production_pool(1,SSpool,transfer)
        ENDIF
      ENDIF

      !>Update state variables
      IF(i_pp>0)THEN
        riverstate%Psed(watertype,i) = tempsed(1)
        CALL new_concentration(PPpool(1),riverstate%water(watertype,i)*1.0E-3,riverstate%conc(i_pp,watertype,i))
      ENDIF
      IF(i_ss>0)THEN
        riverstate%Ssed(watertype,i) = tempsed2(1)
        CALL new_concentration(SSpool(1),riverstate%water(watertype,i)*1.0E-3,riverstate%conc(i_ss,watertype,i))
      ENDIF
    ENDIF

  END SUBROUTINE sedimentation_resuspension

  !>\brief Calculates straight 365-day running average mean of TP
  !>concentration in lake
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Primary production
  !> and mineralization) and Organic carbon (River and Lakes - Primary production and mineralization)
  !-----------------------------------------------------------------------
  SUBROUTINE calculate_lake_tpmean(i,watertype,lakestate)

    USE MODVAR, ONLY : i_sp,i_pp

    !Argument declarations
    INTEGER, INTENT(IN) :: i         !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype !<Lake type (1=local, 2=outlet)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate  !<Lake state

    !Local variables
    REAL meanconc

    IF(lakestate%water(watertype,i)+lakestate%slowwater(watertype,i)>0.)THEN
      meanconc = (lakestate%water(watertype,i)*(lakestate%conc(i_sp,watertype,i)+lakestate%conc(i_pp,watertype,i)) +           &
           lakestate%slowwater(watertype,i)*(lakestate%concslow(i_sp,watertype,i) + lakestate%concslow(i_pp,watertype,i)))/  &
           (lakestate%water(watertype,i)+lakestate%slowwater(watertype,i))
      lakestate%TPmean(watertype,i) = lakestate%TPmean(watertype,i) + (meanconc - lakestate%TPmean(watertype,i))/365.
    ENDIF

  END SUBROUTINE calculate_lake_tpmean

  !>\brief Calculates straight 365-day running average mean of TP
  !>concentration in river
  !>
  !>\b Reference ModelDescription Chapter  Nitrogen and phosphorus processes in rivers and lakes (Primary production
  !> and mineralization) and Organic carbon (River and Lakes - Primary production and mineralization)
  !-------------------------------------------------------------------
  SUBROUTINE calculate_river_tpmean(i,watertype,riverstate)

    USE MODVAR, ONLY : i_sp,i_pp

    !Argument declarations
    INTEGER, INTENT(IN) :: i         !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype !<River type (1=local, 2=main)
    TYPE(riverstatetype),INTENT(INOUT) :: riverstate  !<River states

    riverstate%TPmean(watertype,i) = riverstate%TPmean(watertype,i) + (riverstate%conc(i_sp,watertype,i) &
         + riverstate%conc(i_pp,watertype,i) - riverstate%TPmean(watertype,i))/365.

  END SUBROUTINE calculate_river_tpmean

  !>\brief Calculates and add internal load of phosphorus for lakes
  !!Lake processes in slow turnover lake part
  !>
  !>\b Reference ModelDescription Chapter Nitrogen and phosphorus processes in rivers and lakes (Internal load)
  !-------------------------------------------------------------
  SUBROUTINE internal_lake_load(i,watertype,systemtype,area,lakestate)

    USE MODVAR, ONLY : lakeindex,   &
                       lakedatapar, &
                       lakedataparindex,  &
                       i_sp,i_pp,   &
                       numsubstances
    USE HYPEVARIABLES, ONLY : m_ldprodpp,  &
                              m_ldprodsp

    !Argument declarations
    INTEGER, INTENT(IN) :: i           !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype   !<Lake or river type (1=local, 2=main/outlet)
    INTEGER, INTENT(IN) :: systemtype  !<aquatic system type (1=lake, 2=river)
    REAL, INTENT(IN)    :: area        !<lake surface area/ river bottom area (m2)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate  !<Lake state

    !Local variables
    REAL prodPP, prodSP
    REAL tmpfcn, TPfcn
    REAL vol
    REAL pppar,sppar
    REAL :: sourceP(numsubstances)
    !Local parameters
    INTEGER, PARAMETER :: pooldim = 1

    !>\b Algorithm \n
    !>Check if internal phosphorus load is to be calculated
    IF(systemtype==2) RETURN   !river
    IF(watertype==1) RETURN    !local
    IF(.NOT.ALLOCATED(lakeindex)) RETURN  !no special lakes
    pppar = lakedatapar(lakedataparindex(i,watertype),m_ldprodpp)
    sppar = lakedatapar(lakedataparindex(i,watertype),m_ldprodsp)
    IF(pppar==0 .AND. sppar==0) RETURN
    sourceP = 0.

    !>Calculate pool of P, and concentration and temperature dependent factors
    TPfcn = 0.1
    tmpfcn = 0.86**(ABS(lakestate%temp(watertype,i)-15.))   !laketemp=T20 for olake

    !> Calculate internal load of phosphorus
    prodPP = pppar * TPfcn * tmpfcn * area / 1000.  !kg/d
    prodSP = sppar * TPfcn * tmpfcn * area / 1000.  !kg/d
    sourceP(i_pp) = prodPP
    sourceP(i_sp) = prodSP

    !>Add internal load of phosphorus to lake water
    vol = lakestate%slowwater(watertype,i) * area / 1.0E6
    CALL add_source_to_water(vol,numsubstances,lakestate%concslow(:,watertype,i),sourceP)

  END SUBROUTINE internal_lake_load

  !>Calculate organic carbon processes in river; mineralisation,
  !>primary production
  !>
  !>\b Reference ModelDescription Chapter Organic carbon (River and Lakes - Primary production and mineralization)
  !----------------------------------------------------------------
  SUBROUTINE oc_processes_in_river(i,watertype,area,depth,prodpar,hsatTPpar,limpppar,riverstate)

    USE MODVAR, ONLY : i_oc

    !Argument declarations
    INTEGER, INTENT(IN) :: i          !<index of current subbasin
    INTEGER, INTENT(IN) :: watertype  !<river type (local or main)
    REAL, INTENT(IN)    :: area       !<river area (m2)
    REAL, INTENT(IN)    :: depth      !<river depth (m)
    REAL, INTENT(IN)    :: prodpar    !<model parameter production OC
    REAL, INTENT(IN)    :: hsatTPpar  !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)    :: limpppar   !<limitation of sedimentation parameter (mg/L)
    TYPE(riverstatetype),INTENT(INOUT) :: riverstate  !<River states

    !Local parameters
    INTEGER, PARAMETER :: systemtype = 2    !river system

    IF(i_oc==0)RETURN

    !Calculate the organic carbon processes
    IF(area>0)THEN
      CALL oc_production_mineralisation(systemtype,area,prodpar,hsatTPpar,limpppar,       &
              riverstate%water(watertype,i),riverstate%conc(i_oc,watertype,i),  &
              riverstate%temp(watertype,i),riverstate%TPmean(watertype,i),      &
              riverstate%temp10(watertype,i),riverstate%temp20(watertype,i),depth)
    ENDIF

  END SUBROUTINE oc_processes_in_river

  !>\brief Calculate organic carbon processes in lake
  !!Mineralisation, primary production, sedimentation
  !>
  !>\b Reference ModelDescription Chapter Organic carbon (River and Lakes)
  !------------------------------------------------------------------
  SUBROUTINE oc_processes_in_lake(i,watertype,area,prodpar,hsatTPpar,limpppar,sedocpar,lakestate)

    USE MODVAR, ONLY : i_oc

    !Argument declarations
    INTEGER, INTENT(IN) :: i                        !<current index of subbasin
    INTEGER, INTENT(IN) :: watertype                !<lake type (ilake or olake)
    REAL, INTENT(IN)    :: area                     !<lake area (m2)
    REAL, INTENT(IN)    :: prodpar                  !<model parameter production OC
    REAL, INTENT(IN)    :: hsatTPpar                !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)    :: limpppar                 !<limitation of sedimentation parameter (mg/L)
    REAL, INTENT(IN)    :: sedocpar                 !<OC sedimentation rate  (lakes)
    TYPE(lakestatetype),INTENT(INOUT) :: lakestate  !<Lake state

    !Local parameters
    INTEGER, PARAMETER :: systemtype = 1    !lake

    IF(i_oc==0) RETURN

    !Calculate the nutrient processes
    CALL oc_production_mineralisation(systemtype,area,prodpar,hsatTPpar,limpppar,   &
            lakestate%slowwater(watertype,i),lakestate%concslow(i_oc,watertype,i),  &
            lakestate%temp(watertype,i),lakestate%TPmean(watertype,i),              &
            lakestate%temp10(watertype,i),lakestate%temp20(watertype,i))
    CALL lake_sedimentation(i,i_oc,watertype,sedocpar,0.,lakestate)

  END SUBROUTINE oc_processes_in_lake

  !>\brief Calculates transformation between OC/DIC in water
  !!Simulating the combined processes of primary production and
  !!mineralisation.
  !>
  !>\b Reference ModelDescription Organic carbon (River and lakes - Primary production
  !> and mineralization)
  !----------------------------------------------------------------
  SUBROUTINE oc_production_mineralisation(systemtype,area,prodpar,halfsatTPwater, &
                             limpppar,water,conc,watertemp,waterTPmean,temp10,temp20,depth)

    USE HYPEVARIABLES, ONLY : maxdegradwater, &
                              NCratio

    !Argument declarations
    INTEGER, INTENT(IN)        :: systemtype  !<aquatic system type (1=lake, 2=river)
    REAL, INTENT(IN)           :: area        !<lake surface area/ river bottom area (m2)
    REAL, INTENT(IN)           :: prodpar     !<model parameter production rate OC in water
    REAL, INTENT(IN)           :: halfsatTPwater !<model parameter half saturation TP (mg/L)
    REAL, INTENT(IN)           :: limpppar    !<limitation of sedimentation parameter (mg/L)
    REAL, INTENT(IN)           :: water       !<river or lake water (mm or m3)
    REAL, INTENT(INOUT)        :: conc        !<OC concentration of river or lake
    REAL, INTENT(IN)           :: watertemp   !<water temperature
    REAL, INTENT(IN)           :: waterTPmean !<water TP mean
    REAL, INTENT(IN)           :: temp10      !<10-day water temperature
    REAL, INTENT(IN)           :: temp20      !<20-day water temperature
    REAL, INTENT(IN), OPTIONAL :: depth       !<river depth (m)

    !Local variables
    REAL, DIMENSION(1) :: OCpool, minprodN, minprodC,minC,prodC
    REAL tmpfcn, tmpfcn1, tmpfcn2, TPfcn
    REAL vol
    REAL waterdepth !(m)

    !Local parameter
    INTEGER, PARAMETER :: pooldim = 1

    !>\b Algorithm \n
    !>Calculate pools of organic carbon in the water, water temperature
    !>and fraction of depth of water volume that is active
    IF(systemtype==1) THEN   !lakes
      OCpool = (water * area * conc) /1.0E6  !kg
      waterdepth = water/1000.
    ELSE                     !rivers
      OCpool = (water * conc)/ 1.0E3 !kg
      waterdepth=depth
    ENDIF

    !>Calculate dependency factors (Tot-P and temperature)
    TPfcn = halfsatconcfactor(waterTPmean-limpppar,halfsatTPwater)
    IF(watertemp >= 0.) THEN
      tmpfcn1 = watertemp / 20.
    ELSE
      tmpfcn1 = 0.
    ENDIF
    tmpfcn2 = (temp10 - temp20) / 5.
    tmpfcn = tmpfcn1*tmpfcn2

    !>Calculate production/mineralisation of organic carbon
    minprodN = 0.
    IF(watertemp > 0. ) THEN
      minprodN = prodpar * TPfcn * tmpfcn * waterdepth * area  !kg
      IF(minprodN(1) > 0.) THEN  !production
        minprodC = minprodN * NCratio
      ELSE                       !mineralisation
        minprodC = MAX(-maxdegradwater * OCpool, minprodN * NCratio)
      ENDIF
    ENDIF
    minC = -minprodC
    prodC = minprodC
    IF(minprodC(1)>0.) CALL production_pool(pooldim,OCpool,prodC)
    IF(minprodC(1)<0.) CALL retention_pool(pooldim,OCpool,minC)

    !>Set new concentration due to changes in pools
    IF(systemtype==1) THEN            !lakes
      vol = water * area / 1.0E6
      CALL new_concentration(OCpool(1),vol,conc)
    ELSE                                 !rivers
      IF(water > 0.) THEN
        vol = water / 1.0E3
        CALL new_concentration(OCpool(1),vol,conc)
      ENDIF
    ENDIF

  END SUBROUTINE oc_production_mineralisation

  !>Add load from local diffuse sources to local river inflow
  !>
  !> \b Reference ModelDescription Chapter Nitrogen and phosphorus in land
  !!routines (Nutrient sources - Rural household diffuse source)
  !-----------------------------------------------------------------
  SUBROUTINE add_diffuse_source_to_local_river(i,qin,cin,source,addedflow)

    USE MODVAR, ONLY : i_in,i_on,i_sp,i_pp, &
                       load,                &
                       genpar,              &
                       numsubstances,       &
                       seconds_per_timestep
    USE HYPEVARIABLES, ONLY : m_locsoil

    !Argument declarations
    INTEGER, INTENT(IN) :: i                      !<index of subbasin
    REAL, INTENT(INOUT) :: qin                    !<flow in local river (m3/s)
    REAL, INTENT(INOUT) :: cin(numsubstances)     !<concentration of flow into local river (mg/L)
    REAL, INTENT(OUT)   :: source(numsubstances)  !<local source added to local river (kg/timestep)
    REAL, INTENT(OUT)   :: addedflow              !<added flow (m3/timestep)

    !Local variables
    REAL qhelp
    REAL qadd
    REAL cadd(numsubstances)

    !Initiation
    source = 0.
    cadd = 0.

    !> \b Algorithm \n
    !>Calculate diffuse source from rural households to local river
    addedflow = (1. - genpar(m_locsoil)) * load(i)%volloc   !m3/ts
    qadd = addedflow / seconds_per_timestep   !m3/s !fel vid korta tidsteg!
    IF(qadd>0)THEN
      qhelp = qadd * seconds_per_timestep * 1.E-3   !1000m3/timestep
      cadd = load(i)%locconc
      IF(i_in>0)THEN
        source(i_in) = cadd(i_in) * qhelp    !Diffuse load, ruralB, kg/timestep
        source(i_on) = cadd(i_on) * qhelp    !Diffuse load, ruralB, kg/timestep
      ENDIF
      IF(i_sp>0)THEN
        source(i_sp) = cadd(i_sp) * qhelp    !Diffuse load, ruralB, kg/timestep
        source(i_pp) = cadd(i_pp) * qhelp    !Diffuse load, ruralB, kg/timestep
      ENDIF

      !>Add diffuse source to inflow to local river flow
      IF(qin>0)THEN
        cin = (qin * cin + qadd * cadd)/(qin + qadd)
        qin = qin + qadd
      ELSE
        qin = qadd
        cin = cadd
      ENDIF
    ENDIF

  END SUBROUTINE add_diffuse_source_to_local_river

  !>Add load from point sources to main river inflow
  !>
  !> \b Reference ModelDescription Chapter Water management (Point sources)
  !-----------------------------------------------------------------
  SUBROUTINE add_point_sources_to_main_river(i,qin,cin,source,addedflow)

    USE MODVAR, ONLY : i_in,i_on,i_sp,i_pp, &
                       i_t1,i_t2,           &
                       max_pstype,          &
                       load,                &
                       numsubstances,       &
                       seconds_per_timestep

    !Argument declarations
    INTEGER, INTENT(IN) :: i                        !<index of subbasin
    REAL, INTENT(INOUT) :: qin                      !<flow into main river (m3/s)
    REAL, INTENT(INOUT) :: cin(numsubstances)       !<concentration of flow into main river (mg/L)
    REAL, INTENT(OUT)   :: source(numsubstances,max_pstype)  !<point sources added to main river (kg/timestep)
    REAL, INTENT(OUT)   :: addedflow                !<added flow (m3/timestep)

    !Local variables
    INTEGER k
    REAL divvolps
    REAL qadd
    REAL cadd(numsubstances)

    !Initiation
    source = 0.
    cadd = 0.
    qadd = 0.

    !Calculate source to be added to river
    DO k = 1,max_pstype
      qadd = qadd + load(i)%psvol(k)   !m3/s
    ENDDO
    addedflow = qadd * seconds_per_timestep
    IF(qadd>0)THEN
      divvolps = 1000./qadd/seconds_per_timestep                    !kg/ts,m3/s->mg/L
      DO k = 1,max_pstype
        IF(i_in>0)THEN
          cadd(i_in) = cadd(i_in) + load(i)%psload(k,i_in)
          cadd(i_on) = cadd(i_on) + load(i)%psload(k,i_on)
          source(i_in,k) = load(i)%psload(k,i_in)        !Point source k,IN (kg/timestep)
          source(i_on,k) = load(i)%psload(k,i_on)        !Point source k,ON
        ENDIF
        IF(i_sp>0)THEN
          cadd(i_sp) = cadd(i_sp) + load(i)%psload(k,i_sp)
          cadd(i_pp) = cadd(i_pp) + load(i)%psload(k,i_pp)
          source(i_sp,k) = load(i)%psload(k,i_sp)        !Point source k, SP (kg/timestep)
          source(i_pp,k) = load(i)%psload(k,i_pp)        !Point source k, PP
        ENDIF
        IF(i_t1>0)THEN
          cadd(i_t1) = cadd(i_t1) + load(i)%psload(k,i_t1)
          source(i_t1,k) = load(i)%psload(k,i_t1)        !Point source k, T1
        ENDIF
        IF(i_t2>0)THEN
          cadd(i_t2) = cadd(i_t2) + load(i)%psload(k,i_t2)
          source(i_t2,k) = load(i)%psload(k,i_t2)        !Point source k, T2
        ENDIF
      ENDDO
      cadd(:) = cadd(:) * divvolps    !mg/L

      !Add source to river
      IF(qin>0)THEN
        cin = (qin * cin + qadd * cadd)/(qin + qadd)
        qin = qin + qadd
      ELSE
        qin = qadd
        cin = cadd
      ENDIF
    ENDIF

  END SUBROUTINE add_point_sources_to_main_river

  !>Calculate effect of river wetland constructed for nutrient removal
  !>
  !>\b Reference ModelDescription Chapter Water management (Constructed wetlands)
  !-----------------------------------------------------------------
  SUBROUTINE calculate_river_wetland(i,itype,n,temp5,temp30,qin,cin,cwetland)

    USE MODVAR, ONLY : wetland,     &
                       seconds_per_timestep

    !Argument declarations
    INTEGER, INTENT(IN) :: i           !<index of subbasin
    INTEGER, INTENT(IN) :: itype       !<index of river type (local or main)
    INTEGER, INTENT(IN) :: n           !<number of substances
    REAL, INTENT(IN)    :: temp5       !<temperature (5-day-mean) (degree Celsius)
    REAL, INTENT(IN)    :: temp30      !<temperature (30-day-mean) (degree Celsius)
    REAL, INTENT(IN)    :: qin         !<flow into/out of river wetland (m3/s)
    REAL, INTENT(INOUT) :: cin(n)      !<concentration of flow into/out of river wetland (mg/L)
    REAL, INTENT(INOUT) :: cwetland(n) !<concentration of river wetland (mg/L)

    !Local variables
    REAL wetlandvol     !m3 (constant)
    REAL wetlandinflow  !m3/timestep

    !Start of calculations
    IF(wetland(i,itype)%area==0) RETURN   !no wetland

    wetlandvol = wetland(i,itype)%area * wetland(i,itype)%depth   !m3
    wetlandinflow = qin * wetland(i,itype)%part * seconds_per_timestep     !m3/timestep
    CALL calculate_wetland_np(n,wetlandinflow,cin,wetland(i,itype)%area,wetlandvol,cwetland,temp5,temp30)
    IF(qin>0) cin = cin * (1. - wetland(i,itype)%part) + cwetland * wetland(i,itype)%part           !New concentration

  END SUBROUTINE calculate_river_wetland

  !>\brief Calculate nutrient processes in river wetland.
  !!Retention is limited to 99.9% of the pool.
  !>
  !>\b Reference ModelDescription Chapter Water management (Constructed wetlands)
  !------------------------------------------------------------------------
  SUBROUTINE calculate_wetland_np(n,qin,cin,area,vol,cvol,temp5,temp30)

    USE MODVAR, ONLY : i_in,i_sp,i_pp,  &
                       conductN, conductP

    !Argument declarations
    INTEGER, INTENT(IN) :: n       !<number of substances
    REAL, INTENT(IN)    :: qin     !<flow into wetland (m3/d)
    REAL, INTENT(IN)    :: cin(n)  !<concentration of river flow (mg/l) (before and after wetland processes
    REAL, INTENT(IN)    :: area    !<area of wetland (m2)
    REAL, INTENT(IN)    :: vol     !<volume of wetland (m3)
    REAL, INTENT(INOUT) :: cvol(n) !<concentration of wetland volume (mg/l) (before and after wetland processes
    REAL, INTENT(IN)    :: temp5   !<temperature (5-day-mean) (degree Celsius)
    REAL, INTENT(IN)    :: temp30  !<temperature (30-day-mean) (degree Celsius)

    !Local variables
    REAL wetlandnutrient(n), wetlandconc(n)
    REAL retention(n)
    REAL retention_tp, production_tp
    REAL wetland_tp,srpfrac

    !Local parameters
    REAL, PARAMETER :: teta = 1.2
    REAL, PARAMETER :: tkoeff = 20.   !temperature coefficient (degree Celsius)
    REAL, PARAMETER :: inpar = 2.3    !model parameter for inorganic nitrogen retention (mm/d/degree Celsius)
    REAL, PARAMETER :: sedpar = 0.09  !model parameter for phosphorus sedimentation (m/d)
    REAL, PARAMETER :: uptpar = 0.1   !model parameter for phosphorus uptake (m/d)

    !Calculate the nutrient processes
    wetlandnutrient = vol*cvol+qin*cin         !g
    wetlandconc = wetlandnutrient /(vol+qin)   !mg/l
    retention = 0.
    IF(conductN)THEN
      IF(temp5>0) retention(i_in) = inpar * wetlandconc(i_in) * area * temp5 * 1.E-3         !g/d denitrification
      IF(retention(i_in)<0) retention(i_in) = 0.
      IF(retention(i_in)>0.999*wetlandnutrient(i_in)) retention(i_in) = 0.999 * wetlandnutrient(i_in)
    ENDIF
    IF(conductP)THEN
      retention_tp = sedpar * (wetlandconc(i_pp) + wetlandconc(i_sp)) * area                 !g/d sedimentation
      IF(retention_tp<0) retention_tp = 0.
      production_tp = uptpar * (cin(i_pp) + cin(i_sp)) * (teta ** (temp30 - tkoeff)) * area  !g/d uptake
      IF(production_tp<0) production_tp = 0.
      wetland_tp = wetlandnutrient(i_pp) + wetlandnutrient(i_sp)    !g
      IF(retention_tp - production_tp < 0.999 * wetland_tp)THEN
        srpfrac = wetlandnutrient(i_sp) / wetland_tp
        retention(i_sp) = srpfrac * (retention_tp - production_tp)
        retention(i_pp) = (1.-srpfrac) * (retention_tp - production_tp)
      ELSE
        retention_tp = 0.999 * wetland_tp
        IF(wetland_tp>0)THEN
          srpfrac = wetlandnutrient(i_sp)/wetland_tp
        ELSE
          srpfrac = 0.
        ENDIF
        retention(i_sp) = srpfrac * retention_tp
        retention(i_pp) = (1.-srpfrac) * retention_tp
      ENDIF
    ENDIF
    cvol = (wetlandnutrient - retention)/(vol+qin)    !New concentration of wetland volume

  END SUBROUTINE calculate_wetland_np

END MODULE NPC_SURFACEWATER_PROCESSES

max_stars_repo_name: hydrology-tep/fanfar-postprocessing

max_issues_repo_name: hydrology-tep/fanfar-postprocessing

max_forks_repo_name: hydrology-tep/fanfar-postprocessing


text:
!*******************************************************************
!                                                                  *
      SUBROUTINE EL4INT(JJA,JJB,JA,JB,JC,JD,ICOLBREI,INTERACT)
!                                                                  *
!   --------------  SECTION METWO    SUBPROGRAM 09  -------------  *
!                                                                  *
!     THIS PACKAGE DETERMINES THE VALUES OF MATRIX ELEMENTS        *
!     OF TWO PARTICLE OPERATOR IN CASE :      N'1 = N1 +- 1        *
!                                             N'2 = N2 +- 1        *
!                                             N'3 = N3 -+ 2        *
!                                                                  *
!     SUBROUTINE CALLED: EL41                                      *
!                                                                  *
!   Written by  G. Gaigalas                   NIST, December 2015  *
!                                                                  *
!*******************************************************************
!
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE m_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE el41INT_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER, INTENT(IN) :: JJA,JJB,JA,JB,JC,JD,ICOLBREI
      INTEGER, INTENT(OUT) :: INTERACT
!-----------------------------------------------
      INTERACT=0
      IF(NPEEL.LE.2)RETURN
      IF(JA.EQ.JB) THEN
        CALL EL41INT(JJA,JJB,JC,JD,JA,1,JA,JB,JC,JD,ICOLBREI,INTERACT)
      ELSE IF(JC.EQ.JD) THEN
        CALL EL41INT(JJA,JJB,JA,JB,JC,2,JA,JB,JC,JD,ICOLBREI,INTERACT)
      ELSE
        WRITE(99,100)
        STOP
      END IF
      RETURN
  100 FORMAT(5X,'ERRO IN EL4INT ')
      END SUBROUTINE EL4INT

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
Background: I lived in Davis from 19901994, then moved to Sonoma County where I met Users/EthanKnoop Ethan Knoop. I moved back to Davis in 2000. I am proud to call myself a townie!

Fur Kids: My first born is a 7 year old cockatiel named Sir Peppermint Waddlesworth. Shortly after I started working in a vet hospital, I acquired two calicos named Teeny and Camilla.

Employment: I have had various jobs around town. I used to work at Tecate Grill and Animal Kingdom. I am currently working at one of the Veterinary Hospitals in town as a Registered Veterinary Technician.

Interests: Video games, board games, music, reading, movies, food and sleeping. (Not necessarily in that order!)

Favorite places in town: Stonegate Country Club, Thai Recipes, Peets Coffee and Tea, West Pond, Sweet Briar Books, Let Them Eat Cake.

Favorite departed businesses: Discoveries, Mirai, Crepe House Uni.

20090721 16:58:14 nbsp Do you have a younger sister, Zofia? Users/WilliamLewis

20090721 17:05:58 nbsp No, sorry! Users/EmmyMelton

20120114 16:27:07 nbsp Hi, did you ever end up getting a ticket for running the light on Russel? Users/JLee05

20120114 17:11:39 nbsp No, I didnt. I waited a few months and nothing showed up in the mail. I was paranoid that maybe it got lost in the mail so I went into the police station just to make sure. The lady told me I probably got lucky and told me to never run a red light again. Of course I dont plan on running red lights on purpose! Users/EmmyMelton

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! { dg-do run }
! PR 21912
! We didn't adapt the exit condition to negative steps in array constructors,
! leaving the resulting arrays uninitialized.
integer :: i(5), n, m, l, k

n = 5
i = (/ (m, m = n, 1, -1) /)
if (any (i /= (/ 5, 4, 3, 2, 1 /))) call abort

k = 1

i(5:1:-1) = (/ (m, m = n, k, -1) /)
if (any (i /= (/ 1, 2, 3, 4, 5 /))) call abort

l = -1

i = (/ (m, m = n, 1, l) /)
if (any (i /= (/ 5, 4, 3, 2, 1 /))) call abort

i(5:1:-1) = (/ (m, m = n, k, l) /)
if (any (i /= (/ 1, 2, 3, 4, 5 /))) call abort
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


SUBROUTINE CONTINUITY_MOMENTUM_X(PP1,PP2,HH1,HH2,QQ1,QQ2,T0,DT,DX,DY,I1,I2,IP,UPS)
USE TEST
USE FLOATPRECISION
USE PARAMETERS
USE DEPTH
IMPLICIT NONE
REAL(FP) , DIMENSION(:,:) , POINTER   :: PP1 , PP2 , HH1 , HH2 , QQ1 , QQ2
REAL(FP)                              :: T0  , DT  , DX  , DY
INTEGER                               :: I1  , I2  , IP
LOGICAL                               :: UPS , DWS

REAL(FP) :: FROUDE , WF
REAL(FP) :: ADVE , ADVW
REAL(FP) :: VN , VS , PN , PS , CADL , CADR , CDIFFL , CDIFFR
REAL(FP) :: RESIST , CHEZY , QAV , HAV , GRV , HAV0
REAL(FP) :: HE , HW , HN , HS
REAL(FP) :: DIFFX , DIFFY
INTEGER  :: ML , MU , NL , NU
INTEGER  :: I , J , K , JC , JJC , J0
REAL(FP) :: AMT(6,3001)
REAL(FP) :: SEABND
REAL(FP) :: DTDY , DTDX , DXDY , DTDXDY , DTDXDDY , DTDYDDX , DXDY4 , DT0125 , DT2 , DTDY2
REAL(FP) :: POW1 , GRA05 , GRA0125 , GDM2 , COE1
CHARACTER(LEN=80) :: ERR_MSG
INTEGER :: STATUS

DT2     = 2.D0*DT
DT0125  = 0.125D0*DT
DTDX    = DT*DX
DTDY    = DT*DY
DTDY2   = 2.D0*DTDY
DXDY    = DX*DY
DXDY4   = 4.D0*DXDY
DTDXDY  = DT*DX*DY
DTDXDDY = DT*DX/DY
DTDYDDX = DT*DY/DX
GRA05   = 0.5D0*GRA
GRA0125 = 0.125D0*GRA
GDM2    = GRA/MANNING**2
POW1    = 7.D0/3.D0
COE1    = 4.D0/3.D0

ML = LBOUND(PP1,2)
MU = UBOUND(PP1,2)
NL = LBOUND(PP1,1)
NU = UBOUND(PP1,1)

LPI0:DO I = I1 , I2 , IP
    ! CONTINUITY EQUATION
    JC = -1
    LPJ0:DO J = ML , MU , 1
        JC = JC + 2
        IF(NDD(I,J)>0) THEN
            AMT(1,JC) = 0.D0
            AMT(2,JC) = 0.D0
            AMT(3,JC) = 1.D0
            AMT(4,JC) = 0.D0
            AMT(5,JC) = 0.D0
            AMT(6,JC) = HH1(I,J)
            CYCLE LPJ0
        END IF
        IF(NDD(I,J)==-2) THEN
            AMT(1,JC) = 0.D0
            AMT(2,JC) = 0.D0
            AMT(3,JC) = 1.D0
            AMT(4,JC) = 0.D0
            AMT(5,JC) = 0.D0
            AMT(6,JC) = HH2(I,J)
            CYCLE LPJ0
        END IF
        HE = 0.5D0*(HH1(I,J)+HH1(I,J+1))
        HW = 0.5D0*(HH1(I,J)+HH1(I,J-1))
        HN = 0.5D0*(HH1(I,J)+HH1(I+1,J))
        HS = 0.5D0*(HH1(I,J)+HH1(I-1,J))
        AMT(1,JC) = 0.D0
        AMT(2,JC) = - DTDY*HW
        AMT(3,JC) = + DXDY4
        AMT(4,JC) = + DTDY*HE
        AMT(5,JC) = 0.D0
        AMT(6,JC) = + DXDY4*HH1(I,J) &
                  & - DTDX*(QQ2(I,J)*HN-QQ2(I-1,J)*HS+QQ1(I,J)*HN-QQ1(I-1,J)*HS) &
                  & - DTDY*(PP1(I,J)*HE-PP1(I,J-1)*HW)
    END DO LPJ0

    ! MOMENTUM EQUATION
    JJC = 0
    LPJ1:DO J = ML , MU-1 , 1
        JJC = JJC + 2 ! ID OF BAND MATRIX

        ! CLOSE BOUNDARY CONDITION
        IF(NDD(I,J)>0.OR.NDD(I,J+1)>0.OR.NDDP(I,J)>0) THEN
            AMT(1,JJC) = 0.D0
            AMT(2,JJC) = 0.D0
            AMT(3,JJC) = 1.D0
            AMT(4,JJC) = 0.D0
            AMT(5,JJC) = 0.D0
            AMT(6,JJC) = 0.D0
            CYCLE LPJ1
        END IF

        !  OPEN BOUNDARY CONDITION --- SURFACE LEVEL
        IF(NDD(I,J)==-2.OR.NDD(I,J+1)==-2) THEN
                GRV  = GRA05*(HH1(I,J)+HH1(I,J+1))
                AMT(1,JJC) = 0.D0
                AMT(2,JJC) = - GRV*DTDY
                AMT(3,JJC) = + DXDY*0.5D0*(HH1(I,J)+HH1(I,J+1))
                AMT(4,JJC) = + GRV*DTDY
                AMT(5,JJC) = 0.D0
                AMT(6,JJC) = + DXDY*PP1(I,J)*0.5D0*(HH1(I,J)+HH1(I,J+1)) &
                           & + GRV*DTDY*(DEP0(I,J+1)-DEP0(I,J))
            CYCLE LPJ1
        END IF

        ! FROUDE NUMBER
        FROUDE = DABS(PP1(I,J))/SQRT(GRA0125*(HH1(I,J+1)+HH1(I,J))**3)
        IF(FROUDE>=1.D0) THEN
            WF = 1.D0
        ELSE IF(FROUDE<1.D0.AND.FROUDE>0.25D0) THEN
            WF = COE1*(FROUDE-0.25D0)
        ELSE
            WF = 0.D0
        END IF
        IF(PP1(I,J)<0.D0) WF = -WF

        ! ADVECTION MOMENTUM
        ADVE = 0.25D0*((1.D0-WF)*PP1(I,J+1)+(1.D0+WF)*PP1(I,J)) &
                & / ((1.D0-DABS(WF))*HH1(I,J+1)+DMAX1(WF,0.D0)*HH1(I,J)-DMIN1(WF,0.D0)*HH1(I,J+2))
        ADVW = 0.25D0*((1.D0-WF)*PP1(I,J)+(1.D0+WF)*PP1(I,J-1)) &
                & / ((1.D0-DABS(WF))*HH1(I,J)+DMAX1(WF,0.D0)*HH1(I,J-1)-DMIN1(WF,0.D0)*HH1(I,J+1))

        ! CROSS MOMENTUM AND CROSS DIFFUSION
        IF(NDD(I+1,J)==0.AND.NDD(I+1,J+1)==0) THEN
            VN = 2.D0*(QQ2(I,J)+QQ2(I,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I+1,J)+HH1(I+1,J+1))
            PN = PP1(I+1,J)
        ELSE IF(NDD(I+1,J)/=0.AND.NDD(I+1,J+1)==0) THEN
            VN = 2.D0*(QQ2(I,J)+QQ2(I,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I,J)+HH1(I+1,J+1))
            PN = PP1(I,J)
        ELSE IF(NDD(I+1,J)==0.AND.NDD(I+1,J+1)/=0) THEN
            VN = 2.D0*(QQ2(I,J)+QQ2(I,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I+1,J)+HH1(I,J+1))
            PN = PP1(I,J)
        ELSE
            VN = 0.D0
            PN = PP1(I,J)
        END IF
        IF(NDD(I-1,J)==0.AND.NDD(I-1,J+1)==0) THEN
            VS = 2.D0*(QQ2(I-1,J)+QQ2(I-1,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I-1,J)+HH1(I-1,J+1))
            PS = PP1(I-1,J)
        ELSE IF(NDD(I-1,J)/=0.AND.NDD(I-1,J+1)==0) THEN
            VS = 2.D0*(QQ2(I-1,J)+QQ2(I-1,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I,J)+HH1(I-1,J+1))
            PS = PP1(I,J)
        ELSE IF(NDD(I-1,J)==0.AND.NDD(I-1,J+1)/=0) THEN
            VS = 2.D0*(QQ2(I-1,J)+QQ2(I-1,J+1))/(HH1(I,J)+HH1(I,J+1)+HH1(I-1,J)+HH1(I,J+1))
            PS = PP1(I,J)
        ELSE
            VS = 0.D0
            PS = PP1(I,J)
        END IF
        DIFFY = (VN+VS)**2*DT0125 + EDDYVITY

        IF(.NOT.UPS) THEN
            IF(NDD(I+1,J)<=0.AND.NDD(I+1,J+1)<=0) THEN
                CADL   = - 0.5D0*VS
                CADR   =   0.5D0*(-VN*PP2(I+1,J)+VS*PS-VN*PP1(I,J))
                CDIFFL =   DIFFY
                CDIFFR =   DIFFY*(PP2(I+1,J)-PP1(I,J)+PS)
            ELSE
                CADL   =   0.5D0*(VN-VS)
                CADR   =   0.5D0*(VS*PS-VN*PP1(I,J))
                CDIFFL =   0.D0
                CDIFFR =   DIFFY*(-PP1(I,J)+PS)
            END IF
        ELSEIF(UPS) THEN
            IF(NDD(I-1,J)<=0.AND.NDD(I-1,J+1)<=0) THEN
                CADL   =   0.5D0*VN
                CADR   =   0.5D0*(-VN*PN+VS*PP2(I-1,J)+VS*PP1(I,J))
                CDIFFL =   DIFFY
                CDIFFR =   DIFFY*(PN-PP1(I,J)+PP2(I-1,J))
            ELSE
                CADL =   0.5D0*(VN-VS)
                CADR =   0.5D0*(-VN*PN+VS*PP1(I,J))
                CDIFFL =   0.D0
                CDIFFR =   DIFFY*(PN-PP1(I,J))
            END IF
        END IF

        ! GRAVITY
        GRV  = GRA05*(HH1(I,J+1)+HH1(I,J))

        ! RESISTANCE TERM
        QAV = 0.125D0*(QQ1(I-1,J)+QQ1(I-1,J+1)+QQ2(I-1,J)+QQ2(I-1,J+1)+QQ1(I,J)+QQ1(I,J+1)+QQ2(I,J)+QQ2(I,J+1))
        IF(PP1(I,J)>=0.D0) THEN
            HAV = HH1(I,J)
        ELSE IF(PP1(I,J)<0.D0) THEN
            HAV = HH1(I,J+1)
        END IF
        RESIST = GDM2*SQRT(PP1(I,J)**2+QAV**2)/HAV**POW1

        ! X-DIFFUSION TERM
        DIFFX = (PP1(I,J)/(HH1(I,J+1)+HH1(I,J)))**2*DT2 + EDDYVITY

        AMT(1,JJC) = - (1.D0+WF)*ADVW*DTDY &
                   & - DIFFX*DTDYDDX

        AMT(2,JJC) = - GRV*DTDY

        AMT(3,JJC) = + DXDY &
                   & + ((1.D0+WF)*ADVE-(1.D0-WF)*ADVW)*DTDY &
                   & + CADL*DTDX &
                   & + RESIST*DTDXDY &
                   & + 2.D0*DIFFX*DTDYDDX &
                   & + CDIFFL*DTDXDDY

        AMT(4,JJC) = + GRV*DTDY

        AMT(5,JJC) = + (1.D0-WF)*ADVE*DTDY &
                   & - DIFFX*DTDYDDX

        AMT(6,JJC) = + DXDY*PP1(I,J) &
                   & + CADR*DTDX &
                   & + CDIFFR*DTDXDDY  &
                   & + GRV*DTDY*(DEP0(I,J+1)-DEP0(I,J))

    END DO LPJ1

    JC = MAX0(JC,JJC)
!    IF(I==59) THEN
!OPEN(5,FILE = 'AMT.TXT')
!DO J = 1 , JC
!WRITE(5,'(20E20.10)') (AMT(K,J),K=1,6)
!END DO
!CLOSE(5)
!WRITE(*,*) I
!    END IF
    CALL SOLVE5(JC,AMT)

!IF(I==59) THEN
!OPEN(5,FILE = 'AMT2.TXT')
!DO J = 1 , JC , 1
!WRITE(5,'(20E20.10)') (AMT(K,J),K=1,6)
!END DO
!CLOSE(5)
!PAUSE 'X'
!    END IF

    J0 = 0
    DO J = 2 , JC , 2
        J0 = J0 + 1
        PP2(I,J0) = AMT(6,J)
    END DO
    J0 = 0
    DO J = 1 , JC , 2
        J0 = J0 + 1
        HH2(I,J0) = AMT(6,J)
    END DO
END DO LPI0

RETURN
END SUBROUTINE CONTINUITY_MOMENTUM_X


! ============================================================================
! Name        : f03.f90
! Author      :
! Version     :
! Copyright   : Your copyright notice
! Description : Hello World in Fortran
! ============================================================================

program f03
    implicit none
    print *, "Hello MPI World"
end program