text:
C$Procedure   RDFFPL ( read triangular plate model flatfile )

      SUBROUTINE RDFFPL ( INFILE, PLTTYP, NV, VRTCES, NP, PLATES )
      IMPLICIT NONE

C$ Abstract
C
C     Read a triangular plate model's flatfile's data file.
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     MKPLAT User's Guide
C     SURFACE
C
C$ Keywords
C
C     PLATE
C     FLATFILE
C
C$ Declarations

      INCLUDE               'mkdsk02.inc'
      INCLUDE               'mkdsk.inc'

      CHARACTER*(*)         INFILE
      INTEGER               PLTTYP
      INTEGER               NV
      DOUBLE PRECISION      VRTCES ( 3, * )
      INTEGER               NP
      INTEGER               PLATES ( 3, * )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     INFILE     I   Name of plate data file.
C     PLTTYP     I   Integer code identifying type of data file.
C     NV         O   Number of vertices in model.
C     VRTCES     O   3 x NV array of vertices.
C     NP         O   Number of plates in model.
C     PLATES     O   3 x NP array of plates.
C
C
C$ Detailed_Input
C
C     INFILE     is the full pathname of the plate model flatfile.
C
C     PLTTYP     the integer code identifying the type of data
C                to read from INFILE. Two values are allowed:
C
C                    1    A standard plate-vertex data file
C                         containing vertex coordinates and
C                         the plate-vertex mappings.
C
C                    2    A Gaskel shape data file containing
C                         plate ordered vertex data.
C
C                    3    A vertex-facet table data file
C                         containing vertex coordinates and
C                         the facet listing.
C
C                    4    Rosetta/Osiris ".ver" format file
C                         containing vertex coordinates
C                         and a plate-vertex mapping.
C
C$ Detailed_Output
C
C     NV         is the number of vertices in the plate model.
C
C
C     VRTCES     is an array of the NV vertices given in a body-fixed
C                frame of reference. Elements
C
C                   VRTCES(J,I), J = 1 ... 3
C
C                are, respectively, the X, Y, and Z coordinates of the
C                Ith vertex.
C
C
C     PLATES     is an array of the NP plates. Elements
C
C                   PLATES(J,I), J = 1 ... 3
C
C                are, respectively, the indices of the vertices of the
C                Ith plate. The vertex indices range from 1 to NV.
C
C                The order of the vertices give the direction of the
C                plate's outward normal vector: the vector is parallel
C                to the cross product of the plate edge connecting
C                vertex 1 to vertex 2 and the plate edge connecting
C                vertex 2 to vertex 3, in that order.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If the number of vertices (NV) exceeds the value of MAXVRT,
C        SPICE(TOOMANYVERTICES) signals.
C
C     2) If the number of plates (NP) exceeds the value of MAXPLT,
C        SPICE(TOOMANYPLATES) signals.
C
C     3) If the NPARSD routine fails while parsing a double string,
C        SPICE(BADDOUBLEPRECISION) signals.
C
C     4) If the NPARSI routine fails while parsing an integer string,
C        SPICE(BADINTEGER) signals.
C
C     5) If a data line lacks the expected format (data type 3),
C        SPICE(BADDATALINE) signals.
C
C     6) If the PLTTYP variable has an uncoded value,
C        SPICE(BADDATATYPE) signals.
C
C     7) If the Gaskell ICQ parameter "Q" is too large, the error
C        SPICE(QPARAMOUTOFRANGE) is signaled.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     Input File Format
C     -----------------
C
C     INFILE is assumed to be structured as follows (values
C            are space delimited within each input record).
C
C     -Plate/vertex data file (PLATE_TYPE 1):
C
C        Line  1: NV
C        NV = the number of vertices (INTEGER)
C
C        Lines 2 to (NV+1): VID X Y Z
C        VID         = the vertex ID (INTEGER)
C        X, Y, and Z = the vertex coordinates (DOUBLE PRECISION)
C
C        Line  (NV+2): NP
C        NP = number of triangular plates (INTEGER)
C
C        Lines (NV+3) to (NV+3+NP+1): PID V1 V2 V3
C        PID            = the plate ID (INTEGER)
C        V1, V2, and V3 = the plate vertex ID's (INTEGER)
C
C     -Shape data file (PLATE_TYPE 2):
C
C        Line 1: Q
C        Q = number of vertex coordinates per side of the
C            cube face - not counting the origin. A
C            face has (Q+1)^2 vertex coordinates.
C
C        Lines 2 to (Q+1): X Y Z
C        X, Y, and Z = the vertex coordinates (DOUBLE PRECISION)
C
C     -Vertex/Facet table data file (PLATE_TYPE 3):
C
C        NV = the number of vertices (INTEGER)
C
C        Lines 1 to NV: 'V' X Y Z
C        'V'         = character flag indicating vertex data
C        X, Y, and Z = the vertex coordinates (DOUBLE PRECISION)
C
C        NP = number of triangular plates (INTEGER)
C
C        Lines (NV+1) to (NV+1+NP): 'F' V1 V2 V3
C        'F'            = character flag indicating facet
C                         (plate-vertex) data
C        V1, V2, and V3 = the plate vertex ID's (INTEGER)
C
C$ Examples
C
C  C     The following include file defines MAXVRT and
C  C     MAXPLT, the maximum number of vertices and plates
C  C     used in the plate model software.
C
C        INCLUDE               'mkdsk02.inc'
C
C        INTEGER               NV
C        INTEGER               NP
C        INTEGER               PLATES ( 3, MAXPLT )
C        INTEGER               PLTTYP
C
C        DOUBLE PRECISION      VRTCES ( 3, MAXVRT )
C
C        CHARACTER*(80)        INFILE
C
C        CALL PROMPT ( 'Enter name of plate data file : ', INFILE )
C        CALL PROMPT ( 'Enter type of plate data      : ', PLTTYP )
C
C        CALL RDFFPL ( INFILE,  PLTTYP, NV, VERTCES, NP, PLATES )
C
C$ Restrictions
C
C     It is the user's responsibility to properly dimension
C     arrays in the calling routine large enough to accept data
C     from the given file.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     J.A. Bytof      (JPL)
C
C$ Version
C
C-    MKDSK Version 5.1.0, 21-MAR-2017 (NJB)
C
C        Now references non-portable include file
C
C           mkdsk02.inc
C
C        instead of
C
C           dsk02.inc
C
C        This enables use of smaller maximum vertex and plate
C        counts on platforms where the default values consume
C        excessive memory.
C
C        Added checks for excessive vertex and plate counts.
C
C        Added FAILED checks following RDFFDI, RDNBL, and
C        PRSINT calls.
C
C-    MKDSK Version 5.0.0, 25-APR-2016 (NJB)
C
C        Bug fix: added check-out and return to input format 2
C        branch, in the case where a Q value out of range is found.
C
C
C
C        07-APR-2015 (NJB)
C
C           Increased string length used to parse numeric tokens.
C           Updated message indicating completion of input file read.
C           Cleaned up debugging code.
C
C        05-AUG-2014 (NJB)
C
C           Argument list change: arguments PID and VID have been
C           removed; arguments X, Y, Z and V1, V2, V3 have been
C           replaced by arrays VRTCES and PLATES respectively.
C
C           Local array VEC is no longer used.
C
C-    MKDSK Version 4.0.0, 08-JUN-2010 (NJB)
C
C        Added capability of reading Rosetta Osiris ".ver"
C        format file.
C
C-    MKDSK Version 3.1.0, 04-MAY-2010 (NJB)
C
C        Changed INCLUDE file from platmax.inc to dsk02.inc.
C        Added INCLUDE statement referencing mkdsk.inc to
C        declare parameter MAXQ.
C
C-    MKDSK Version 3.0.1, 08-OCT-2009 (NJB)
C
C        Re-ordered header sections.
C
C-    MKDSK Version 3.0.0, 25-OCT-2004 (EDW)
C
C        Added capability to process Bob Gaskell shape files.
C        Added MAXQ parameter to 'pltmax.inc'.
C
C-    MKDSK Version 2.0.0, 22-OCT-1998 (JAB)
C
C-    MKDSK Version 1.0.0, 09-APR-1997 (JAB)
C
C-&

C$ Index_Entries
C
C     read triangular plate model flatfile
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION      VNORM

      LOGICAL               EQSTR
      LOGICAL               EVEN
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               MAXWDS
      PARAMETER           ( MAXWDS = 3 )

      INTEGER               SHORT
      PARAMETER           ( SHORT = 80 )

      INTEGER               TOKLEN
      PARAMETER           ( TOKLEN = 40 )
C
C     Local variables
C
      CHARACTER*(TOKLEN)    ARC     ( 6 )
      CHARACTER*(160)       ERROR
      CHARACTER*(10)        FORMAT
      CHARACTER*(LNSIZE)    ILINE
      CHARACTER*(SHORT)     TOKENS ( MAXWDS )

      DOUBLE PRECISION      ARD    ( 3 )
      DOUBLE PRECISION      VEC    ( 3 )
      DOUBLE PRECISION      W1     ( 3 )
      DOUBLE PRECISION      W2     ( 3 )

      INTEGER               ARI    ( 4 )
      INTEGER               FACE
      INTEGER               I
      INTEGER               IX1
      INTEGER               IX2
      INTEGER               J
      INTEGER               N      ( 0:MAXQ, 0:MAXQ, 6)
      INTEGER               N0
      INTEGER               NC
      INTEGER               ND
      INTEGER               NI
      INTEGER               NREC
      INTEGER               NTOK
      INTEGER               PTR
      INTEGER               Q

      LOGICAL               EOF

C
C     Saved variables
C
      SAVE

C
C     Standard SPICE error handling.
C

      IF ( RETURN () ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'RDFFPL' )

C
C     Initialize record counter.
C
      NREC = 1

C
C     Read the data file corresponding to the PLTTYP ID.
C
      IF( PLTTYP .EQ. 1 ) THEN

C
C        Plate data type 1.
C
C        First line, a single integer, number of vertices.
C
         FORMAT = 'I'

         CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                 ARD, ARI, ARC, EOF  )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RDFFPL' )
            RETURN
         END IF

         NV = ARI(1)

C
C        Check the routine can handle the vertices.
C
         IF ( NV .GT. MAXVRT ) THEN
            CALL SETMSG ( 'Number of vertices # exceeds limit #.' )
            CALL ERRINT ( '#',  NV                                )
            CALL ERRINT ( '#',  MAXVRT                            )
            CALL SIGERR ( 'SPICE(TOOMANYVERTICES)'                )
            CALL CHKOUT ( 'RDFFPL'                                )
            RETURN
         END IF

C
C        Read data: integer double double double.
C
         FORMAT = 'I D D D'

         DO I = 1, NV

            NREC   = NREC + 1

            CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                    ARD, ARI, ARC, EOF  )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'RDFFPL' )
               RETURN
            END IF

            IF ( EOF ) THEN
               CALL EXTMSI ( 'End of file after line #.', '#', NREC )
            END IF

            VRTCES(1,I) = ARD(1)
            VRTCES(2,I) = ARD(2)
            VRTCES(3,I) = ARD(3)

         END DO

C
C        Plate data type 1 reads the plat-vertex information from
C        the plate file.
C
C        Read number of plates.
C
         NREC   = NREC + 1
         FORMAT = 'I'

         CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                 ARD, ARI, ARC, EOF  )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RDFFPL' )
            RETURN
         END IF

         IF ( EOF ) THEN
            CALL EXTMSI ( 'End of file after line #.', '#', NREC )
         END IF

C
C        Check we can process 'NP' plates.
C
         NP = ARI(1)

         IF ( NP .GT. MAXPLT ) THEN
            CALL SETMSG ( 'Number of plates # exceeds limit #.' )
            CALL ERRINT ( '#', NP                               )
            CALL ERRINT ( '#', MAXPLT                           )
            CALL SIGERR ( 'SPICE(TOOMANYPLATES)'                )
            CALL CHKOUT ( 'RDFFPL'                              )
            RETURN
         END IF

C
C        Read each plate's ID and corresponding vertex set defining
C        the plate.
C
         FORMAT = 'I I I I'

         DO I = 1, NP

            NREC   = NREC + 1

            CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                    ARD, ARI, ARC, EOF  )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'RDFFPL' )
               RETURN
            END IF

            IF ( EOF ) THEN
               CALL EXTMSI ( 'End of file after line #.', '#', NREC )
            END IF

            PLATES(1,I) = ARI(2)
            PLATES(2,I) = ARI(3)
            PLATES(3,I) = ARI(4)

         END DO

      ELSE IF( PLTTYP .EQ. 2 ) THEN

C
C        Plate data type 2.
C
C        Read a Gaskell shape file. Shape files contain only vertex
C        data, the plate-vertex mappings implicitly known from the
C        ordering of the data. The Gaskell model uses quadralaterals
C        plates, a plate defined by the vector set:
C
C        [ VEC( *, I, J  , FACE), VEC( *, I+1, J  , FACE),
C          VEC( *, I, J+1, FACE), VEC( *, I+1, J+1, FACE) ]
C
C        The file lists the vertex data in terms cube faces. Six
C        faces, each with Q+1 x Q+1 vertices.
C
         FORMAT = 'I'

         CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                 ARD, ARI, ARC, EOF  )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RDFFPL' )
            RETURN
         END IF

         Q  = ARI(1)
         NV = 6*( (Q+1)**2)

C
C        Does the read Q exceed the maximum value.
C
         IF ( Q .GT. MAXQ ) THEN

            CALL SETMSG ( 'Shape parameter Q = #, exceeds maximum '
     .      //            'value #. This error may indicate that '
     .      //            'the input file format is something other '
     .      //            'than the Gaskell ICQ format.'              )
            CALL ERRINT ( '#', Q                                      )
            CALL ERRINT ( '#', MAXQ                                   )
            CALL SIGERR ( 'SPICE(QPARAMOUTOFRANGE)'                   )
            CALL CHKOUT ( 'RDFFPL'                                    )
            RETURN

         END IF
C
C        Check the routine can handle the vertices.
C
         IF ( NV .GT. MAXVRT ) THEN
            CALL SETMSG ( 'Number of vertices # exceeds limit #.' )
            CALL ERRINT ( '#', NV                                 )
            CALL ERRINT ( '#', MAXVRT                             )
            CALL SIGERR ( 'SPICE(TOOMANYVERTICES)'                )
            CALL CHKOUT ( 'RDFFPL'                                )
            RETURN
         END IF

C
C        Read data: double double double
C
         FORMAT = 'D D D'

         N0 = 0

         DO FACE = 1, 6

            DO J = 0, Q

               DO I = 0, Q

                  NREC = NREC + 1

C
C                 Read the vertex data. The 3-Vector
C                    _                 _
C                   |  VEC(1,I,J,FACE)  |
C                   |  VEC(2,I,J,FACE)  |
C                   |_ VEC(3,I,J,FACE) _|
C
C                 contains the vertex coordinates.
C
                  CALL RDFFDI ( INFILE, NREC, FORMAT,
     .                          ND,     NI,   NC,
     .                          VEC,    ARI,  ARC,
     .                          EOF )

                  IF ( FAILED() ) THEN
                     CALL CHKOUT ( 'RDFFPL' )
                     RETURN
                  END IF

                  IF ( EOF ) THEN
                     CALL EXTMSI ( 'End of file after line #.',
     .                             '#', NREC)
                  END IF

C
C                Store the vertex ID to the 'N' array using the same
C                indexing as 'VEC'.
C
                 N0          = N0+1
                 N(I,J,FACE) = N0

C
C                Copy the vertex coordinates to the corresponding
C                plate model arrays. This operation is somewhat
C                repetitive, but clarity is most important.
C
                 VRTCES(1,N0) = VEC(1)
                 VRTCES(2,N0) = VEC(2)
                 VRTCES(3,N0) = VEC(3)

               END DO

            END DO

         END DO

C
C        Plate data type 2 calculates the plate-vertex mappings
C        from the ordering of the vertex data.

C
C        Number of plates based on the 'Q' value.
C
         NP = 12*(Q**2)

C
C        Check we can process 'NP' plates.
C
         IF ( NP .GT. MAXPLT ) THEN
            CALL SETMSG ( 'Number of plates # exceeds limit #.' )
            CALL ERRINT ( '#', NP                               )
            CALL ERRINT ( '#', MAXPLT                           )
            CALL SIGERR ( 'SPICE(TOOMANYPLATES)'                )
            CALL CHKOUT ( 'RDFFPL'                              )
            RETURN
         END IF

C
C        Expand a Gaskell shape model vertex ordering to
C        form a plate-vertex map.
C

C
C        Cube edges. Associate later IDs with
C        the first occurrence.
C
         DO I=1,Q-1

            N(I,Q,6) = N(Q-I, Q  , 4)
            N(I,0,6) = N(I  , Q  , 2)
            N(I,0,5) = N(Q  , Q-I, 1)
            N(I,0,4) = N(Q-I, 0  , 1)
            N(I,0,3) = N(0  , I  , 1)
            N(I,0,2) = N(I  , Q  , 1)

         END DO

         DO J=1,Q-1

            N(Q,J,6) = N(J  , Q, 5)
            N(Q,J,5) = N(0  , J, 4)
            N(Q,J,4) = N(0  , J, 3)
            N(Q,J,3) = N(0  , J, 2)
            N(0,J,6) = N(Q-J, Q, 3)
            N(0,J,5) = N(Q  , J, 2)

         END DO

C
C        Cube corners. Associate later IDs with
C        the first occurrence.
C
         N(0,0,3) = N(0,0,1)
         N(Q,0,4) = N(0,0,1)
         N(0,0,2) = N(0,Q,1)
         N(Q,0,3) = N(0,Q,1)
         N(0,0,4) = N(Q,0,1)
         N(Q,0,5) = N(Q,0,1)
         N(0,0,5) = N(Q,Q,1)
         N(Q,0,2) = N(Q,Q,1)
         N(0,0,6) = N(0,Q,2)
         N(Q,Q,3) = N(0,Q,2)
         N(0,Q,5) = N(Q,Q,2)
         N(Q,0,6) = N(Q,Q,2)
         N(Q,Q,4) = N(0,Q,3)
         N(0,Q,6) = N(0,Q,3)
         N(Q,Q,5) = N(0,Q,4)
         N(Q,Q,6) = N(0,Q,4)

         N0 = 0

         DO FACE = 1, 6

            DO I = 0, Q-1

               DO J = 0, Q-1

C
C                 As mentioned, the Gaskell shape model uses
C                 quadrilateral plates as opposed to triangular used by
C                 the NAIF plate system. We can reduce a quadrilateral
C                 into two triangles by adding a vector connecting
C                 opposite facing vertices creating two triangles
C                 (1 and 2):
C
C                               V1
C                               /\
C                              / 1\
C                          V4 /____\ V2
C                             \  2 /
C                              \  /
C                               \/
C                               V3
C
C                 or
C                                V1
C                               /|\
C                              / | \
C                         V4  /  |2 \ V2
C                             \ 1|  /
C                              \ | /
C                               \|/
C                                V3
C
C                 We chose the connecting vector to minimize the cross
C                 product magnitude of the connected vertex vectors.
C
C                    connection = MIN( V1 x V3, V2 X V4)
C

C
C                 Connection 1.
C

                  IX1 = (FACE-1)*((Q+1)**2) + ((Q+1)* J   ) +  I    + 1
                  IX2 = (FACE-1)*((Q+1)**2) + ((Q+1)*(J+1)) + (I+1) + 1

                  CALL VCRSS ( VRTCES(1,IX1), VRTCES(1,IX2), W1 )

C
C                 Connection 2.
C

                  IX1 = (FACE-1)*((Q+1)**2) + ((Q+1)* J   ) + (I+1) + 1
                  IX2 = (FACE-1)*((Q+1)**2) + ((Q+1)*(J+1)) +  I    + 1

                  CALL VCRSS ( VRTCES(1,IX1), VRTCES(1,IX2), W2 )

C
C                 Calculate the magnitudes of the cross products;
C                 branch based on the minimum value.
C
C                 Fill two entries of the plate-vertex mapping
C                 from the divided quadralateral.
C
                  IF( VNORM(W1) .LE. VNORM(W2) ) THEN
C
C                    Apply connection 1.
C
                     N0           = N0+1
                     PLATES(1,N0) = N(I ,  J,   FACE)
                     PLATES(2,N0) = N(I+1, J+1, FACE)
                     PLATES(3,N0) = N(I+1, J,   FACE)

                     N0           = N0+1
                     PLATES(1,N0) = N(I,   J,   FACE)
                     PLATES(2,N0) = N(I,   J+1, FACE)
                     PLATES(3,N0) = N(I+1, J+1, FACE)

                  ELSE
C
C                    Apply connection 2.
C
                     N0           = N0+1
                     PLATES(1,N0) = N(I,   J,   FACE)
                     PLATES(2,N0) = N(I,   J+1, FACE)
                     PLATES(3,N0) = N(I+1, J,   FACE)

                     N0           = N0+1
                     PLATES(1,N0) = N(I+1, J,   FACE)
                     PLATES(2,N0) = N(I,   J+1, FACE)
                     PLATES(3,N0) = N(I+1, J+1, FACE)

                  END IF

               END DO

            END DO

         END DO

      ELSE IF( PLTTYP .EQ. 3 ) THEN
C
C        Plate data type 3.
C
C        Read data. The data format depends on the
C        prefix marker either 'v' (vertex) or 'f' (facet).
C        Since we don't know the prefix a priori,
C        extract the four data elements, then parse.
C
         NP     = 0
         NV     = 0
         FORMAT = 'C C C C'
         EOF    = .FALSE.

         DO WHILE ( .NOT. EOF )

            CALL RDFFDI ( INFILE, NREC, FORMAT, ND, NI, NC,
     .                    ARD, ARI, ARC, EOF  )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'RDFFPL' )
               RETURN
            END IF

            IF ( .NOT. EOF ) THEN

C
C              Branch based on the type of data, 'v' or 'f'.
C
               IF(  EQSTR( ARC(1), 'v' )  ) THEN

C
C                 Vertex data. Process a line of three doubles.
C                 The algorithm includes an implicit assumption that
C                 the data order defines the vertex IDs.  If not,
C                 the output plate files will be useless.
C
                  NV = NV + 1

                  IF ( NV .GT. MAXVRT ) THEN

                     CALL SETMSG ( 'Number of vertices # exceeds '
     .               //            'limit #.'                    )
                     CALL ERRINT ( '#',  NV                      )
                     CALL ERRINT ( '#',  MAXVRT                  )
                     CALL SIGERR ( 'SPICE(TOOMANYVERTICES)'      )
                     CALL CHKOUT ( 'RDFFPL'                      )
                     RETURN

                  END IF

                  DO J = 2, 4

                     CALL NPARSD ( ARC(J), ARD(J-1), ERROR, PTR )

                     IF ( PTR .NE. 0 ) THEN
                        CALL SETMSG ( 'D.P. error (#) in line #.' )
                        CALL ERRCH  ( '#', ERROR                  )
                        CALL ERRINT ( '#', NREC                   )
                        CALL SIGERR ( 'SPICE(BADDOUBLEPRECISION)' )
                        CALL CHKOUT ( 'RDFFPL'                    )
                        RETURN
                     END IF

                  END DO

                 VRTCES(1,NV) = ARD(1)
                 VRTCES(2,NV) = ARD(2)
                 VRTCES(3,NV) = ARD(3)

               ELSE IF ( EQSTR( ARC(1), 'f') ) THEN

C
C                 Plate-vertex data. Process a line of three integers.
C                 The algorithm includes an implicit assumption that
C                 the data order defines the plate IDs. If not,
C                 the output plate files will be useless.
C
                  NP = NP + 1

                 IF ( NP .GT. MAXPLT ) THEN

                    CALL SETMSG( 'Number of plates # exceeds limit #.' )
                    CALL ERRINT( '#', NP                               )
                    CALL ERRINT( '#', MAXPLT                           )
                    CALL SIGERR( 'SPICE(TOOMANYPLATES)'                )
                    CALL CHKOUT( 'RDFFPL'                              )
                    RETURN

                 END IF

                  DO J =2, 4

                     CALL NPARSI ( ARC(J), ARI(J-1), ERROR, PTR )

                     IF ( PTR .NE. 0 ) THEN
                        CALL SETMSG ( 'Integer error (#) in line #.' )
                        CALL ERRCH  ( '#', ERROR                     )
                        CALL ERRINT ( '#', NREC                      )
                        CALL SIGERR ( 'SPICE(BADINTEGER)'            )
                        CALL CHKOUT ( 'RDFFPL'                       )
                        RETURN
                     END IF

                  END DO

                  PLATES(1,NP) = ARI(1)
                  PLATES(2,NP) = ARI(2)
                  PLATES(3,NP) = ARI(3)

               ELSE

C
C                 This block executes if the first characters of
C                 a type 3 data lacks a 'v' or 'f'.
C
                  CALL SETMSG ( 'Bad data line at record #.' )
                  CALL ERRINT ( '#', NREC                    )
                  CALL SIGERR ( 'SPICE(BADDATALINE)'         )
                  CALL CHKOUT ( 'RDFFPL'                     )
                  RETURN

               END IF

               NREC = NREC + 1

            END IF

         END DO

C
C        All type 3 data read. We should have a non zero NP and NV,
C        if not, something failed.
C
         IF ( (NP .EQ. 0) .OR. (NV .EQ. 0) ) THEN
            CALL SETMSG ( 'Read error, type 3 data file.'
     .                 // ' Num plates = #1, num verts = #2.'
     .                 // 'Both should be non-zero'         )
            CALL ERRINT ( '#1', NP                          )
            CALL ERRINT ( '#2', NV                          )
            CALL SIGERR ( 'SPICE(DATAREADFAILED)'           )
            CALL CHKOUT ( 'RDFFPL'                          )
            RETURN
         END IF

      ELSE IF ( PLTTYP .EQ. 4 ) THEN
C
C        We have a Rosetta Osiris style ".ver" file.
C
C        Get the vertex and plate count from the first non-blank
C        input line.
C
         CALL RDNBL ( INFILE, ILINE, EOF )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RDFFPL' )
            RETURN
         END IF

         CALL LPARSM ( ILINE, ' ,', 2, NTOK, TOKENS )

         IF ( NTOK .NE. 2 ) THEN

            CALL SETMSG ( 'Vertex and plate count were expected '
     .      //            'on first line of input file. Line was '
     .      //            '#.'                                    )
            CALL ERRCH  ( '#', ILINE                              )
            CALL SIGERR ( 'SPICE(UNRECOGNIZEDFORMAT)'             )
            CALL CHKOUT ( 'RDFFPL'                                )
            RETURN

         END IF

         CALL PRSINT ( TOKENS(1), NV )
         CALL PRSINT ( TOKENS(2), NP )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'RDFFPL' )
            RETURN
         END IF

         IF ( NV .GT. MAXVRT ) THEN

            CALL SETMSG ( 'Number of vertices # exceeds limit #.' )
            CALL ERRINT ( '#',  NV                                )
            CALL ERRINT ( '#',  MAXVRT                            )
            CALL SIGERR ( 'SPICE(TOOMANYVERTICES)'                )
            CALL CHKOUT ( 'RDFFPL'                                )
            RETURN

         END IF

         IF ( NP .GT. MAXPLT ) THEN

            CALL SETMSG ( 'Number of plates # exceeds limit #.' )
            CALL ERRINT ( '#', NP                               )
            CALL ERRINT ( '#', MAXPLT                           )
            CALL SIGERR ( 'SPICE(TOOMANYPLATES)'                )
            CALL CHKOUT ( 'RDFFPL'                              )
            RETURN

         END IF

C
C        Read the vertex data and store it in the output vertex array.
C
         DO I = 1, NV

            CALL RDNBL ( INFILE, ILINE, EOF )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'RDFFPL' )
               RETURN
            END IF

            IF ( EOF ) THEN

               CALL SETMSG ( 'Expected to find # vertices in input '
     .         //            'file # but ran out of data after reading '
     .         //            '# lines of vertex data.'                 )
               CALL ERRINT ( '#', NV                                   )
               CALL ERRCH  ( '#', INFILE                               )
               CALL ERRINT ( '#', I                                    )
               CALL SIGERR ( 'SPICE(FILETRUNCATED)'                    )
               CALL CHKOUT ( 'RDFFPL'                                  )
               RETURN

            END IF

C
C           Parse the input line; we expect it to contain the
C           components of the Ith vertex.
C
            CALL LPARSM ( ILINE, ' ,', 3, NTOK, TOKENS )

            IF ( NTOK .NE. 3 ) THEN

               CALL SETMSG ( 'Three vertex components were expected '
     .         //            'on current line of input file. Line was '
     .         //            '#.'                                    )
               CALL ERRCH  ( '#', ILINE                              )
               CALL SIGERR ( 'SPICE(UNRECOGNIZEDFORMAT)'             )
               CALL CHKOUT ( 'RDFFPL'                                )
               RETURN

            END IF

            CALL PRSDP ( TOKENS(1), VRTCES(1,I) )
            CALL PRSDP ( TOKENS(2), VRTCES(2,I) )
            CALL PRSDP ( TOKENS(3), VRTCES(3,I) )

            IF ( FAILED() ) THEN
               CALL CHKOUT( 'RDFFPL' )
               RETURN
            END IF

         END DO

C
C        Read the plate indices. Discard the first
C        and every second non-blank line of plate data; these
C        lines all contain the same vertex count (3).
C
         J = 0

         DO I = 1,  2 * NP

            CALL RDNBL ( INFILE, ILINE, EOF )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'RDFFPL' )
               RETURN
            END IF

            IF ( EOF ) THEN

               CALL SETMSG ( 'Expected to find # plates in input '
     .         //            'file # but ran out of data after reading '
     .         //            '# lines of plate data.'                  )
               CALL ERRINT ( '#', NP                                   )
               CALL ERRCH  ( '#', INFILE                               )
               CALL ERRINT ( '#', I                                    )
               CALL SIGERR ( 'SPICE(FILETRUNCATED)'                    )
               CALL CHKOUT ( 'RDFFPL'                                  )
               RETURN

            END IF

            IF ( EVEN(I) ) THEN
C
C              This line should contain the components of the Jth plate.
C
               J = J + 1

C
C              Parse the current line and store the plate components.
C
               CALL LPARSM ( ILINE, ' ,', 3, NTOK, TOKENS )

               IF ( NTOK .NE. 3 ) THEN

                  CALL SETMSG ( 'Three plate components were expected '
     .            //            'on current line of input file. Line '
     .            //            'was #.'                              )
                  CALL ERRCH  ( '#', ILINE                            )
                  CALL SIGERR ( 'SPICE(UNRECOGNIZEDFORMAT)'           )
                  CALL CHKOUT ( 'RDFFPL'                              )
                  RETURN

               END IF

               CALL PRSINT ( TOKENS(1), PLATES(1,J) )
               CALL PRSINT ( TOKENS(2), PLATES(2,J) )
               CALL PRSINT ( TOKENS(3), PLATES(3,J) )

               IF ( FAILED() ) THEN
                  CALL CHKOUT( 'RDFFPL' )
                  RETURN
               END IF

            END IF

         END DO

      ELSE

         CALL SETMSG ( 'Unkown plate data type: #.' )
         CALL ERRINT( '#', PLTTYP                   )
         CALL SIGERR ( 'SPICE(BADDATATYPE)'         )
         CALL CHKOUT ( 'RDFFPL'                     )
         RETURN

      END IF

      CALL TOSTDO ( '...Done reading plate model input file.' )
      CALL TOSTDO ( ' ' )

      CALL CHKOUT ( 'RDFFPL' )

      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:
c iritest.for, version number can be found at the end of this comment.
c-----------------------------------------------------------------------
c
c test program for the iri_web subroutine
c
c-version-mm/dd/yy ----------corrections--------------------------
c 2000.01 05/07/01 initial version
c 2000.02 07/11/01 line 210: do i=1,100 instead of i=2,100 (K. Tokar)
c 2000.03 28/12/01 output oar(39) for IG12 (R. Conkright, NGDC, NOAA)
c 2000.04 28/10/02 replace TAB/6 blanks, enforce 72/line (D. Simpson)
c 2000.05 02/06/03 Ne(Te) only 300,400; foF1 and hmF1 output corr.
c 2000.06 01/19/05 (.not.jf(20)) instead of (..jf(2)) (G. Schiralli)
c 2005.01 05/06/06 included spread-F (jf(28)) and topside (jf(29)) options
C 2007.00 05/18/07 Release of IRI-2007
c 2007.02 10/31/08 outf(100) -> outf(500), numhei=numstp=500
c 2007.03 02/12/09 added new D-region option (h=-3)
c 2007.11 04/19/10 correct TEC for normal output  [Shunrong Zhang]
c
C 2012.00 10/05/11 IRI-2012: bottomside B0 B1 model (SHAMDB0D, SHAB1D),
C 2012.00 10/05/11    bottomside Ni model (iriflip.for), auroral foE
C 2012.00 10/05/11    storm model (storme_ap), Te with PF10.7 (elteik),
C 2012.00 10/05/11    oval kp model (auroral_boundary), IGRF-11(igrf.for),
C 2012.00 10/05/11    NRLMSIS00 (cira.for), CGM coordinates, F10.7 daily
C 2012.00 10/05/11    81-day 365-day indices (apf107.dat), ap->kp (ckp),
C 2012.00 10/05/11    array size change jf(50) outf(20,1000), oarr(100).
C 2012.00 03/21/12    PIKTAB=4 output for D-region
C 2012.01 09/16/12 Corrected UT output (hour-25)
C 2012.03 09/18/14 input JF(18): FIELDG not UT_LT          (A. Mazzella)
C 2012.03 09/18/14 rzino, igino logical not real; and more (A. Mazzella)
C 2012.03 09/18/14 jf defaults and special for piktab=3    (A. Mazzella)
C 2012.03 09/23/14 jf(22): output option Ni in [m-3]/1.e9
C 2012.03 09/23/14 jf(36:38) include in input choices
C 2012.03 09/23/14 added output options for parameters oar(59:86)
C 2012.05 04/27/15 delete double line CHARACTER*9  pname(6)
C 2015.01 06/30/15 scid=1.0E-9                              (A. Charisi)
C 2015.01 06/30/15 jino-outf(9)&jio2-outf(8) write.jio2,jino(A. Charisi)
C 2015.01 07/12/15 output of ion density/composition now with 3 digits
C 2015.01 07/12/15 call read_ig_rz readapf107
C 2015.02 08/13/15 delete COMMON/CONST2
C 2015.02 08/13/15 ursifo=jf(5) just before IRI_SUB call
C 2016.01 06/01/16 User specified B0 when jf(43)=false
C 2016.02 08/22/16 Allow user input options for F10.7D and F10.7_81
C
      Program IRItest

      INTEGER           pad1(6),jdprof(77),piktab
      DIMENSION         outf(20,1000),oar(100,1000),jfi(6)
      LOGICAL		    jf(50),rzino,igino
      CHARACTER*2       timev(2)
      CHARACTER*3       uni(86),sopt,seopt
      CHARACTER*4       IMZ(8),MAP,xtex,coorv(2)
      CHARACTER*5       ITEXT(8),tsopt
      CHARACTER*6       pna(86)
      CHARACTER*7       popt
      CHARACTER*8       bopt,topt
      CHARACTER*9       pname(7)
      CHARACTER*10      dopt,hopt
      CHARACTER*11      iopt,rzopt,igopt,f8opt,fdopt
      CHARACTER*16      f1opt

      DATA  IMZ  /' km ','GEOD','GEOD','yyyy',' mm ',' dd ','YEAR',
     &      'L.T.'/, ITEXT/'  H  ',' LATI',
     &      ' LONG',' YEAR','MONTH',' DAY ','DAYOF',' HOUR'/

      data pna/'NmF2','hmF2','NmF1','hmF1','NmE','hmE','NmD','hmD',
     &  'h05','B0','NVmin','hVtop','Tpeak','hTpek','T300','T400','T600',
     &  'T1400','T3000','T120','Ti450','hTeTi','sza','sundec','dip',
     &  'diplat','modip','Lati','Srise','Sset','season','Longi',
     &  'Rz12','cov','B1','M3000','TEC','TECtop','IG12','F1prb','F107d',
     &  'C1','daynr','vdrft','foF2r','F10781','foEr','sprd_F','MLAT',
     &  'MLON','Ap_t','Ap_d','invdip','MLTinv','CGMlat','CGMlon',
     &  'CGMmlt','CGM_AB','CGMm0','CGMm1','CGMm2','CGMm3','CGMm4',
     &  'CGMm5','CGMm6','CGMm7','CGMm8','CGMm9','CGMm10','CGMm11',
     &  'CGMm12','CGMm13','CGMm14','CGMm15','CGMm16','CGMm17','CGMm18',
     &  'CGMm19','CGMm20','CGMm21','CGMm22','CGMm23','kp_t','dec','L',
     &  'DIMO'/
      data uni/'m-3','km','m-3','km','m-3','km','m-3','km','km','km',
     &   'm-3','km','K','km',7*'K','km',6*'deg',2*'h',' ','deg',4*' ',
     &   'm-2','%',5*' ','m/s',4*' ',2*'deg',2*' ','deg','h',2*'deg',
     &   'h',25*'deg',' ','deg',' ','Gau'/,
     &   timev/'LT','UT'/,coorv/'geog','geom'/

      data jfi/8,9,13,14,15,16/

c	  COMMON/const2/icalls,montho,nmono,iyearo,idaynro,ursifo,rzino,
c     &	            igino,ut0
c
c		icalls=0
c		montho=-1
c		nmono=-1
c		iyearo=-1
c		idaynro=-1
c		rzino=.true.
c		igino=.true.
c		ut0=-1

		call read_ig_rz
        call readapf107

        nummax=1000

        do 6249 i=1,100
6249    oar(i,1)=-1.0

c user input of IRI input parameters
c
1       print *,'jmag(=0/1,geog/geom),lati/deg,long/deg'
        read(5,*) jm,xlat,xlon
        print *,'year(yyyy),mmdd(or -ddd),iut(=0/1,LT/UT),hour'
        read(5,*) iy,imd,iut,hour
        print *,'height/km'
        read(5,*) hx

        print *,'output-option'
        print *,'(enter 0 for standard table of IRI parameters)'
        print *,'(enter 1 for list of peak heights and densities)'
        print *,'(enter 2 for plasma frequencies, B0, M3000, ',
     &                        'valley, width and depth,)'
        print *,'(enter 3 for 6 parameters of your choice)'
        print *,'(enter 4 for D-region models at 60,65,..,110 km)'
        print *,'(enter 5 special test output)'
        read(5,*) piktab

        print *,'upper height [km] for TEC integration (0 for no TEC)'
        read(5,*) htec_max

        print *,'variable? (1/2/../8 for height/lat/long/year/month/',
     &                        'day/day of year/hour)'
        read(5,*) ivar
        print *,'begin, end, and stepsize for the selected variable'
        read(5,*) vbeg,vend,vstp

        print *,'Options: t(rue) or f(alse)'
        print *,'Enter 0 to use standard or 1 to enter your own'
        read(5,*) jchoice
          do i=1,50
                jf(i)=.true.
                enddo
          if(piktab.eq.4) jf(24)=.false.
        if(jchoice.eq.0) then
c defaults for jf(1:50)
c          jf(1)=.false.      ! f=no electron densities
c          jf(2)=.false.      ! f=no temperatures
c          jf(3)=.false.      ! f=no ion composition
          jf(4)=.false.      ! t=B0table f=other models (f)
          jf(5)=.false.      ! t=CCIR  f=URSI foF2 model (f)
          jf(6)=.false.      ! t=DS95+DY85   f=RBV10+TTS03 (f)
c          jf(7)=.false.      ! t=tops f10.7<188 f=unlimited (t)
          jf(21)=.false.     ! f=ion drift not computed (f)
c          jf(22)=.false.     ! ion densities in m-3 (t)
          jf(23)=.false.     ! t=AEROS/ISIS f=TTS Te with PF10.7
c          jf(24)=.false.     ! t=D-reg-IRI-1990 f=FT-2001 (t)
c          jf(26)=.false.	  ! f=STORM model turned off (t)
          jf(28)=.false.	  ! f=spread-F not computed (f)
          jf(29)=.false.     ! t=old  f=New Topside options (f)
          jf(30)=.false.     ! t=corr f=NeQuick topside (f)
c          jf(31)=.false.     ! t=B0ABT f=Gulyaeva (t)
          jf(33)=.false. 	  ! f=auroral boundary off (f)
c          jf(34)=.false. 	  ! t=messages on
          jf(35)=.false. 	  ! f=auroral E-storm model off
c          jf(36)=.false. 	  ! t=hmF2 w/out foF2_storm f=with
c          jf(37)=.false. 	  ! t=topside w/out foF2_storm f=with
c          jf(38)=.false. 	  ! t=WRITEs off in IRIFLIP f=on
          jf(39)=.false. 	  ! new hmF2 models
c          jf(40)=.false. 	  ! t=AMTB-model, f=Shubin-COSMIC model
c          jf(41)=.false. 	  ! COV=f(IG12) (IRI before Oct 2015)
c          jf(42)=.false. 	  ! Te w/o PF10.7 dependance
c          jf(43)=.false. 	  ! B0 user input
        else
          print *,'Compute Ne, T, Ni? (enter: t,t,t  if you want all)'
          read(5,*) jf(1),jf(2),jf(3)
        if(jf(1)) then
              print *,'LAY version: t=standard ver., f=LAY version.',
     &              ' {standard:t}'
              read(5,*) jf(11)
              print *,'Ne Topside: t=IRI-2001/h0.5, f=new options {f}'
              read(5,*) jf(29)
              print *,'Ne Topside: t=IRI01_corrt, f=NeQuick/h0.5 {f}'
              read(5,*) jf(30)
              print *,'Ne Topside: t=F10.7<188, f=unlimited {t}'
              read(5,*) jf(7)
              print *,'Ne topside: t=w/o foF2 storm model, f=with {t}'
              read(5,*) jf(37)
              print *,'foF2 model: t=CCIR, f=URSI-88 {f}'
              read(5,*) jf(5)
              print *,'foF2: t=with storm model, f=without {t}'
              read(5,*) jf(26)
              print *,'hmF2: t=f(M3000F2), f=new models {f}'
              read(5,*) jf(39)
              print *,'hmF2: t=AMTB-model, f=Shubin-COSMIC model {t}'
              read(5,*) jf(40)
              if(jf(39)) then
                 print *,'hmF2: t=w/o foF2 storm model, f=with {t}'
                 read(5,*) jf(36)
                 endif
              print *,'F2 peak density or foF2: t=model, ',
     &              'f=user input {t}'
              read(5,*) jf(8)
              print *,'F2 peak height or M3000F2: t=model, ',
     &              'f=user input {t}'
              read(5,*) jf(9)
              print *,'Auroral boundary model: t=on, f=off {f}'
              read(5,*) jf(33)
              print *,'Bottomside thickness B0: t=Bil-2000, ',
     &            'f=other options {f}.'
              read(5,*) jf(4)
              print *,'Bottomside thickness B0: t=ABT-2009, ',
     &            'f= Gul-1987 {t}.'
              read(5,*) jf(31)
              print *,'Bottomside thickness B0: t=model, ',
     &            'f=user input {t}.'
              read(5,*) jf(43)
              print *,'F1 peak density or foF1: t=model, ',
     &            'f=user input {t}'
              read(5,*) jf(13)
            if(.not.jf(11)) then
              print *,'F1 peak height: t=model, f=user input {t}'
              read(5,*) jf(14)
            endif
              print *,'F1: t=with probability model, f=without   {t}'
              read(5,*) jf(19)
              print *,'F1: t=standard probability, f=with L ',
     &            'condition {t}'
              read(5,*) jf(20)
              print *,'E peak density or foE: t=model, f=user input {t}'
              read(5,*) jf(15)
              print *,'E peak height: t=model, f=user input {t}'
              read(5,*) jf(16)
              print *,'E peak auroral storm model: t=on, f=off {f}'
              read(5,*) jf(35)
              print *,'D: t=IRI-1990, f= FT-2001 {t}'
              read(5,*) jf(24)
          endif
        if(jf(2)) then
              print *,'Te(Ne) model: t=not used, f=correlation is',
     &          ' used. {t}'
              read(5,*) jf(10)
              print *,'Te: t=Bil-1985, f=TBT-2012 {f}'
              read(5,*) jf(23)
              print *,'Te: t=TBT-2012 with PF107 dep., f=w/o {t}'
              read(5,*) jf(42)
          endif
        if(jf(3)) then
              print *,'Ion comp. model: t=DS95/DY85, f=RBV10/TTS05 {f}'
              read(5,*) jf(6)
              if(.not.jf(6)) then
              	print *,'IRIFLIP: t=messages off, f=on {t}'
              	read(5,*) jf(38)
              	endif
              print *,'Ni: t=ion composition in %, f=ion densities',
     &             'in cm-3 {t}'
              read(5,*) jf(22)
              endif
           print *,'Equat. Vert. Ion Drift: t=computed, ',
     &            'f=not computed {t}'
              read(5,*) jf(21)
           print *,'Spread-F probability: t=computed, ',
     &            'f=not computed {t}'
              read(5,*) jf(28)
           print *,'COV: t: COV=F10.7_365, f: COV=func(IG12).  {t}'
              read(5,*) jf(41)
           print *,'Sunspot index: t=from file, f=user input.  {t}'
              read(5,*) jf(17)
           print *,'Ionospheric index: t=from file, f=user input. {t}'
              read(5,*) jf(27)
           print *,'F10.7D Index: t=from file, f=user input {t}'
              read(5,*) jf(25)
           print *,'F10.7_81 Index: t=from file, f=user input {t}'
              read(5,*) jf(32)
           print *,'dip, magbr, modip: t=IGRF, f=old FIELDG using ',
     &             'POGO68/10 for 1973 {t}'
              read(5,*) jf(18)
           print *,'Messages on (t) off (f) {t}'
              read(5,*) jf(34)
           print *,'Message output unit: t=(UNIT=6), f=(UNIT=11). {t}'
              read(5,*) jf(12)
       endif

c       if(piktab.gt.3) jf(24)=.false.
c option to enter six additional parameters
c

      if(PIKTAB.eq.3) then
        print *,'6 Parameters of your choice (number:1-86)'
        print *,(pna(j),j=1,10)
        print *,(pna(j),j=11,20)
        print *,(pna(j),j=21,30)
        print *,(pna(j),j=31,40)
        print *,(pna(j),j=41,50)
        print *,(pna(j),j=51,60)
        print *,(pna(j),j=61,70)
        print *,(pna(j),j=71,80)
        print *,(pna(j),j=81,86)
        print *,'or 0,0,0,0,0,0 for default:'
        print *,'      spread-F probability [48]'
        print *,'      equatorial vertical ion drift [44]'
        print *,'      foF2_storm/foF2_quiet [45]'
        print *,'      foE_storm/foE_quiet [47]'
        print *,'      eqward auroral boundy CGM-Lat [58]'
c        print *,'      solar zenith angle [23]'
c        print *,'      modified dip latitude [27]'
        print *,'      Ap for current UT [51]'
        read(5,*) (pad1(j),j=1,6)
c change defaults for computation of specific parameters
        	jf(21)=.true.  ! spread-F prob. computed
        	jf(28)=.true.  ! vertical ion drift computed
        	jf(33)=.true.  ! auroral boundary computed
        	jf(35)=.true.  ! foE_storm computed
        if(pad1(1).eq.0) then
        	pad1(1)=48     ! spread-F probability
            pad1(2)=44     ! equatorial vertical ion drift
            pad1(3)=45     ! fof2_storm/foF2_quiet
            pad1(4)=47     ! foE_storm/foE_quiet
            pad1(5)=58     ! CGM_lat auroral boundary
            pad1(6)=51     ! ap for current UT
            endif
      endif

c option to enter measured values for NmF2, hmF2, NmF1, hmF1, NmE, hmE,
c B0, N(300), N(400), N(600) if available;
c
          print *,' '
          print *,' '
          print *,' '
          numstp=int((vend-vbeg)/vstp)+1
              if(ivar.eq.1) numstp=1
       if(jf(1)) then
         if(.not.jf(8).or..not.jf(9).or..not.jf(13).or..not.jf(14).or.
     &      .not.jf(15).or..not.jf(16).or..not.jf(43)) then
            var=vbeg
            i=1
2234        if(.not.jf(8)) then
              jf(26)=.false.    ! storm model off, if user input
              print *,'foF2/Mhz or NmF2/m-3 for ',itext(ivar),
     &             '=',var
              read(5,*) oar(1,i)
              pname(1)='foF2/MHz'
              if(oar(1,i).gt.30.) pname(1)='NmF2/m-3'
              endif
            if(.not.jf(9)) then
              print *,'hmF2/km or M3000F2 for ',itext(ivar),
     &              '=',var
              read(5,*) oar(2,i)
              pname(2)='M(3000)F2'
              if(oar(2,i).gt.50.) pname(2)='hmF2/km'
              endif
            if(.not.jf(13)) then
              print *,'foF1/MHz or NmF1/m-3 for ',itext(ivar),
     &               '=',var
              read(5,*) oar(3,i)
              pname(3)='foF1/MHz'
              if(oar(3,i).gt.30.) pname(3)='NmF1/m-3'
              endif
            if(.not.jf(14)) then
              print *,'hmF1/km for ',itext(ivar),'=',var
              read(5,*) oar(4,i)
              pname(4)='hmF1/km'
              endif
            if(.not.jf(15)) then
              print *,'foE/MHz or NmE/m-3 for ',itext(ivar),
     &                '=',var
              read(5,*) oar(5,i)
              pname(5)='foE/MHz'
              if(oar(5,i).gt.30.) pname(5)='NmE/m-3'
              endif
            if(.not.jf(16)) then
              print *,'hmE/km for ',itext(ivar),'=',var
              read(5,*) oar(6,i)
              pname(6)='hmE/km'
              endif
            if(.not.jf(43)) then
              print *,'B0/km for ',itext(ivar),'=',var
              read(5,*) oar(10,i)
              pname(7)='B0/km '
              endif
            i=i+1
            var=var+vstp
            if(ivar.gt.1.and.var.le.vend) goto 2234
         endif
       endif

c option to enter Ne for Te-Ne relationship
c
        if(jf(2).and..not.jf(10)) then
          var=vbeg
          do 1235 i=1,numstp
            print *,'Ne(300km),Ne(400km)/m-3',
     &         ' for ',itext(ivar),'=',var,' [-1 if not]'
            read(5,*) oar(15,i),oar(16,i)
1235        var=var+vstp
          endif

c option to enter F107D and/or PF107
c
            if(.not.jf(25)) then
                        print *,'User input for F107D:'
                        read(5,*) f107d
                        do i=1,100
                                    oar(41,i)=f107d
                                    enddo
                        endif

            if(.not.jf(32)) then
                        print *,'User input for PF107:'
                        read(5,*) pf107d
                        do i=1,100
                                    oar(46,i)=pf107d
                                    enddo
                        endif

c option to enter Rz12 and/or IG12
c
            if(.not.jf(17)) then
                        print *,'User input for Rz12'
                        read(5,*) oar(33,1)
                        do i=2,100
                                    oar(33,i)=oar(33,1)
                                    enddo
                        endif

            if(.not.jf(27)) then
                        print *,'User input for IG12'
                        read(5,*) oar(39,1)
                        do i=2,100
                                    oar(39,i)=oar(39,1)
                                    enddo
                        endif

            if(.not.jf(25)) then
                        print *,'User input for F10.7D'
                        read(5,*) oar(41,1)
                        do i=2,100
                                    oar(41,i)=oar(41,1)
                                    enddo
                        endif

            if(.not.jf(32)) then
                        print *,'User input for F10.7_81d'
                        read(5,*) oar(46,1)
                        do i=2,100
                                    oar(46,i)=oar(46,1)
                                    enddo
                        endif

c end of user input
c

        num1=(vend-vbeg)/vstp+1
        numstp=iabs(num1)
        if(numstp.GT.nummax) numstp=nummax

        if(jf(29)) then
             popt='IRI2001'
        else
             if(jf(30)) then
                   popt='IRIcorr'
             else
                   popt='NeQuick'
             endif
        endif
        map='URSI'
        if(jf(5)) map='CCIR'

        if(jf(39)) then
             hopt='CCIR-M3000'
        else
             if(jf(40)) then
                   hopt='AMTB-2013'
             else
                   hopt='Shubin2015'
             endif
        endif

        if(jf(4)) then
             bopt='BIl-2000'
        else
             if(jf(31)) then
                   bopt='ABT-2009'
             else
                   bopt='Gul-1987'
             endif
        endif

        iopt='RBV10+TTS05'
        if(jf(6)) iopt='DS95 + DY85'

        dopt='FT01+DRS95'
        if(jf(24)) dopt='IRI-1990'

        sopt='off'
        if(jf(26)) sopt='on '
        seopt='off'
        if(jf(35)) seopt='on '

        topt='TBT-2012'
        if(jf(23)) topt='BIl-1985'
        tsopt=' with'
        if(jf(23)) tsopt='  w/o'

        if(jf(19)) then
              f1opt='Scotto-97 no L'
              if(.not.jf(20)) f1opt='Scotto-97 with L'
        else
              f1opt='IRI-95'
        endif

        rzopt=' user input'
        if(jf(17)) rzopt=' '
        igopt=' user input'
        if(jf(27)) igopt=' '
        fdopt=' user input'
        if(jf(25)) fdopt=' '
        f8opt=' user input'
        if(jf(32)) f8opt=' '

        hxx=hx
        jmag=jm
        mmdd=imd

c calling IRI subroutine
c
        phour=hour
        call iri_web(jmag,jf,xlat,xlon,iy,mmdd,iut,hour,
     &          hxx,htec_max,ivar,vbeg,vend,vstp,outf,oar)

c preparation of results page
c
        write(7,3991) iy,mmdd,phour,timev(iut+1),
     &        coorv(jmag+1),xlat,xlon,hxx
        if(jf(1)) then
       		write(7,3314) popt
            if(jf(8)) then
            	write(7,301) map
                write(7,3291) sopt
                endif
            if(jf(9)) write(7,303) hopt
            write(7,309) bopt
            write(7,3295) f1opt
            write(7,3299) seopt
            write(7,3081) dopt
            numi=numstp
            if(ivar.eq.1) numi=1
            do j=1,6
                ij=jfi(j)
                if(.not.jf(ij)) then
                    write(7,302) pname(j)
                    write(7,402) (oar(j,i),i=1,numi)
                    endif
                enddo
                if(.not.jf(43)) then
                    write(7,302) pname(7)
                    write(7,402) (oar(10,i),i=1,numi)
                    endif
            endif

        if(jf(2)) write(7,3292) topt,tsopt
        if(jf(3)) write(7,329) iopt

        if(ivar.eq.1) then
                if(oar(3,1).lt.1.) oar(4,1)=0.
                yp2=0
                if(oar(3,1).gt.0.0) yp2=oar(3,1)/1.e6
                write(7,213) oar(1,1)/1.E6,yp2,oar(5,1)/1.E6
                write(7,214) oar(2,1),oar(4,1),oar(6,1)
        else
                write(7,307)
        endif

        write(7,211) oar(23,1),oar(25,1),oar(27,1)

        write(7,223) oar(33,1),rzopt
        write(7,2231) oar(39,1),igopt
        write(7,2238) oar(41,1),fdopt
        write(7,2237) oar(46,1),f8opt

        if(htec_max.gt.50.0) write(7,3914) htec_max

3991    format(///'yyyy/mmdd(or -ddd)/hh.h):',I4,'/',I4,'/',F4.1,
     &      A2,2X,A4,' Lat/Long/Alt=',F5.1,'/',F6.1,'/',F6.1/)
3914    format(/'TEC [1.E16 m-2] is obtained by numerical integ',
     &     'ration in 1km steps'/'  from 50 to ',f6.1,' km. ',
     &     't is the percentage of TEC above the F peak.')
3916    format(/'M3000F2: Propagation factor related to hmF2'/
     &     'B0: bottomside thickness parameter.')
301     format(A4,' maps are used for the F2 peak density (NmF2)')
302     format(A9,' provided by user:')
402     format(7(1PE10.3))
303     format(A10,'model is used for F2 peak height (hmF2)')
304     format(A25,'=',F5.1,') provided by user')
307     format(1x/'Solar and magnetic parameter for the 1st profile',
     &          ' point:')
309     format(A8,' option is used for the bottomside thickness ',
     &          'parameter B0')
3081    format(A10,' option is used for D-region')
329     format(A11,' option is used for ion composition')
3291    format('foF2 STORM model is turned ',A3)
3299    format('foE auroral storm model is turned ',A3)
3292    format(A8,A5,' solar dependence is used for the electron',
     &          ' temperature')
3293    format(A8,' option is used for the D-region Ne')
3295    format(A16,' option is used for the F1 occurrence probability')
211     format('Solar Zenith Angle/degree',28X,F6.1/
     &          'Dip (Magnetic Inclination)/degree',21X,
     &          F6.2/'Modip (Modified Dip)/degree',27X,F6.2)
223     format('Solar Sunspot Number (12-months running mean) Rz12',
     &          4X,F5.1,A11)
2231    format('Ionospheric-Effective Solar Index IG12',16X,
     &          F5.1,A11)
2238    format('Solar radio flux F10.7 (daily)',24X,
     &          F5.1,A11)
2237    format('Solar radio flux F10.7 (81-day average)',15X,
     &          F5.1,A11)
213     format(/'Peak Densities/cm-3: NmF2=',F9.1,'   NmF1=',F9.1,
     &          '   NmE=',F9.1)
214     format('Peak Heights/km:     hmF2=',F9.2,'   hmF1=',F9.2,
     &          '   hmE=',F9.2/)
3314    format(A7,' is used for topside Ne profile')
c
c table head .......................................................
c

        agnr=7          !output unit number
        xtex=imz(ivar)
        if(jmag.gt.0.and.(ivar.eq.2.or.ivar.eq.3)) xtex='GEOM'
        if(iut.gt.0.and.ivar.eq.8) xtex='U.T.'

        IF(PIKTAB.EQ.4) WRITE(7,8199)
        IF(PIKTAB.EQ.3) WRITE(7,8191) ITEXT(IVAR),
     &    (pna(pad1(j)),j=1,6),xtex,(uni(pad1(j)),j=1,6)
        IF(PIKTAB.EQ.2) WRITE(7,8194) ITEXT(IVAR),xtex
        IF(PIKTAB.EQ.1) WRITE(7,8192) ITEXT(IVAR),xtex

        IF(PIKTAB.EQ.0) THEN
        	if(jf(22)) then
        		WRITE(7,8193) ITEXT(IVAR),xtex
        	else
        		WRITE(7,9193) ITEXT(IVAR),xtex
        	endif
			ENDIF
8191  FORMAT(/'-'/2X,A5,6A10/3X,A4,6A10)
8192  FORMAT(/'-'/2X,A5,6X,'PEAK ALTITUDES IN KM',8X,'PEAK DEN',
     &  'SITIES IN cm-3  TEC top/%'/3X,A4,'    hmF2  hmF1   hmE   ',
     &  'hmD      NmF2   NmF1    NmE    NmD  1E16m-2')
8194  FORMAT(/'-'/2X,A5,3X,'M3000   B0',3X,'B1   E-VALLEY',7X,
     &  'PLASMA ','FREQUENCIES / MHz'/3X,A4,11X,'km       W/km ',
     &  ' Depth',5X,'foF2   foF1   foE   foD')
8193  FORMAT(/'-'/1X,A5,' ELECTRON DENSITY   TEMPERATURES ',
     &  5X,'ION PERCENTAGES[%]*10',4X,'1E16m-2'/2X,A4,' Ne/cm-3 Ne/',
     &  'NmF2 Tn/K  Ti/K  Te/K  O+  N+  H+ He+ O2+ NO+ Clust TEC t/%')
9193  FORMAT(/'-'/1X,A5,' ELECTRON DENSITY   TEMPERATURES ',
     &  4X,'ION DENSITIES[cm-3]/100',4x,'1E16m-2'/2X,A4,' Ne/cm-3 Ne/',
     &  'NmF2 Tn/K  Ti/K  Te/K  O+  N+  H+ He+ O2+ NO+ Clust TEC t/%')
8199  FORMAT(/'-'/1X,'h',8X,' D-REGION ELECTRON DENSITY IN CM-3'/
     &  1X,'km',18X,'DRS-95: Stratos Warming/Winter Anomaly'/5X,
     &  'IRI-07',4x,'FIRI  SW/WA=0/0  0.5/0   1/0    0/0.5    0/1')

c
c output: D-region PIKTAB=4
c
c D-REGION ELECTRON DENSITY IN CM-3:
c    IRI-07    FIRI  Danilov:SW/WA=0/0  0.5/0   1/0    0/0.5    0/1
c                    DRS-95: Stratos Warming/Winter Anomaly
c

		if(piktab.eq.4) then
            do 2591 lix=1,77
            	jdprof(lix)=-1
            	dichte=outf(14,lix)
2591            if(dichte.gt.0.) jdprof(lix)=int(dichte/1.e6+0.5)
			do 2592 lix=1,11
				ihtemp=55+lix*5
            	WRITE(7,3810) ihtemp,jdprof(lix),jdprof(lix+11),
     &  			jdprof(lix+22),jdprof(lix+33),jdprof(lix+44),
     &  			jdprof(lix+55),jdprof(lix+66)
2592		    continue
3810    FORMAT(I3,7I8)
			goto 2357
		 	endif

        xcor=vbeg

        do 1234 li=1,numstp

c
c output: peak densities and altitudes PIKTAB=1
c

      IF(PIKTAB.eq.1) THEN
        if(oar(3,li).lt.1.) oar(4,li)=0.
        iyp1=int(oar(1,li)/1.e6+.5)
        iyp2=0
        if(oar(3,li).gt.0.0) iyp2=int(oar(3,li)/1.e6+.5)
        iyp3=int(oar(5,li)/1.e6+.5)
        iyp4=int(oar(7,li)/1.e6+.5)
            tec=oar(37,li)
        if(tec.gt.0.0) then
            tec=tec/1.e16
            itopp=int(oar(38,li)+.5)
        else
            tec=-1.0
            itopp=-1
        endif
        WRITE(7,3910) XCOR,oar(2,li),oar(4,li),oar(6,li),oar(8,li),
     &    iyp1,iyp2,iyp3,iyp4,tec,itopp
3910    FORMAT(F7.1,2X,4F6.1,1X,I9,3I7,1X,F6.2,I4)
        GOTO 1234
      ENDIF
c
c output: plasma frequencies and profile parameters  PIKTAB=2
c

      IF(PIKTAB.eq.2) THEN
        if(oar(3,li).lt.1.) oar(4,li)=0.
        fyp1=SQRT(oar(1,li)/1.24E10)
        fyp2=0
        if(oar(3,li).gt.0.0) fyp2=SQRT(oar(3,li)/1.24E10)
        fyp3=SQRT(oar(5,li)/1.24E10)
        fyp4=SQRT(oar(7,li)/1.24E10)
            tec=oar(37,li)
        if(tec.gt.0.0) then
            tec=tec/1.e16
            itopp=int(oar(38,li)+.5)
        else
            tec=-1.0
            itopp=-1
        endif
        wvalley=oar(12,li)-oar(6,li)
        dvalley=0.0
        if(oar(5,li).gt.0.0) dvalley=oar(11,li)/oar(5,li)
        WRITE(7,3950) XCOR,oar(36,li),oar(10,li),oar(35,li),
     &    wvalley,dvalley,fyp1,fyp2,fyp3,fyp4
3950    FORMAT(F7.1,2X,F6.4,F6.1,F4.1,F6.1,F8.4,1X,4F7.3)
        GOTO 1234
      ENDIF
c
c output: 6 parameters of your choice    PIKTAB=3
c

      IF(PIKTAB.eq.3) THEN
c        if(pad1.eq.45.and.oar(pad1,li).le.0.0) oar(pad1,li)=-1.
c        if(pad2.eq.45.and.oar(pad2,li).le.0.0) oar(pad2,li)=-1.
c        if(pad3.eq.45.and.oar(pad3,li).le.0.0) oar(pad3,li)=-1.
        WRITE(7,3919) XCOR,oar(pad1(1),li),oar(pad1(2),li),
     &        oar(pad1(3),li),oar(pad1(4),li),oar(pad1(5),li),
     &        oar(pad1(6),li)
3919    FORMAT(F7.1,6(1X,1PE9.2))
        GOTO 1234
      ENDIF
c
c output: special for test purposes    PIKTAB=5
c

      IF(PIKTAB.eq.5) THEN
c ----------- B0, B1 ----------------
c        WRITE(8,4919) XCOR,oar(10,li),oar(35,li)
c4919    FORMAT(F7.1,2X,F6.2,2X,F5.3)
c ----------- Te ----------------
        WRITE(8,4919) XCOR,outf(2,li),outf(3,li),outf(4,li)
4919    FORMAT(F7.1,2X,F6.1,2X,F6.1,2X,F6.1)
c ----------- Ne, TEC ----------------
c        WRITE(8,4919) XCOR,outf(1,li),oar(37,li)
c4919    FORMAT(F7.1,2X,E12.5,2X,E12.5)
c ----------- Ni ----------------
c        type*,XCOR,outf(1,li),outf(5,li),outf(6,li),
c     &   outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li)
c        WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li),
c     &   outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li)
c4919    FORMAT(F7.1,2X,E12.5,7F10.4)
c ----------- hmF2 ----------------
c        type*,XCOR,oar(26,li),oar(2,li),oar(36,li),
c     &   sqrt(oar(1,li)/oar(5,li)),oar(33,li),oar(39,li)
c        WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li),
c     &   outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li)
c4919    FORMAT(F7.1,2X,E12.5,7F10.4)
c ----------- auroral boundary ----------------
c        print *,XCOR,oar(55,li),oar(56,li),oar(54,li),
c     &     oar(57,li),oar(54,li)-oar(57,li),oar(58,li)
cc        WRITE(8,4919) XCOR,outf(1,li),outf(5,li),outf(6,li),
cc     &   outf(7,li),outf(8,li),outf(9,li),outf(10,li),outf(11,li)
cc4919    FORMAT(F7.1,2X,E12.5,7F10.4)

        GOTO 1234
      ENDIF
c
c output: standard
c
        if(ivar.eq.1) then
                oar(1,li)=oar(1,1)
                oar(37,li)=oar(37,1)
                oar(38,li)=oar(38,1)
                endif
        jne=int(outf(1,li)/1.e6+.5)
        xner=outf(1,li)/oar(1,li)
        jtn=int(outf(2,li)+.5)
        jti=int(outf(3,li)+.5)
        jte=int(outf(4,li)+.5)
        scid=1.0E-8
        if(jf(22)) scid=10.
        jio=INT(OUTF(5,li)*scid+.5)
        jih=INT(OUTF(6,li)*scid+.5)
        jihe=INT(OUTF(7,li)*scid+.5)
        jio2=INT(OUTF(8,li)*scid+.5)
        jino=INT(OUTF(9,li)*scid+.5)
        jicl=INT(OUTF(10,li)*scid+.5)
        jin=INT(OUTF(11,li)*scid+.5)
c        print *,'O+,N+,O2+,NO+=',outf(5,li),
c     &    outf(11,li),outf(9,li),outf(8,li)
        if(outf(1,li).lt.0) jne=-1
        if(outf(1,li).lt.0) xner=-1.
        if(outf(2,li).lt.0) jtn=-1
        if(outf(3,li).lt.0) jti=-1
        if(outf(4,li).lt.0) jte=-1
        if(outf(5,li).lt.0) jio=-1
        if(outf(6,li).lt.0) jih=-1
        if(outf(7,li).lt.0) jihe=-1
        if(outf(8,li).lt.0) jio2=-1
        if(outf(9,li).lt.0) jino=-1
        if(outf(10,li).lt.0) jicl=-1
        if(outf(11,li).lt.0) jin=-1
            tec=oar(37,li)
        if(tec.gt.0.0) then
            tec=tec/1.e16
            itopp=int(oar(38,li)+.5)
        else
            tec=-1.0
            itopp=-1
        endif
c        print *, XCOR,jne,xner,jtn,jti,jte,jio,jin,
c     &        jih,jihe,jino,jio2,jicl,tec,itopp
        WRITE(7,7117) XCOR,jne,xner,jtn,jti,jte,jio,jin,
     &        jih,jihe,jio2,jino,jicl,tec,itopp
7117    FORMAT(F6.1,I8,1x,F6.3,3I6,7I4,f6.1,i4)

1234    xcor=xcor+vstp

2357    print *,'Enter 0 to exit or 1 to generate another profile?'
        read(5,*) icontinue
        if (icontinue.gt.0) goto 1

            stop
            end

ext: for


text:
      SUBROUTINE MB03YA( WANTT, WANTQ, WANTZ, N, ILO, IHI, ILOQ, IHIQ,
     $                   POS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To annihilate one or two entries on the subdiagonal of the
C     Hessenberg matrix A for dealing with zero elements on the diagonal
C     of the triangular matrix B.
C
C     MB03YA is an auxiliary routine called by SLICOT Library routines
C     MB03XP and MB03YD.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     WANTT   LOGICAL
C             Indicates whether the user wishes to compute the full
C             Schur form or the eigenvalues only, as follows:
C             = .TRUE. :  Compute the full Schur form;
C             = .FALSE.:  compute the eigenvalues only.
C
C     WANTQ   LOGICAL
C             Indicates whether or not the user wishes to accumulate
C             the matrix Q as follows:
C             = .TRUE. :  The matrix Q is updated;
C             = .FALSE.:  the matrix Q is not required.
C
C     WANTZ   LOGICAL
C             Indicates whether or not the user wishes to accumulate
C             the matrix Z as follows:
C             = .TRUE. :  The matrix Z is updated;
C             = .FALSE.:  the matrix Z is not required.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices A and B. N >= 0.
C
C     ILO     (input) INTEGER
C     IHI     (input) INTEGER
C             It is assumed that the matrices A and B are already
C             (quasi) upper triangular in rows and columns 1:ILO-1 and
C             IHI+1:N. The routine works primarily with the submatrices
C             in rows and columns ILO to IHI, but applies the
C             transformations to all the rows and columns of the
C             matrices A and B, if WANTT = .TRUE..
C             1 <= ILO <= max(1,N); min(ILO,N) <= IHI <= N.
C
C     ILOQ    (input) INTEGER
C     IHIQ    (input) INTEGER
C             Specify the rows of Q and Z to which transformations
C             must be applied if WANTQ = .TRUE. and WANTZ = .TRUE.,
C             respectively.
C             1 <= ILOQ <= ILO; IHI <= IHIQ <= N.
C
C     POS     (input) INTEGER
C             The position of the zero element on the diagonal of B.
C             ILO <= POS <= IHI.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the upper Hessenberg matrix A.
C             On exit, the leading N-by-N part of this array contains
C             the updated matrix A where A(POS,POS-1) = 0, if POS > ILO,
C             and A(POS+1,POS) = 0, if POS < IHI.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,N)
C             On entry, the leading N-by-N part of this array must
C             contain an upper triangular matrix B with B(POS,POS) = 0.
C             On exit, the leading N-by-N part of this array contains
C             the updated upper triangular matrix B.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1,N).
C
C     Q       (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
C             On entry, if WANTQ = .TRUE., then the leading N-by-N part
C             of this array must contain the current matrix Q of
C             transformations accumulated by MB03XP.
C             On exit, if WANTQ = .TRUE., then the leading N-by-N part
C             of this array contains the matrix Q updated in the
C             submatrix Q(ILOQ:IHIQ,ILO:IHI).
C             If WANTQ = .FALSE., Q is not referenced.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.  LDQ >= 1.
C             If WANTQ = .TRUE., LDQ >= MAX(1,N).
C
C     Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C             On entry, if WANTZ = .TRUE., then the leading N-by-N part
C             of this array must contain the current matrix Z of
C             transformations accumulated by MB03XP.
C             On exit, if WANTZ = .TRUE., then the leading N-by-N part
C             of this array contains the matrix Z updated in the
C             submatrix Z(ILOQ:IHIQ,ILO:IHI).
C             If WANTZ = .FALSE., Z is not referenced.
C
C     LDZ     INTEGER
C             The leading dimension of the array Z.  LDZ >= 1.
C             If WANTZ = .TRUE., LDZ >= MAX(1,N).
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The method is illustrated by Wilkinson diagrams for N = 5,
C     POS = 3:
C
C           [ x x x x x ]       [ x x x x x ]
C           [ x x x x x ]       [ o x x x x ]
C       A = [ o x x x x ],  B = [ o o o x x ].
C           [ o o x x x ]       [ o o o x x ]
C           [ o o o x x ]       [ o o o o x ]
C
C     First, a QR factorization is applied to A(1:3,1:3) and the
C     resulting nonzero in the updated matrix B is immediately
C     annihilated by a Givens rotation acting on columns 1 and 2:
C
C           [ x x x x x ]       [ x x x x x ]
C           [ x x x x x ]       [ o x x x x ]
C       A = [ o o x x x ],  B = [ o o o x x ].
C           [ o o x x x ]       [ o o o x x ]
C           [ o o o x x ]       [ o o o o x ]
C
C     Secondly, an RQ factorization is applied to A(4:5,4:5) and the
C     resulting nonzero in the updated matrix B is immediately
C     annihilated by a Givens rotation acting on rows 4 and 5:
C
C           [ x x x x x ]       [ x x x x x ]
C           [ x x x x x ]       [ o x x x x ]
C       A = [ o o x x x ],  B = [ o o o x x ].
C           [ o o o x x ]       [ o o o x x ]
C           [ o o o x x ]       [ o o o o x ]
C
C     REFERENCES
C
C     [1] Bojanczyk, A.W., Golub, G.H., and Van Dooren, P.
C         The periodic Schur decomposition: Algorithms and applications.
C         Proc. of the SPIE Conference (F.T. Luk, Ed.), 1770, pp. 31-42,
C         1992.
C
C     NUMERICAL ASPECTS
C
C     The algorithm requires O(N**2) floating point operations and is
C     backward stable.
C
C     CONTRIBUTORS
C
C     D. Kressner, Technical Univ. Berlin, Germany, and
C     P. Benner, Technical Univ. Chemnitz, Germany, December 2003.
C
C     REVISIONS
C
C     V. Sima, June 2008 (SLICOT version of the HAPACK routine DLADFB).
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO
      PARAMETER          ( ZERO = 0.0D+0 )
C     .. Scalar Arguments ..
      LOGICAL            WANTQ, WANTT, WANTZ
      INTEGER            IHI, IHIQ, ILO, ILOQ, INFO, LDA, LDB, LDQ, LDZ,
     $                   N, POS
C     .. Array Arguments ..
      DOUBLE PRECISION   A(LDA,*), B(LDB,*), Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      INTEGER            I1, I2, J, NQ
      DOUBLE PRECISION   CS, SN, TEMP
C     .. External Subroutines ..
      EXTERNAL           DLARTG, DROT, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN
C
C     .. Executable Statements ..
C
C     Check the scalar input parameters.
C
      INFO = 0
      NQ = IHIQ - ILOQ + 1
      IF ( N.LT.0 ) THEN
         INFO = -4
      ELSE IF ( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF ( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
         INFO = -6
      ELSE IF ( ILOQ.LT.1 .OR. ILOQ.GT.ILO ) THEN
         INFO = -7
      ELSE IF ( IHIQ.LT.IHI .OR. IHIQ.GT.N ) THEN
         INFO = -8
      ELSE IF ( POS.LT.ILO .OR. POS.GT.IHI ) THEN
         INFO = -9
      ELSE IF ( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -11
      ELSE IF ( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -13
      ELSE IF ( LDQ.LT.1 .OR. WANTQ .AND. LDQ.LT.N ) THEN
         INFO = -15
      ELSE IF ( LDZ.LT.1 .OR. WANTZ .AND. LDZ.LT.N ) THEN
         INFO = -17
      END IF
C
C     Return if there were illegal values.
C
      IF ( INFO.NE.0 ) THEN
         CALL XERBLA( 'MB03YA', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( N.EQ.0 )
     $   RETURN
C
      IF ( WANTT ) THEN
         I1 = 1
         I2 = N
      ELSE
         I1 = ILO
         I2 = IHI
      END IF
C
C     Apply a zero-shifted QR step.
C
      DO 10  J = ILO, POS-1
         TEMP = A(J,J)
         CALL DLARTG( TEMP, A(J+1,J), CS, SN, A(J,J) )
         A(J+1,J) = ZERO
         CALL DROT( I2-J, A(J,J+1), LDA, A(J+1,J+1), LDA, CS, SN )
         CALL DROT( MIN(J,POS-2)-I1+2, B(I1,J), 1, B(I1,J+1), 1, CS,
     $              SN )
         IF ( WANTQ )
     $      CALL DROT( NQ, Q(ILOQ,J), 1, Q(ILOQ,J+1), 1, CS, SN )
   10 CONTINUE
      DO 20  J = ILO, POS-2
         TEMP = B(J,J)
         CALL DLARTG( TEMP, B(J+1,J), CS, SN, B(J,J) )
         B(J+1,J) = ZERO
         CALL DROT( I2-J, B(J,J+1), LDB, B(J+1,J+1), LDB, CS, SN )
         CALL DROT( J-I1+2, A(I1,J), 1, A(I1,J+1), 1, CS, SN )
         IF ( WANTZ )
     $      CALL DROT( NQ, Z(ILOQ,J), 1, Z(ILOQ,J+1), 1, CS, SN )
   20 CONTINUE
C
C     Apply a zero-shifted RQ step.
C
      DO 30  J = IHI, POS+1, -1
         TEMP = A(J,J)
         CALL DLARTG( TEMP, A(J,J-1), CS, SN, A(J,J) )
         A(J,J-1) = ZERO
         SN = -SN
         CALL DROT( J-I1, A(I1,J-1), 1, A(I1,J), 1, CS, SN )
         CALL DROT( I2 - MAX( J-1,POS+1 ) + 1, B(J-1,MAX( J-1,POS+1 )),
     $              LDB, B(J,MAX(J-1,POS+1)), LDB, CS, SN )
         IF ( WANTZ )
     $      CALL DROT( NQ, Z(ILOQ,J-1), 1, Z(ILOQ,J), 1, CS, SN )
   30 CONTINUE
      DO 40  J = IHI, POS+2, -1
         TEMP = B(J,J)
         CALL DLARTG( TEMP, B(J,J-1), CS, SN, B(J,J) )
         B(J,J-1) = ZERO
         SN = -SN
         CALL DROT( J-I1, B(I1,J-1), 1, B(I1,J), 1, CS, SN )
         CALL DROT( I2-J+2, A(J-1,J-1), LDA, A(J,J-1), LDA, CS, SN )
         IF ( WANTQ )
     $      CALL DROT( NQ, Q(ILOQ,J-1), 1, Q(ILOQ,J), 1, CS, SN )
   40 CONTINUE
      RETURN
C *** Last line of MB03YA ***
      END

ext: f

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


text:
*dk,reconloop3d
      subroutine reconloop3d(iobtuse,toldamage)
C
C #####################################################################
C
C     PURPOSE -
C
C        None
C
C     INPUT ARGUMENTS -
C
C        None
C
C     OUTPUT ARGUMENTS -
C
C        None
C
C     CHANGE HISTORY -
C
C        $Log: reconloop3d.f,v $
C        Revision 2.00  2007/11/09 20:04:00  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.1   15 May 2001 14:45:54   kuprat
CPVCS    We now pass TOLDAMAGE to RECON2.
CPVCS
CPVCS       Rev 1.0   Fri Apr 17 15:34:38 1998   dcg
CPVCS    Initial revision.
C
C ######################################################################
C
      implicit none
C
      include 'local_element.h'
      include 'neibor.h'
C
      integer iobtuse,ierror,nsdtopo,length,icmotype,i,npoints,ilen
     &   ,itype,ntets,ierrmat
      real*8 toldamage

      character*132 logmess
      character*32 cmo, isubname
C
      isubname='reconloop3d'
C
C  get info from mesh object and make sure it is a 2d mesh
      call cmo_get_name(cmo,ierror)
      if(ierror.ne.0) then
         write(logmess,'(a)')'RECON3D: no 3d mesh object'
         call writloga('default',0,logmess,0,ierror)
         return
      endif
      call cmo_get_info('ndimensions_topo',cmo,
     *                  nsdtopo,length,icmotype,ierror)
      if(nsdtopo.ne.3) then
         write(logmess,'(a)')'RECON3D: bad 3d mesh object'
         call writloga('default',0,logmess,0,ierror)
         return
      endif
C
C
      do i=1,1
         write(logmess,'(a,i5)') 'RECON3D reconnection loop: ',i
         call writloga('default',0,logmess,0,ierror)
         call cmo_get_info('nnodes',cmo,npoints,ilen,itype,ierror)
         call cmo_get_info('nelements',cmo,ntets,ilen,itype,ierror)
         irecnall=1
         nrecon=ntets
         call recon2(npoints,ntets,toldamage)
         call cmo_get_info('nnodes',cmo,npoints,ilen,itype,ierror)
         call cmo_get_info('nelements',cmo,ntets,ilen,itype,ierror)
         call cmo_newlen(cmo,ierror)
C*****   call matbld3d(cmo,ierrmat)
         if(iobtuse.eq.1) then
            call refine_edge_3d(cmo,ierror)
C*****      if(ierror.gt.0) goto 100
            call cmo_get_info('nnodes',cmo,npoints,ilen,itype,ierror)
            call cmo_get_info('nelements',cmo,ntets,ilen,itype,ierror)
            irecnall=1
            nrecon=ntets
            call recon2(npoints,ntets,toldamage)
            call cmo_get_info('nnodes',cmo,npoints,ilen,itype,ierror)
            call cmo_get_info('nelements',cmo,ntets,ilen,itype,ierror)
            call cmo_newlen(cmo,ierror)
            ierrmat=-1
C*****      call matbld3d(cmo,ierrmat)
            if(ierrmat.lt.0) goto 110
         else
            goto 100
         endif
      enddo
 100  continue
 110  continue
C
      goto 9999
 9999 continue
C
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


      subroutine hru_lte_control (isd)

      use hru_lte_module
      use hydrograph_module
      use output_landscape_module
      use basin_module
      use climate_module
      use time_module
      use plant_data_module
      use conditional_module

      implicit none

      integer :: isd                    !             |
      !real :: timeint(1000)
      real :: a1 = .2                   !             |
      real :: a2 = .8                   !             |
      integer :: ihlt_db                !             |
      integer :: iwgn                   !             |
      integer :: iplt                   !none         |counter
      real :: precip                    !mm           |precipitation
      real :: tmax                      !deg C        |maximum average monthly temperature
      real :: tmin                      !deg C        |minimum average monthly temperature
      real :: raobs                     !             |
      real :: rmx                       !             |
      real :: tave                      !             |
      real :: yield                     !             |
      real :: ws                        !             |
      real :: strsair                   !             |
      real :: snowfall                  !             |
      real :: snowmelt                  !             |
      real :: runoff                    !             |
      real :: xx                        !             |
      real :: exp                       !             |
      real :: r2                        !             |
      real :: amax1                     !             |
      real :: cn_sd                     !             |
      real :: precipeff                 !             |
      real :: xxi                       !             |
      real :: xsd                       !             |
      real :: ch                        !             |
      real :: tan                       !             |
      real :: h                         !             |
      real :: acos                      !             |
      real :: ramm                      !MJ/m2        |extraterrestrial radiation
      real :: pet                       !             |
      real :: tstress                   !             |sum of temperature stress
      real :: tk                        !deg C        |mean air temperature
      real :: alb                       !none         |albedo when soil is moist
      real :: d                         !cm           |displacement height for plant type
      real :: gma                       !kPa/deg C    |psychrometric constant
      real :: ho                        !none         |variable to hold intermediate calculation
                                        !             |result
      real :: aph                       !             |
      real :: aet                       !mm            |sum of actual et during growing season (for hi water stress)
      real :: b1                        !             |
      real :: delg                      !             |
      real :: parad                     !             |
      real :: drymat                    !             |
      real :: satco                     !             |
      real :: pl_aerfac                 !             |
      integer :: iend                   !             |points to rule set in d_table
      integer :: istart                 !none         |points to rule set in d_table
      real :: scparm                    !             |
      real :: air                       !             |
      real :: amin1                     !             |
      real :: tgx                       !             |
      real :: rto                       !none         |cloud cover factor
      real :: reg                       !             |
      real :: deltalai                  !             |
      real :: sw_excess                 !mm H2O       |amount of water in excess of field capacity
                                        !             |stored in soil layer on the current day
      real :: swf                       !cfu          |fraction of manure containing active colony forming units
      real :: flowlat                   !             |
      real :: f                         !             |
      real :: ff                        !             |
      real :: flow_tile                 !             |
      real :: perc                      !             |
      real :: revap                     !mm           |revap
      real :: percdeep                  !             |
      real :: chflow                    !             |
      real :: chflow_m3                 !m^3/s        |Runoff in CMS
      real :: runoff_m3                 !             |
      real :: bf_m3                     !             |
      real :: peakr                     !m^3/s        |peak runoff rate in channel
      real :: peakrbf                   !             |
      real :: sedin                     !             |
      real :: qssubconc                 !             |
      real :: qssub                     !             |
      real :: cnv                       !none         |conversion factor (mm => m^3)
      real :: wndspd                    !none         |windspeed
      real :: rhum                      !none         |relative humidity

      ihlt_db = ob(icmd)%props
      iwst = ob(icmd)%wst
      iwgn = wst(iwst)%wco%wgn
      iplt = hlt(isd)%iplant
      precip = wst(iwst)%weat%precip
      tmax = wst(iwst)%weat%tmax
      tmin = wst(iwst)%weat%tmin
      raobs = wst(iwst)%weat%solrad
      rmx = wst(iwst)%weat%solradmx

      tave  = (tmax + tmin) / 2.
      yield = 0.
      ws = 0.
      strsair = 1.
      tstress = 0.
      snowfall = 0.
      snowmelt = 0.
      air = 0.

          !! compute curve number
          xx = hlt(isd)%wrt1 - hlt(isd)%wrt2 * hlt(isd)%sw
          if (xx < -20.) xx = -20.
          if (xx > 20.) xx = 20.
          if ((hlt(isd)%sw + Exp(xx)) > 0.001) then
            r2 = hlt(isd)%smx * (1. - hlt(isd)%sw / (hlt(isd)%sw + Exp(xx)))
          end if
          r2 = amax1(3.,r2)
          cn_sd = 25400. / (r2 + 254.)

          IF (tave .lt.0.) THEN
            ! IF ave temp < 0  compute snowfall
            snowfall = precip
            hlt(isd)%snow = hlt(isd)%snow + precip
            runoff = 0.
          ELSE
            snowfall = 0.
            ! IF ave temp > 0  compute runoff
            snowmelt = 4.57 * tave
            IF (snowmelt > hlt(isd)%snow) THEN
              snowmelt = hlt(isd)%snow
              hlt(isd)%snow = 0.
            ELSE
              hlt(isd)%snow = hlt(isd)%snow - snowmelt
            END IF

            precipeff = precip + snowmelt
            xx = precipeff - a1 * r2
            IF (xx.gt.0.) THEN
              runoff = xx ** 2 / (precipeff + a2 * r2)
            ELSE
              runoff = 0.
            END IF
            hlt(isd)%sw = hlt(isd)%sw + (precipeff - runoff)
          END IF

          xxi = 30. * time%mo - 15.
          xsd = .4102 * SIN((xxi-80.25)/58.13)
          ch = -hlt(isd)%yls * tan(xsd) / hlt(isd)%ylc
          IF (ch.lt.1.) THEN
            IF (ch.le.-1.) THEN
              h = 3.1415
            ELSE
              h = acos(ch)
            END IF
          ELSE
            h = 0.
          END IF

          IF (hlt_db(ihlt_db)%ipet .eq. "harg") THEN
            ! compute potential et with Hargreaves Method (harg)
            ramm = rmx / (2.5 - .0022 * tave)
            pet = .0032 * ramm * (tave +17.8) * (tmax - tmin) ** .6
          ELSE
            ! compute potential et with Preistley-Taylor Method (p_t)
            tk = tave  + 273.
            alb = .23
            d = EXP(21.255 - 5304. / tk) * 5304. / tk ** 2
            gma = d / (d +.68)
            ho = 23.9 * raobs * (1. - alb) / 58.3
            aph = 1.28
            pet = aph * ho * gma
          END IF
          pet = hlt(isd)%etco * pet
!
!         compute actual et
!
          xx = 1. - hlt(isd)%sw / hlt(isd)%awc
          IF (xx.lt.0.0001) xx = 0.0001
          aet = pet * EXP(-xx)                  !at sw=0, exp(-xx)=.368
          aet = pet * hlt(isd)%sw / hlt(isd)%awc  !try linear- may need s-curve

!         begin growth for plants
          if (hlt(isd)%gro == "n") then
            ! istart points to rule set in d_table
            istart = hlt(isd)%start
            d_tbl => dtbl_lum(istart)
            call conditions (isd, istart)
            call actions (isd, icmd, istart)
          end if

!         end growth for plants
          if (hlt(isd)%gro == "y") then
            ! sum aet and pet for water stress on hi
            hlt(isd)%aet = hlt(isd)%aet + aet
            hlt(isd)%pet = hlt(isd)%pet + pet
            ! iend points to rule set in d_table
            iend = hlt(isd)%end
            d_tbl => dtbl_lum(iend)
            call conditions (isd, icmd)
            call actions (isd, icmd, iend)
          end if

         ! calc yield, print max lai, dm and yield
          if (pco%mgtout == "y") then
            write (4700,*) isd, time%day, time%yrc, pldb(iplt)%plantnm, hlt(isd)%alai, hlt(isd)%dm, yield
            if (pco%csvout == "y") then
              write (4701,*) isd, time%day, time%yrc, pldb(iplt)%plantnm, hlt(isd)%alai, hlt(isd)%dm, yield
            end if
          end if

!
!         compute plant growth - b1=et adjustment factor b1=1 during growing season b1=.6 IF no
!
          b1 = hlt(isd)%etco - .4        !evap coef ranges from .4-.8 when etco ranges from .8-1.2
          IF (hlt(isd)%gro == "y") THEN
            b1 = hlt(isd)%etco
            delg = (tave - pldb(iplt)%t_base) / hlt(isd)%phu
            IF (delg.lt.0.) THEN
              delg = 0.
            END IF
            hlt(isd)%g = hlt(isd)%g + delg
            parad = .5 * raobs * (1.-EXP(-.65 * (hlt(isd)%alai + .05)))
            drymat = parad * pldb(iplt)%bio_e * hlt(isd)%stress
            ws = aet / pet

            !compute aeration stress
            if (hlt(isd)%sw .gt. hlt(isd)%awc) THEN
              satco = (hlt(isd)%sw - hlt(isd)%awc) / (hlt(isd)%por - hlt(isd)%awc)
              pl_aerfac = .85
              scparm = 100. * (satco - pl_aerfac) / (1.0001 - pl_aerfac)
              if (scparm > 0.) then
                strsair = 1. - (scparm / (scparm + Exp(2.9014 - .03867 * scparm)))
              else
                strsair = 1.
              end if
            end if

            !irrigate IF water stress is < 0.8
            IF (hlt_db(ihlt_db)%irr > "no_irr") THEN
              IF (ws < 0.8) THEN
                air = hlt(isd)%awc - hlt(isd)%sw
                air = amin1 (76.2, air)     ! max application = 3 inches

                ! take from shallow aquifer
                IF (hlt_db(ihlt_db)%irrsrc == "shal_aqu") THEN
                  hlt(isd)%gw = hlt(isd)%gw - air
                  IF (hlt(isd)%gw < 0.) THEN
                    air = air + hlt(isd)%gw
                    hlt(isd)%gw = 0.
                  END IF
                END IF

                ! take from deep aquifer
                IF (hlt_db(ihlt_db)%irrsrc == "deep_aqu") THEN
                  hlt(isd)%gwdeep = hlt(isd)%gwdeep - air
                  IF (hlt(isd)%gwdeep < 0.) THEN
                    air = air + hlt(isd)%gwdeep
                    hlt(isd)%gwdeep = 0.
                  END IF
                END IF
                hlt(isd)%sw = hlt(isd)%sw + air
              END IF
            END IF

            if (tave .gt.pldb(iplt)%t_base) THEN
              tgx = 2. * pldb(iplt)%t_opt - pldb(iplt)%t_base - tave
              rto = ((pldb(iplt)%t_opt - tave ) /(tgx + 1.e-6))**2
              IF (rto.le.200.) THEN
                tstress = EXP(-0.1054*rto)
              ELSE
                tstress = 0.
              END IF
            ELSE
              tstress = 0.
            END IF

!
!         compute boimass and leaf area
!
            reg = amin1(ws,tstress,strsair)
            hlt(isd)%dm = hlt(isd)%dm + reg * drymat
            f = hlt(isd)%g / (hlt(isd)%g + EXP(plcp(iplt)%leaf1 - plcp(iplt)%leaf2 * hlt(isd)%g))
            ff = f - hlt(isd)%hufh
            hlt(isd)%hufh = f
            deltalai = ff * pldb(iplt)%blai * (1.0 - EXP(5.0 *(hlt(isd)%alai - pldb(iplt)%blai))) * sqrt(reg)
            hlt(isd)%alai = hlt(isd)%alai + deltalai
          END IF

!         adjust actual et for growing season
          aet = b1 * aet

          !compute lateral soil flow
          sw_excess = hlt(isd)%sw - hlt(isd)%awc
          if (sw_excess > 0.) then
            swf = (hlt(isd)%sw - hlt(isd)%awc) / (hlt(isd)%por - hlt(isd)%awc)
            flowlat = .024 * swf * hlt(isd)%sc * hlt_db(ihlt_db)%slope / hlt_db(ihlt_db)%slopelen
            flowlat = amin1(hlt(isd)%sw, flowlat)
            hlt(isd)%sw = hlt(isd)%sw - flowlat
          else
            flowlat = 0.
          end if

          !compute tile flow
          sw_excess = hlt(isd)%sw - hlt(isd)%awc
          if (sw_excess > 0. .and. hlt(isd)%tdrain > 0.) then
            flow_tile = sw_excess * (1. - Exp(-24. / hlt(isd)%tdrain))
            flow_tile = amin1(flow_tile, 10.)     !assume a drainage coefficient of 12.5 mm
          else
            flow_tile = 0.
          end if
          flow_tile = amin1(hlt(isd)%sw, flow_tile)
          hlt(isd)%sw = hlt(isd)%sw - flow_tile

          !compute percolation from bottom of soil profile
          sw_excess = hlt(isd)%sw - hlt(isd)%awc * hlt(isd)%perco
          if (sw_excess > 0.) then
            perc = sw_excess * (1. - Exp(-24. / hlt(isd)%hk))
          else
            perc = 0.
          end if
          if (perc.lt.0.) perc = 0.
          perc = amin1(hlt(isd)%sw, perc)
          hlt(isd)%sw = hlt(isd)%sw - perc
          !limit perc with perco_co

          aet = amin1(hlt(isd)%sw, aet)
          hlt(isd)%sw = hlt(isd)%sw - aet

          hlt(isd)%gw = hlt(isd)%gw + perc
          revap = aet * hlt_db(ihlt_db)%revapc
          percdeep = perc * hlt_db(ihlt_db)%percc
          hlt(isd)%gwflow = hlt(isd)%gwflow * hlt_db(ihlt_db)%abf + perc * (1. - hlt_db(ihlt_db)%abf)
          hlt(isd)%gwflow = amin1(hlt(isd)%gwflow, hlt(isd)%gw)
          hlt(isd)%gw = hlt(isd)%gw - hlt(isd)%gwflow

          revap = amin1(revap, hlt(isd)%gw)
          hlt(isd)%gw = hlt(isd)%gw - revap

          percdeep = amin1(percdeep, hlt(isd)%gw)
          hlt(isd)%gw = hlt(isd)%gw - percdeep

          hlt(isd)%gwdeep = hlt(isd)%gwdeep + percdeep

          chflow = runoff + flowlat + flow_tile + hlt(isd)%gwflow

!!        compute channel peak rate using SCS triangular unit hydrograph
          chflow_m3 = 1000. * chflow * hlt_db(ihlt_db)%dakm2
	      runoff_m3 = 1000. * runoff * hlt_db(ihlt_db)%dakm2
	      bf_m3 = 1000. * (flowlat + hlt(isd)%gwflow)*hlt_db(ihlt_db)%dakm2
          peakr = 2. * runoff_m3 / (1.5 * hlt_db(ihlt_db)%tc)
	      peakrbf = bf_m3 / 86400.
          peakr = (peakr + peakrbf)     !* prf

!!        compute sediment yield with MUSLE
          sedin = (runoff * peakr * 1000. * hlt_db(ihlt_db)%dakm2) ** .56 * hlt(isd)%uslefac

	    !! add subsurf sediment - t=ppm*mm*km2/1000.
	    qssubconc = 500.
	    qssub = qssubconc * (flowlat + hlt(isd)%gwflow) * hlt_db(ihlt_db)%dakm2 / 1000.
	    sedin = sedin + qssub

          cnv = hlt_db(ihlt_db)%dakm2 * 1000.

         !   output_waterbal - SWAT-DEG0140
        hltwb_d(isd)%precip = precip
        hltwb_d(isd)%snofall = snowfall
        hltwb_d(isd)%snomlt = snowmelt
        hltwb_d(isd)%surq_gen = runoff
        hltwb_d(isd)%latq = flowlat + hlt(isd)%gwflow
        hltwb_d(isd)%wateryld = chflow
        hltwb_d(isd)%perc = perc
        hltwb_d(isd)%et = aet
        hltwb_d(isd)%tloss = 0.
        hltwb_d(isd)%eplant = 0.
        hltwb_d(isd)%esoil = 0.
        hltwb_d(isd)%surq_cont = 0.
        hltwb_d(isd)%cn = cn_sd
        hltwb_d(isd)%sw = hlt(isd)%sw
        hltwb_d(isd)%sw_final = hlt(isd)%sw
        hltwb_d(isd)%sw_300 = hlt(isd)%sw / 300.
        hltwb_d(isd)%snopack = hlt(isd)%snow
        hltwb_d(isd)%pet = pet
        hltwb_d(isd)%qtile = flow_tile
        hltwb_d(isd)%irr = air

!    output_nutbal - no nutrients currently in SWAT-DEG

!    output_losses - SWAT-DEG
        hltls_d(isd)%sedyld = sedin / (100. * hlt_db(ihlt_db)%dakm2) !! sedyld(isd) / hru_ha(isd)
        hltls_d(isd)%sedorgn = 0.   !! sedorgn(isd)
        hltls_d(isd)%sedorgp = 0.   !! sedorgp(isd)
        hltls_d(isd)%surqno3 = 0.   !! surqno3(isd)
        hltls_d(isd)%latno3 = 0.    !! latno3(isd)
        hltls_d(isd)%surqsolp = 0.  !! surqsolp(isd)
        hltls_d(isd)%usle = 0.      !! usle
        hltls_d(isd)%sedmin = 0.    !! sedminpa(isd) + sedminps(isd)
        hltls_d(isd)%tileno3 = 0.   !! tileno3(isd)

!    output_plantweather - SWAT-DEG
        hltpw_d(isd)%lai =  hlt(isd)%alai     !! lai
        hltpw_d(isd)%bioms =  hlt(isd)%dm     !! total biomass
        hltpw_d(isd)%yield =  yield          !! crop yield
        hltpw_d(isd)%residue =  0.           !! residue
        hltpw_d(isd)%sol_tmp =  0.           !! soil(isd)%phys(2))%tmp
        hltpw_d(isd)%strsw = 1. - ws         !! (1.-strsw_av(isd))
        hltpw_d(isd)%strsa = 1. - strsair    !! (1.-strsw_av(isd))
        hltpw_d(isd)%strstmp = 1. - tstress  !! (1.-strstmp_av)
        hltpw_d(isd)%strsn = 0.              !! (1.-strsn_av)
        hltpw_d(isd)%strsp = 0.              !! (1.-strsp_av)
        hltpw_d(isd)%nplnt = 0.              !! nplnt(isd)
        hltpw_d(isd)%percn = 0.              !! percn(isd)
        hltpw_d(isd)%pplnt = 0.              !! pplnt(isd)
        hltpw_d(isd)%tmx = tmax              !! tmx(isd)
        hltpw_d(isd)%tmn = tmin              !! tmn(isd)
        hltpw_d(isd)%tmpav = tave            !! tmpav(isd)
        hltpw_d(isd)%solrad = raobs          !! hru_ra(isd)
        hltpw_d(isd)%wndspd = wst(iwst)%weat%windsp         !! windspeed(isd)
        hltpw_d(isd)%rhum = wst(iwst)%weat%rhum             !! relative humidity(isd)
        hltpw_d(isd)%phubase0 = wst(iwst)%weat%phubase0     !! base zero potential heat units

         !! set values for outflow hydrograph
         !! storage locations set to zero are not currently used
         ob(icmd)%peakrate = peakr
         ob(icmd)%hd(1)%temp = 5. + .75 * tave        !!wtmp
         ob(icmd)%hd(1)%flo = chflow * cnv            !!qdr m3/d
         ob(icmd)%hd(1)%sed = sedin                   !!sedyld
         ob(icmd)%hd(1)%orgn = 0.
         ob(icmd)%hd(1)%sedp = 0.
         ob(icmd)%hd(1)%no3 = 0.
         ob(icmd)%hd(1)%solp = 0.
         ob(icmd)%hd(1)%chla = 0.
         ob(icmd)%hd(1)%nh3 = 0.                         !! NH3
         ob(icmd)%hd(1)%no2 = 0.                         !! NO2
         ob(icmd)%hd(1)%cbod = 0.
         ob(icmd)%hd(1)%dox = 0.
         ob(icmd)%hd(1)%san = 0.                             !! det sand
         ob(icmd)%hd(1)%sil = 0.                             !! det silt
         ob(icmd)%hd(1)%cla = 0.                             !! det clay
         ob(icmd)%hd(1)%sag = 0.                             !! det sml ag
         ob(icmd)%hd(1)%lag = 0.                             !! det lrg ag

         !! set values for recharge hydrograph
         !ob(icmd)%hd(2)%flo = perc

       return

      end subroutine hru_lte_control