text:
C$Procedure      SCENCD ( Encode spacecraft clock )

      SUBROUTINE SCENCD ( SC, SCLKCH, SCLKDP )

C$ Abstract
C
C     Encode character representation of spacecraft clock time into a
C     double precision number.
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     SCLK
C
C$ Keywords
C
C     CONVERSION
C     TIME
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE               'sclk.inc'

      INTEGER               SC
      CHARACTER*(*)         SCLKCH
      DOUBLE PRECISION      SCLKDP

C$ Brief_I/O
C
C     Variable  I/O  Description
C     --------  ---  --------------------------------------------------
C     SC         I   NAIF spacecraft identification code.
C     SCLKCH     I   Character representation of a spacecraft clock.
C     SCLKDP     O   Encoded representation of the clock count.
C     MXPART     P   Maximum number of spacecraft clock partitions.
C
C$ Detailed_Input
C
C     SC         is the standard NAIF ID of the spacecraft whose clock's
C                time is being encoded.
C
C     SCLKCH     is the character representation of some spacecraft's
C                clock count.
C
C                SCLKCH will have the following general format:
C
C                             'pp/sclk_string', or just
C                                'sclk_string'
C
C                'pp' is an integer greater than or equal to one
C                and is called the partition number.
C
C                Each mission is divided into some number of partitions.
C                A new partition starts when the spacecraft clock
C                resets, either to zero, or to some other
C                value. Thus, the first partition for any mission
C                starts with launch, and ends with the first clock
C                reset. The second partition starts immediately when
C                the first stopped, and so on.
C
C                In order to be completely unambiguous about a
C                particular time, you need to specify a partition number
C                along with the standard clock string.
C
C                Information about when partitions occur for different
C                missions is contained in a spacecraft clock kernel
C                file, which needs to be loaded into the kernel pool,
C                using the routines CLPOOL and FURNSH.
C
C                The routine SCPART is used to read the partition
C                start and stop times, in encoded units of SCLK (called
C                "ticks" -- see SCLKDP below) from the kernel file.
C
C                If the partition number is included, it must be
C                separated from the rest of the string by a '/'.
C                Any number of spaces may separate the partition number,
C                the '/', and the rest of the clock string.
C
C
C                If the partition number is omitted, a default partition
C                will be assumed. The default partition is the lowest-
C                numbered partition that contains the given clock time.
C                If the clock time does not fall in any of the
C                partition boundaries then an error is signaled.
C
C
C                'sclk_string' is a spacecraft specific clock string.
C                Using Galileo as an example, the full format is
C
C                               wwwwwwww:xx:y:z
C
C                where z is a mod-8 counter (values 0-7) which
C                increments approximately once every 8 1/3 ms., y is a
C                mod-10 counter (values 0-9) which increments once
C                every time z turns over, i.e., approximately once every
C                66 2/3 ms., xx is a mod-91 (values 0-90) counter
C                which increments once every time y turns over, i.e.,
C                once every 2/3 seconds. wwwwwwww is the Real-Time Image
C                Count (RIM), which increments once every time xx turns
C                over, i.e., once every 60 2/3 seconds. The roll-over
C                expression for the RIM is 16777215, which corresponds
C                to approximately 32 years.
C
C                wwwwwwww, xx, y, and z are referred to interchangeably
C                as the fields or components of the spacecraft clock.
C                SCLK components may be separated by any of these
C                five characters: ' '  ':'  ','  '-'  '.'
C                Any number of spaces can separate the components and
C                the delimiters. The presence of the RIM component
C                is required. Successive components may be omitted, and
C                in such cases are assumed to represent zero values.
C
C                Values for the individual components may exceed the
C                maximum expected values. For instance, '0:0:0:9' is
C                an acceptable Galileo clock string, and will convert
C                to the same number of ticks as '0:0:1:1'.
C
C                Consecutive delimiters containing no intervening digits
C                are treated as if they delimit zero components.
C
C                Trailing zeros should always be included to match the
C                length of the counter.  For example, a Galileo clock
C                count of '25684.90' should not be represented as
C                '25684.9'.
C
C                Some spacecraft clock components have offset, or
C                starting, values different from zero.  For example,
C                with an offset value of 1, a mod 20 counter would
C                cycle from 1 to 20 instead of from 0 to 19.
C
C                See the SCLK required reading for a detailed
C                description of the Voyager and Mars Observer clock
C                formats.
C
C
C$ Detailed_Output
C
C     SCLKDP     is the double precision encoding of SCLKCH.
C
C                The encoding is such that order and proximity will be
C                preserved. That is, if t1, t2, and t3 are spacecraft
C                clock times, and t1*, t2*, and t3* are their encodings,
C                then if
C
C                              t1 < t2 < t3, and
C
C                t2 is closer to t1 than to t3, you will have the result
C                that
C
C                             t1* < t2* < t3*, and
C
C                t2* is closer to t1* than to t3*.
C
C                The units of encoded SCLK are "ticks since the start of
C                the mission", where a "tick" is defined to be the
C                shortest time increment expressible by a particular
C                spacecraft's clock.
C
C                Each clock string without partition number represents
C                a certain number of ticks, but you need to include
C                partition information to determine the relative
C                position of that time in relation to the start of the
C                mission.
C
C                Since the end time of one partition is coincident
C                with the begin time of the next, there are two
C                different representations for this instant, and they
C                will both yield the same encoding.
C
C                For example, if partition 1 has an end time of t1, and
C                partition 2 has a begin time of t2, then if we did
C
C                   CALL SCENCD ( '1/t1', SC, X ) and
C                   CALL SCENCD ( '2/t2', SC, Y ), then
C
C                                  X = Y.
C
C                The individual routines TIKSnn, where nn is the
C                clock type code, contain more detailed information
C                on the conversion process.
C
C$ Parameters
C
C     MXPART     is the maximum number of spacecraft clock partitions
C                expected in the kernel file for any one spacecraft.
C                See the INCLUDE file sclk.inc for this parameter's
C                value.
C
C$ Exceptions
C
C     1) If the number of partitions in the kernel file for spacecraft
C        SC exceeds the parameter MXPART, the error
C        'SPICE(TOOMANYPARTS)' is signaled.
C
C
C     If a partition number is included in the SCLK string, the
C     following exceptions may occur:
C
C     2) If the partition number cannot be parsed as an integer, the
C        error 'SPICE(BADPARTNUMBER)' is signaled.
C
C     3) If the partition number is not in the range of the number of
C        partitions found in the kernel pool, the error
C        'SPICE(BADPARTNUMBER)' is signaled.
C
C     4) If the clock count does not fall in the boundaries of the
C        specified partition, the error 'SPICE(NOTINPART)' is
C        signaled.
C
C
C     If a partition number is not included in the SCLK string, the
C     following exception may occur.
C
C     5) If the clock count does not fall in the boundaries of any
C        partition found in the kernel pool, the error
C        'SPICE(NOPARTITION)' is signaled.
C
C     The following error is signaled by a routine called by SCENCD
C
C     6)  If any of the extracted clock components cannot be parsed as
C         integers, or the string has too many components, or the value
C         of one of the components is less than the offset value, then
C         the error SPICE(INVALIDSCLKSTRING) is signaled.
C
C$ Files
C
C     A kernel file containing spacecraft clock partition information
C     for the desired spacecraft must be loaded, using the routines
C     CLPOOL and FURNSH, before calling this routine.
C
C$ Particulars
C
C     In general, it is difficult to compare spacecraft clock counts
C     numerically since there are too many clock components for a
C     single comparison.  This routine provides a method of assigning a
C     single double precision number to a spacecraft's clock count,
C     given one of its character representations.
C
C     The routine SCDECD performs the inverse operation to SCENCD,
C     converting an encoded double precision number to character format.
C
C     To convert the string to ticks since the start of the mission,
C     SCENCD
C
C        1) Converts the non-partition portion of the string to
C           ticks, using the routine SCTIKS.
C
C        2) Determines the partition number for the clock time,
C           either by getting it directly from the input string, or
C           determining the default partition if none was specified.
C
C        3) Includes partition start and stop times, which are also
C           measured in ticks, to compute the number of ticks
C           since the beginning of the mission of the clock time.
C
C$ Examples
C
C      Double precision encodings of spacecraft clock counts are used to
C      tag pointing data in the C-kernel.
C
C      In the following example, pointing for a sequence of images from
C      the Voyager 2 narrow angle camera is requested from the C-kernel
C      using an array of character spacecraft clock counts as input.
C      The clock counts attached to the output are then decoded to
C      character and compared with the input strings.
C
C            CHARACTER*(25)     SCLKIN   ( 4 )
C            CHARACTER*(25)     SCLKOUT
C            CHARACTER*(25)     CLKTOL
C
C            DOUBLE PRECISION   TIMEIN
C            DOUBLE PRECISION   TIMOUT
C            DOUBLE PRECISION   CMAT     ( 3, 3 )
C
C            INTEGER            NPICS
C            INTEGER            SC
C
C            DATA  NPICS     /  4                   /
C
C            DATA  SCLKIN    / '2 / 20538:39:768',
C           .                  '2 / 20543:21:768',
C           .                  '2 / 20550:37',
C           .                  '2 / 20561:59'       /
C
C            DATA  CLKTOL   /  '      0:01:000'     /
C
C      C
C      C     The instrument we want pointing for is the Voyager 2
C      C     narrow angle camera.  The reference frame we want is
C      C     J2000. The spacecraft is Voyager 2.
C      C
C            INST = -32001
C            REF  = 'J2000'
C            SC   = -32
C
C      C
C      C     Load the appropriate files. We need
C      C
C      C     1) CK file containing pointing data.
C      C     2) Spacecraft clock kernel file, for SCENCD and SCDECD.
C      C
C            CALL CKLPF  ( 'VGR2NA.CK' )
C            CALL CLPOOL
C            CALL FURNSH ( 'SCLK.KER'  )
C
C      C
C      C     Convert the tolerance string to ticks.
C      C
C            CALL SCTIKS ( SC, CLKTOL, TOL )
C
C            DO I = 1, NPICS
C
C               CALL SCENCD ( SC, SCLKIN( I ), TIMEIN )
C
C               CALL CKGP   ( INST, TIMEIN, TOL, REF, CMAT, TIMOUT,
C           .                 FOUND )
C
C               CALL SCDECD ( SC, TIMOUT, SCLKOUT )
C
C               WRITE (*,*)
C               WRITE (*,*) 'Input  s/c clock count: ', SCLKIN( I )
C               WRITE (*,*) 'Output s/c clock count: ', SCLKOUT
C               WRITE (*,*) 'Output C-Matrix:        ', CMAT
C               WRITE (*,*)
C
C            END DO
C
C     The output from such a program might look like:
C
C
C            Input  s/c clock count:  2 / 20538:39:768
C            Output s/c clock count:  2/20538:39:768
C            Output C-Matrix:  'first C-matrix'
C
C            Input  s/c clock count:  2 / 20543:21:768
C            Output s/c clock count:  2/20543:22:768
C            Output C-Matrix:  'second C-matrix'
C
C            Input  s/c clock count:  2 / 20550:37
C            Output s/c clock count:  2/20550:36:768
C            Output C-Matrix:  'third C-matrix'
C
C            Input  s/c clock count:  2 / 20561:59
C            Output s/c clock count:  2/20561:58:768
C            Output C-Matrix:  'fourth C-matrix'
C
C$ Restrictions
C
C     None.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman (JPL)
C     J.M. Lynch   (JPL)
C     R.E. Thurman (JPL)
C
C$ Version
C
C-    SPICELIB Version 1.2.0, 28-FEB-2014 (BVS)
C
C        Added FAILED checks to prevent passing uninitialized values to
C        ANINT, which can causing numeric exceptions on some
C        environments.
C
C-    SPICELIB Version 1.1.0, 05-FEB-2008 (NJB)
C
C        The values of the parameter MXPART is now
C        provided by the INCLUDE file sclk.inc.
C
C-    SPICELIB Version 1.0.2, 22-AUG-2006 (EDW)
C
C        Replaced references to LDPOOL with references
C        to FURNSH.
C
C-    SPICELIB Version 1.0.1, 10-MAR-1992 (WLT)
C
C        Comment section for permuted index source lines was added
C        following the header.
C
C-    SPICELIB Version 1.0.0, 03-SEP-1990 (JML) (RET)
C
C-&

C$ Index_Entries
C
C     encode spacecraft_clock
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               FAILED
      LOGICAL               RETURN
      INTEGER               CPOS

C
C     Local variables
C
      INTEGER               POS
      INTEGER               NPARTS
      INTEGER               PNTER
      INTEGER               PART
      INTEGER               I

      DOUBLE PRECISION      TICKS
      DOUBLE PRECISION      PSTART   ( MXPART )
      DOUBLE PRECISION      PSTOP    ( MXPART )
      DOUBLE PRECISION      PTOTLS   ( MXPART )

      CHARACTER*(25)        ERROR

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

C
C     Convert the non-partition portion of the clock string to ticks.
C
      POS = CPOS ( SCLKCH, '/', 1 )

      CALL SCTIKS (  SC, SCLKCH( POS+1 : ), TICKS )

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

      TICKS = ANINT ( TICKS )

C
C     Read the partition start and stop times (in ticks) for this
C     mission. Error if there are too many of them.
C
      CALL SCPART ( SC, NPARTS, PSTART, PSTOP )

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

      IF ( NPARTS .GT. MXPART ) THEN

         CALL SETMSG ( 'The number of partitions, #, for spacecraft ' //
     .                 '# exceeds the value for parameter MXPART, #.' )
         CALL ERRINT ( '#', NPARTS           )
         CALL ERRINT ( '#', SC               )
         CALL ERRINT ( '#', MXPART           )
         CALL SIGERR ( 'SPICE(TOOMANYPARTS)' )
         CALL CHKOUT ( 'SCENCD'              )
         RETURN

      END IF

C
C     PSTART and PSTOP represent integers but are read from the
C     kernel pool as double precision numbers. Make them whole
C     numbers so that logical tests may be performed with them.
C
      DO I = 1, NPARTS
         PSTOP (I) = ANINT ( PSTOP(I)  )
         PSTART(I) = ANINT ( PSTART(I) )
      END DO

C     For each partition, compute the total number of ticks in that
C     partition plus all preceding partitions.
C
      PTOTLS( 1 ) = ANINT ( PSTOP( 1 ) - PSTART( 1 ) )

      DO I = 2, NPARTS
         PTOTLS( I ) = ANINT ( PTOTLS( I-1 ) + PSTOP( I ) - PSTART( I ))
      END DO

C
C     Determine the partition number for the input clock string:
C
C        If it was included in the string make sure it's valid for
C        this mission.
C
C           Error if
C
C           1) The partition number can't be parsed.
C           2) The partition number is not in the range 1 to the number
C              of partitions.
C           3) The clock count does not fall in the boundaries of the
C              specified partition.
C
C        If it wasn't included, determine the default partition for
C        this clock count.
C
C           Error if
C
C           1) The clock count does not fall in the boundaries of any
C              of the partitions.
C
C
      IF ( POS .EQ. 1 ) THEN

         CALL SETMSG ( 'Unable to parse the partition number '   //
     .                 'from SCLK string #.'  )
         CALL ERRCH  ( '#', SCLKCH            )
         CALL SIGERR ( 'SPICE(BADPARTNUMBER)' )
         CALL CHKOUT ( 'SCENCD'               )
          RETURN
      END IF

      IF ( POS .GT. 1 ) THEN

         PART = 0

         CALL NPARSI ( SCLKCH( :POS-1 ), PART, ERROR, PNTER )

         IF ( ERROR .NE. ' ' ) THEN

            CALL SETMSG ( 'Unable to parse the partition number '   //
     .                    'from SCLK string #.'  )
            CALL ERRCH  ( '#', SCLKCH            )
            CALL SIGERR ( 'SPICE(BADPARTNUMBER)' )
            CALL CHKOUT ( 'SCENCD'               )
            RETURN

         ELSE IF ( PART .LE. 0 .OR. PART .GT. NPARTS ) THEN

            CALL SETMSG ( 'Partition number # taken from SCLK string '//
     .                    '# is not in acceptable range 1 to #.'     )
            CALL ERRINT ( '#', PART              )
            CALL ERRCH  ( '#', SCLKCH            )
            CALL ERRINT ( '#', NPARTS            )
            CALL SIGERR ( 'SPICE(BADPARTNUMBER)' )
            CALL CHKOUT ( 'SCENCD'               )
            RETURN

         ELSE IF ( TICKS .LT. PSTART( PART )   .OR.
     .             TICKS .GT. PSTOP ( PART ) )  THEN

            CALL SETMSG ( 'SCLK count # does not fall in the '        //
     .                    'boundaries of partition number #.'  )
            CALL ERRCH  ( '#', SCLKCH        )
            CALL ERRINT ( '#', PART          )
            CALL SIGERR ( 'SPICE(NOTINPART)' )
            CALL CHKOUT ( 'SCENCD'           )
            RETURN

         END IF

      ELSE

         PART = 1

         DO WHILE (   PART  .LE. NPARTS             .AND.
     .              ( TICKS .LT. PSTART( PART )     .OR.
     .                TICKS .GT. PSTOP ( PART ) )        )

            PART = PART + 1

         END DO

         IF ( PART .GT. NPARTS ) THEN

            CALL SETMSG ( 'SCLK count # does not fall in the '        //
     .                    'boundaries of any of the partitions for '  //
     .                    'spacecraft #.'      )
            CALL ERRCH  ( '#', SCLKCH          )
            CALL ERRINT ( '#', SC              )
            CALL SIGERR ( 'SPICE(NOPARTITION)' )
            CALL CHKOUT ( 'SCENCD'             )
            RETURN

         END IF

      END IF

C
C     Now we have a valid partition number, and the number of ticks for
C     the clock string. To convert to ticks since the start of the
C     mission, add in the total number of ticks in preceding partitions
C     and subtract off the starting ticks value for this partition.
C
      IF ( PART .GT. 1 ) THEN
         SCLKDP = TICKS - PSTART( PART ) + PTOTLS( PART - 1 )
      ELSE
         SCLKDP = TICKS - PSTART( PART )
      END IF

      CALL CHKOUT ( 'SCENCD' )
      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:
        SUBROUTINE  MA_S4DG ( marrpt, ipt, iret )
C************************************************************************
C* MA_S4ED                                                              *
C*                                                                      *
C* This subroutine decodes the drifting buoy section 4 group            *
C* 9i(d)Z(d)Z(d)Z(d), which contains the drogue type and the length of  *
C* the cable at which the drogue is attached.                           *
C*                                                                      *
C* MA_S4ED ( MARRPT, IPT, IRET )                                        *
C*                                                                      *
C* Input parameters:                                                    *
C*      MARRPT         CHAR*             Report array                   *
C*	IPT	       INTEGER	       	 Points to the space preceding  *
C*                                       the group 9iZZZZ               *
C* Output parameters:                                                   *
C*      RIVALS(IRDROT) REAL              Drogue type                    *
C*      RIVALS(IRCALN) REAL              Cable length in meters         *
C*      IRET           INTEGER           Return code                    *
C*                                         0 = normal return            *
C*                                         1 = problems                 *
C*                                                                      *
C**                                                                     *
C* Log:                                                                 *
C* R. Hollern/NCEP       9/99                                           *
C************************************************************************
        INCLUDE         'macmn.cmn'
C*
        CHARACTER*(*)   marrpt
C*
        CHARACTER       fld1*1,  fld3*3
C------------------------------------------------------------------------
        iret  = 0
C
        ipt  = ipt + 2
        fld1 = marrpt ( ipt:ipt )
        CALL  ST_INTG ( fld1, ival, ier )
        IF ( ier .eq. 0 ) THEN
            rivals ( irdrot )  =  FLOAT ( ival )
          ELSE
            iret = 1
            RETURN
        END IF
C
        ipt  = ipt + 1
        fld3 = marrpt ( ipt:ipt+2 )
        CALL  ST_INTG ( fld3, ival, ier )
        IF ( ier .eq. 0 ) THEN
            rivals ( ircaln )  =  FLOAT ( ival )
        END IF
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Sacramento Valley Historical Railways is a nonprofit, allvolunteer organizations organization dedicated to educating the public about railroads and railroad history and to the preservation and operation of historic railroad equipment and structures in Northern California. Located at the former Southern Pacific Depot in Woodland, the men who tirelessly work there to preserve the Southern Pacific 060 SixCoupled Switching Locomotive, the Southern Pacific Caboose 1156, and former Spreckels Sugar #1, exYolo Shortline #50 50 ton General Electric Switcher, are friendly, open to talks about the railroads, and more.

1233

Steam locomotive 1233 was built in July, 1918, by the Baldwin Locomotive Works of Eddystone, Pennsylvania as part of an order of 15 060 switch engines. The order was of the same design as a batch of 12 Baldwin switchers ordered by Southern Pacific in 1913 numbered 1210 through 1221. Of the order which included the 1233, five were ordered for the Central Pacific and one was ordered for the Arizona Eastern. The latter was designated by Southern Pacific as Class S51 and the remainder were designated Class S10. Once delivered, they became part of the more than 450 sixwheeled steam locomotives owned by the Southern Pacific system during the 100plus years that steam reigned as king on the railroads.

Weighing in at 154,600 pounds on its 51 inch drivers, the 1233 occupied 30 feet, eight and onequarter inches of rail and stood 13 feet, eight inches at its highest point. Oilfired, the 1233s boiler operated at a maximum of 175 psi and developed a tractive effort of 27,376 pounds. The locomotives in 1233s buildorder were delivered with acetylene headlamps and oilfired lamps in the cab. Each engine cost $35,000.

The 1233 lived an unremarkable life, switching freight cars in the yards at Tracy and Oakland for most of its career. About six months after it was delivered, the 1233 underwent several boiler modifications at the Tracy shops. When it was returned to service on February 18, 1919, the 1233 had an increased maximum boiler pressure of 190 psi and its tractive effort had increased by 2,350 pounds.

In 1926, the entire firebox was replaced, although the reason was not noted on the record. In 1934, the 1233 was converted (at Oakland) to operate with electric lights. A dynamo was installed on the boiler just ahead of the cab, the carbon tank was removed from the firemans cab floor, and the oil lamps in the cab and acetylene headlamp were all replaced with a 32 volt electrical system. Sometime in 1936, a new crown sheet was installed. The 1233 had its last major overhaul in October 1951, when it received new journal brasses and tires, the boiler flues were reset, and the locomotive received a thorough goingover.

Near the end of its service, the 1233 was used to shuttle passenger cars through the wash rack at Oakland. In March of 1956, the 1233 was retired from service. Alco S2 dieselelectric locomotives had been purchased to replace the aging steamers and the 1233 joined many of its sisters in the dead line at Oakland. Its date with the scrappers torch was apparently sealed when it was vacated from the company roster on July 24, 1956.

At about the same time, the Woodland Chamber of Commerce was looking for a steam locomotive to display at the fairgrounds. In the mid1950s, the glut of deadlined steam locomotives had caused the price of steel and iron to drop to nearly nothing. Railroads, rather than hold the locomotives and wait for the market to improve, wanted to move the locomotives off their property as soon as possible. Like many other railroads, Southern Pacific offered steam locomotives to cities for display. And Woodland, like many other cities, wanted a steam locomotive, but not a big one...just a small reminder of the steam era.

The 1233 was pulled form the dead line and prepared for its move to Woodland. The tender retired behind the 1233, a Class 52C1, was apparently deemed unfit to make the journey, so the tender was recovered from another locomotive. In a photo provided by rail photographer Stan Kistler, the 1233 sits in Oakland on September 9, 1956, with a fresh coat of paint and the replacement tender, a class 70C10. This tender was apparently also not roadworthy, as it was replaced prior to the trip with tender #6970, another Class 70C10. In an ironic twist, the 6970 was retired behind the 1234!

Final arrangements were made and the 1233 moved in a train from Oakland to Woodland near the end of December, 1956. Interestingly enough, when the locomotive arrived in Davis, it was deemed to be facing the wrong way. According to the Woodland Daily Democrat newspaper, some debate had ensued days earlier as to whether the 1233 should be displayed facing Highway 113 from the fairgrounds or show its hind end to all who passed. By the time it was decided to point the locomotive to the west, it was already en route. So, upon arrival at Davis, the locomotive was run around the wye which surrounds the depot there, then sent on to Woodland. It was stored on the Adams Grain siding across the highway from the fairgrounds.

Just how the locomotive would be moved across the street was the source of much discussion. In fact, the Chamber of Commerce sponsored a contest to solicit ideas. Winning the $10 first prize was Paul Oxley of Woodland, who suggested the obvious: lay track across the street and roll it into the fairgrounds.

On Saturday, March 30, 1957, work began to move the 1233 across East Street into the Yolo County Fairgrounds. Southern Pacific Roadmaster Frank Pucci directed the track laying and supervised the move. A local trucking company provided two truck tractors. An abundance of help showed up, and the entire scene took on a carnivallike atmosphere. The move took the entire weekend and was not without incident. Close examination of the tender front deck and locomotive apron will reveal the scars of a hasty torch cutting. Apparently the track curvature into the fairgrounds was sharper than the minimum radius the 1233 could negotiate, so the offending parts had to be temporarily removed.

Once the move was completed, small wedges were welded on the rails ahead and behind one set of locomotive wheels. A dedication was held during the Yolo County Fair in August of that year.

The 1233 was moved about five years later when East Street was widened to four lanes. Track was laid and the 1233 was pushed back about 50 feet. Over the years, it suffered the indignation of serving as a jungle gym for kids, a target for rocktoting delinquents, and a habitat for squirrels, moss and rust. It was eventually fenced, and remained so until June of 1982 when initial restoration work began.

Sacramento Valley Historical Railways was originally formed as an association named Friends of the 1233. The founders and other volunteers began the initial cosmetic restoration of the 1233 in 1982. That year, the 1233 was painted for the Yolo County Fair. After the fair, the group decided that the 1233 could be made operational and went to work. In the succeeding years, the group incorporated as Sacramento Valley Historical Railways and obtained title to the 1233. SVHR also acquired other railroad equipment and the Southern Pacific Depot in Woodland.

Progress on the restoration of the 1233 was made each year and in the spring of 1989 it ran for the first time since the mid 1950s on a short piece of panel track built at the Yolo County Fairgrounds. However, it was apparent that further restoration and a complete boiler tube job would be required before the 1233 could be operated reliably. Also, the panel track had to be taken up for the fair, and the fairgrounds would not be a place where the 1233 could be operated in any significant manner.

After much discussion followed by negotiations with the (then) Yolo Shortline Railroad, it was decided that the best opportunity to operate the 1233 was to place it into excursion operations on the Yolo Shortline. While this resulted in some departures from perfect historically accurate preservation (e.g. the 1233 was lettered for the Yolo Shortline), the maintenance and care would be increased, and the 1233 would have a far greater public visibility.

Other Notes

The firebox on the 1233 was replaced again much later than 1936. We think that most of the firebox was replaced in the 1951 shopping. Most of the firebox sheets are full thickness and the stays look new, except for being encrusted with scale (probably due to lack of boiler treatment in the last days of SP steam operations). There is no wearing of welds, stay bolt threads, peened edges, etc., inside the firebox, which would occur through the normal course of use and firing of the boiler.

A power reverser replaced the original Johnson Bar. Federal regulations required this, we believe in the 1930s.

The current tender behind the 1233 is an ALCO tender, not Baldwin built. It started life, we think, behind SP consolidations. It is at least the fourth tender to operate behind the 1233 (we have pictures showing three other tenders during the 1233s working life). The 1233s original tender was a sausage type of tender. The number on the back of the current tender, which was uncovered during repainting, was the 1267, not the 1234. Perhaps it saw service behind both of those engines.

In January 1995, the Yolo Shortline trucked the 1233 from the fairgrounds approximately 1 mile to its yard on East Main Street. The engine was trucked separately from the tender. It rained on both days on which the moves occured. Temporary panel and ramp tracks were built to load the locomotive and tender onto the low bed trailers. Each move took approximately 45 hours from start to finish.

The planned work on the 1233 to place it into operation did not get started until early winter of 19951996. All of the 1233s new tubes were reset and two superheater flues were replaced. The FRA inspections were made and a successful hydro test was accomplished. The crosshead guides were remachined and remounted and new crosshead slippers were made. The rods were put on and a new boiler jacket was fabricated. Many other small items were attended to, such as fixing the airpump, getting piston rings unstuck (from carbonization) and installing the smokebox appliances. The 1233 was painted and lettered in a style patterned after the pre1940s SP style.

The first fireup of the 1233 occured in April of 1996 and several test runs were made in the Woodland yards to check new bearing and fix minor problems that arose. The first steam trip was made with a full excursion train in May when the 1233 ran from Woodland to Clarksburg and back. Most people riding that first trip never new that the 1233 broke its main water line to the boiler injectors that day. It happened while the 1233 and its train was awaiting the return of its passengers from a lunch break. The crew, with few tools, managed to make a hasty fix.

While in operation, the 1233 proved capable for both freight and passenger service. The 1233 saw some freight service, since it was going that way anyway, so why not take the freight too? The primary operational problem was foaming in the boiler and priming (carrying water over into the throttle and to the cylinders). This was due primarily to old scale remaining in the boiler from SP days. Frequent boiler washes helped the problem and water treatment proved essential. The cylinder lubricant was also changed to help alleviate the problems caused by the foaming and priming.

Woodland Depot

The California Pacific Railroad, a predecessor of the Southern Pacific Railroad, built the first railroad depot in Woodland located between Main Street and Oak Street. It was a two story, box shaped structure with “living apartments” on the second floor for the station agent. It had a large freight house, but a small waiting room for passengers.

Realizing that such a structure was unsuitable for the growing city of Woodland, the Chamber of Commerce, on February 2, 1910, sent three of its members to San Francisco to talk to Southern Pacific officials about a new depot. When they returned to Woodland on the midnight train, they were able to report that a new passenger depot would be built, but they had no firm date as to when that event would occur. Rumors on the status of the new depot persisted until February 10, 1911, when Division Superintendent Sheridan arrived in Woodland on a special train. He informed Station Agent Fingland that General Manager’s Order No. 2 had been issued authorizing the construction of the new depot. The recent news that the Northern Electric Railway was building a rail line into Woodland with a new depot at Second and Main Streets probably hastened this decision.

After some controversy, the City of Woodland allowed the Southern Pacific Railroad to locate the new depot south of the old depot. The new location straddled Lincoln Avenue and forced the closure of Lincoln Avenue to East Street but allowed long trains to stop in front of the new depot without having to uncouple at Main and Oak Streets for roadway traffic. Construction of the new Woodland Depot began and it was completed in 1911. The original blueprints, dated February, 1910, were revised on February 13, 1911. An additional blueprint for the cement work in the waiting rooms and restrooms is dated April 21, 1911. A Sanborn Insurance map dated January, 1912, shows the completed depot and the relocation of the old depot to a spot south of Lincoln Avenue and east of Sixth Street. The old depot was then designated as a freight only facility.

Constructed mostly of prime redwood and Douglas fir, the new depot was equipped with “every modern accommodation” for its time. It was a standard Southern Pacific colonnade style design with Colonial Revival and Craftsman elements. The 156’ long arcade on the east side was constructed to protect the passengers from the weather and was 40’ longer than the old depot to accommodate longer trains. It had a baggage room, station agent’s office, restrooms, and separate waiting rooms for men and women. The women’s waiting room provided a private room for unaccompanied women travelers and a place where mothers could nurse their children.

During the early years of its operation, this depot was heavily patronized and was a stop for many scheduled trains on the West Valley Line. Historic photos regularly show large crowds in front of the depot. The depot was also witness to several historic events. On May 20, 1942, 418 Yolo County citizens of Japanese ancestry were required to board a train to the Merced Assembly Center. The majority were then sent on to the Granada Relocation Center in Colorado for the remainder of World War II. Famed Depression era photographer Dorothea Lange captured this emotional moment in a series of black and white photographs, one of which is shown here. The photos also show many of the architectural elements of the depot at that time.

In the early 1950s, Southern Pacific passenger service to Woodland declined and was eventually suspended. In 1955, the Greyhound Bus Company leased the two waiting rooms. They installed a ticket counter on the south wall and it is believed that they removed the wall between the two waiting rooms and put double doors on the north and west walls at this time. A Greyhound Bus sign was put on the north end of the roof.

There were additional changes; In 1956, the depot was painted “an attractive shade of seafoam green” replacing its original colonial yellow paint scheme which was no longer popular with local residents. The next year, freight operations were moved from the old freight depot to the baggage room of the passenger depot. A door was installed between the agent’s office and the baggage room and the trunk platform was removed from the west wall of the baggage room. The sliding door above the trunk platform was raised to accommodate a truck loading dock that was installed in the western two thirds of the baggage room and extended through the doorway. The screen partitions were also removed from the counter in the office. Later, a restroom for train crews was constructed in the southeast corner of the baggage room accessed through a new exterior door in the south wall. By this time, the only railroad employee remaining at the depot was the local freight agent.

By the 1980’s, the depot had fallen into disrepair and had not had a passenger train stop there for many years. On April 24, 1982, the last regularly scheduled passenger train through Woodland, Amtrak’s Train 11the Starlight, rolled past the depot without stopping. The worst news came in 1988, when Southern Pacific announced plans to raze the depot. Having just installed a new computer system, Southern Pacific decided to transfer the station agent to Sacramento, where he could accomplish all of the Woodland work by telephone and computer.

With that announcement, a group of local preservationists and rail fans went before the Woodland City Council. The Council was encouraged to try and save the depot, but cautioned that no local taxpayers’ money could be expected. A Woodland S.P. Depot Committee was then formed which consisted of representatives of Sacramento Valley Historical Railways, the Yolo County Historical Society, the City of Woodland, and members of the general public. Sacramento Valley Historical Railways, a local group that was restoring Southern Pacific steam locomotive 1233 at the Yolo Fair Grounds, agreed to lead the preservation movement.

The committee approached the Southern Pacific Company and secured an agreement to have the depot donated to the Sacramento Valley Historical Railways with the provision that it be moved away from the railroads main line. The railroad also agreed to sell its property located about 200 feet southwest of the depot to Sacramento Valley Historical Railways so that the depot could be relocated there. The original depot, converted to a freight depot, had been moved to this location but was burned down by transients in the 1970s. The committee applied for a Community Development Block Grant and raised matching funds to purchase the property. In 1990, the grant was approved, and the property was purchased from Southern Pacific in 1991. At the same time, the owner of the property adjacent to the new depot site at Sixth Street and Lincoln Avenue donated it to Sacramento Valley Historical Railways. A grain storage building that stood on this property was subsequently dismantled. This provided additional space to display locomotive 1233 and other full sized railroad equipment owned by Sacramento Valley Historical Railways.

After many numerous delays, the depot was moved during the summer of 1992 and set on its new foundation. The cost of the move was donated by a individual and many Woodland businesses made donations of material and labor for the foundation. A chain link fence was erected around the depot so that repair and restoration could begin. Since the depot was moved, volunteers from Sacramento Valley Historical Railways have spent hundreds of hours restoring it to what it looked like in the early 1920s.

Officers
 Also see http://www.svhr.org/contact/contact.htm contact page

  Mike Adams, President
  Al Van Hoosear, Vice President
  Mike Evans, Secretary
  Terry Schmidt, Treasurer
  Richard Buzz Jones, Railhead Editor

2012 Meetings

Board Meetings
 Jan 31
 Feb 28
 Apr 03
 May 01
 May 29
 Jul 03
 Jul 31
 Aug 28
 Oct 02
 Oct 30
 Nov 27
 Jan 08

Work Days
 Feb 04
 Mar 03
 Apr 07
 May 05
 Jun 02
 Jul 07
 Aug 04
 Sep 01
 Oct 06
 Nov 03
 Dec 01
 Jan 12

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine qqb_wgam_z(p,z)
      implicit none
c---Integrated Matrix element SUBTRACTION squared averaged
C---over initial colors and spins
c     q(-p1)+qbar(-p2) -->  W^+(nu(p3)+e+(p4))+a(p5)+g(p6)
c     q(-p1)+qbar(-p2) -->  W^-(e-(p3)+nubar(p4))+a(p5)+g(p6)
      include 'constants.f'
      include 'qcdcouple.f'
      include 'scale.f'
      include 'PR_new.f'
      include 'agq.f'
      integer is
      double precision z,xl12,p(mxpart,4),dot,ii_qq,ii_qg,tempqq,tempqg

      xl12=dlog(two*dot(p,1,2)/musq)

c----contributions for one leg

      do is=1,3
      tempqq=+ason2pi*cf*ii_qq(z,xl12,is)
      tempqg=+ason2pi*tr*ii_qg(z,xl12,is)

      Q1(q,q,a,is)=tempqq
      Q2(a,a,q,is)=tempqq
      Q1(a,a,q,is)=tempqq
      Q2(q,q,a,is)=tempqq

      Q2(q,g,q,is)=tempqg
      Q2(a,g,q,is)=tempqg
      Q2(q,g,a,is)=tempqg
      Q2(a,g,a,is)=tempqg
      Q1(q,g,q,is)=tempqg
      Q1(a,g,q,is)=tempqg
      Q1(q,g,a,is)=tempqg
      Q1(a,g,a,is)=tempqg
      enddo
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


module variables
  real(kind=8) :: mbh=0.0     ! central black hole mass
  real(kind=8) :: radius  ! radius
  real(kind=8) :: m200=0.0  ! virial mass
  real(kind=8) :: mvir=0.0  ! virial mass
  real(kind=8) :: mhalo=0.0 ! halo mass
  real(kind=8) :: v200=0.0  ! virial velocity
  integer(kind=4) :: halotype ! Type of halo
  real(kind=8) :: rs      ! NFW scale radius
  real(kind=8) :: r200    ! virial radius
  real(kind=8) :: c200    ! NFW concentration
  real(kind=8) :: sigma0  ! central surface density
  real(kind=8) :: rdisc=0.0   ! disc scale length
  real(kind=8) :: bdisc=0.0   ! disc shape parameter (Miyamoto-Nagai disc)
  real(kind=8) :: mbulge=0.0  ! bulge mass
  real(kind=8) :: fbulge=0.0  ! bulge extent
  real(kind=8) :: rbulge=0.0  ! bulge extent
  real(kind=8) :: mstar=0.0   ! stellar mass
  real(kind=8) :: mdisc=0.0   ! disc mass
  real(kind=8) :: fdisc=0.0   ! stellar mass
  real(kind=8) :: beta0
  real(kind=8) :: r_a

  integer(kind=4) :: nrbin
  real(kind=8), allocatable :: rtab(:)
  real(kind=8), allocatable :: sigma_r(:)
  real(kind=8), allocatable :: sigma_p(:)
  real(kind=8), allocatable :: sigma_m(:)
end module variables


text:
program test
  real :: f = 1.2
  real(kind=8) :: d = 3.4

  print *, int(f)
  print *, int(d)
end program test

max_stars_repo_name: clementval/fc

max_issues_repo_name: clementval/fc

max_forks_repo_name: clementval/fc


text:
      program t
      integer::i(4)
      i(1)=1
      i(2)=2
      i(3)=3
      i(4)=4
      print '(i0)',i(2:3)
      end program t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


module variables
  implicit none
  integer n_el, sp_deg, n_bas, n_weyl
  real*8  nuc_repul
  real*8,  dimension(:,:), allocatable  :: S, H
  real*8,  dimension(:,:,:,:), allocatable  :: EE
  integer AllocateStatus
end module


text:
LOGICAL :: L = .FALSE.

SELECT CASE (L)
   CASE (.TRUE.)
      CALL abort
   CASE (.FALSE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

SELECT CASE (L)
   CASE (.TRUE., .FALSE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

SELECT CASE (L)
   CASE (.FALSE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

SELECT CASE (L)
   CASE (.NOT. .TRUE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

SELECT CASE (.NOT. L)
   CASE (.TRUE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

SELECT CASE (Truth_or_Dare() .OR. L)
   CASE (.TRUE.)
      CONTINUE
   CASE DEFAULT
      CALL abort
END SELECT

CONTAINS

   FUNCTION Truth_or_Dare ()
      LOGICAL Truth_or_Dare
      Truth_or_Dare = .TRUE.
   END FUNCTION

END

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
program t
end program


text:
Describe Users/PLLC here.

20090208 17:41:34 nbsp Welcome to the Wiki. Please dont Creative Hair delete pages or peoples comments. If you have a comment or minor change for the page, please do so, but eliminating the entire page without explaining why is disruptive. Users/robinlaughlin

20090208 17:43:40 nbsp Howdy! I keep reverting your edits because the page is already doing fine on its own and doesnt need to be remade. By the way, you share the same IP address as Users/KC, which is pretty solid evidence that you and they are the same person. That sort of thing is generally frowned on around here, because pretending to be more than one person is kinda unethical. Hope all this helps you understand why there have been reverts and such on the Creative Hair page. Users/JoePomidor

20090208 17:45:56 nbsp Ah, also, I have been making sure that the comments put on the page have the proper attribution of who said them. This is important since the wiki allows people to change pretty much any content on a page, and so the only real stamp of veracity is ones name. This means a comment without a name means nothing, and so will either be attributed to the person who actually said it, or deleted. Users/JoePomidor

20090208 17:48:12 nbsp It looks like KC and PLLC were both editing from the same IP address. It could be that they are both the same pseudonymous person? If that is the case they may just be removing their own comment. Situations like this can be avoided by Importance of using your RealName using your real name. Users/JasonAller

20090208 18:00:36 nbsp Hallo again, The best summation I can give as far as comments go is, its ok to edit or delete your own, no matter what page. The big thing is just not to edit other peoples comments, as that tends to offend. However, it is OK to delete other peoples comments on your user page if you so desire. I deleted your comments off my page, since that seemed to be what you were asking for, but in the future you can do that sort of thing too by clicking the edit button at the top of any given page, next to the title. Users/JoePomidor


! Run with `-np 2` and then with `-np 3` or more.
!
! Demonstrate first what happens if you just use first and last: rank 0
! deadlocks, draw on the blackboard the diagram, as rank 1 sends to itself, but
! rank 0 tries to recieve from itself.
!
program deadlock
    use mpi_f08
    implicit none
    integer, parameter :: N = 5
    integer :: my_rank, n_ranks, prev, next, tag
    integer, allocatable :: a(:), b(:)
    type(MPI_Status) :: status

    call MPI_Init()
    call MPI_Comm_rank(MPI_COMM_WORLD, my_rank)
    call MPI_Comm_size(MPI_COMM_WORLD, n_ranks)

    prev = modulo(my_rank - 1, n_ranks)
    next = modulo(my_rank + 1, n_ranks)

    allocate(a(N), source=my_rank)
    allocate(b(N), source=-1)

    tag = 0

    print "(a, i0, a, *(1x, i2))", "before at rank=", my_rank, ", a(:)=", a(:)
    print "(a, i0, a, *(1x, i2))", "before at rank=", my_rank, ", b(:)=", b(:)

    call MPI_Sendrecv(                                                          &
        sendbuf=a, sendcount=N, sendtype=MPI_INTEGER, dest=next,   sendtag=tag, &
        recvbuf=b, recvcount=N, recvtype=MPI_INTEGER, source=prev, recvtag=tag, &
        comm=MPI_COMM_WORLD, status=status                                      )

    print "(a, i0, a, *(1x, i2))", "after  at rank=", my_rank, ", a(:)=", a(:)
    print "(a, i0, a, *(1x, i2))", "after  at rank=", my_rank, ", b(:)=", b(:)

    deallocate(a)
    deallocate(b)

    call MPI_Finalize()
end program deadlock


! { dg-do compile }
!
! PR fortran/44347 - arguments of SELECTED_REAL_KIND shall be scalar
! Testcase contributed by Vittorio Zecca <zeccav AT gmail DOT com>
!

  dimension ip(1), ir(1)
  i = selected_real_kind(ip, i)      ! { dg-error "must be a scalar" }
  j = selected_real_kind(i, ir)      ! { dg-error "must be a scalar" }
end


text:
      SUBROUTINE NASTIM (IHR, IMN, ISC, CPUSEC)
      REAL ARRAY(2)
      REAL RESULT
      CALL ETIME(ARRAY,RESULT)
      SECS   = ARRAY(2)
      IHR    = SECS / 3600.
      IMN    = ( SECS - 3600.*IHR ) / 60.
      ISC    = SECS - ( 3600.*IHR ) - ( 60.*IMN )
      CPUSEC = SECS
      RETURN
      END


text:
program da_verif_obs
   !---------------------------------------------------------------------------
   ! History:
   !
   !  Abstract:
   !  Program to read diagnostics written in fort.50 by WRFVAR
   !  and write in proper format to get ploted using PC-XL utility
   !
   !  Author:   Syed RH Rizvi     NCAR/MMM         05/25/2006
   !  Updates:
   !            Hui Shao          NCAR/RAL/DATC    05/02/2007
   !                      Diagnositics for GPSREF
   !            Syed RH Rizvi     NCAR/MMM         05/08/2007
   !            Significance test & error bars are added
   !---------------------------------------------------------------------------

   use da_verif_obs_control, only : surface_type, upr_type, gpspw_type, &
      gpsref_type, record1, record2, record3, &
      record4, record5, record6, stats_value, exp_dirs, out_dirs, nstd,nstdh, &
      rmiss, diag_unit_out, nml_unit, alpha, &
      diag_unit_in, info_unit, exp_num, end_date, file_path_string, &
      if_plot_bias, if_plot_airsret, if_plot_airep,if_plot_abias, &
      if_plot_buoy, if_plot_gpspw, if_plot_gpsref, if_plot_pilot, &
      if_plot_profiler, if_plot_polaramv, if_plot_qscat, if_plot_rmse, &
      if_plot_sound, if_plot_sonde_sfc, if_plot_synop, if_plot_surface, &
      if_plot_upr, if_plot_ships, if_plot_metar, if_plot_tamdar,interval, stdp, start_date, &
      if_plot_geoamv, stdh, num_miss, &
      wrf_file, istart, iend, jstart, jend
   use da_verif_obs_init, only : initialize_surface_type, initialize_upr_type, &
      initialize_gpspw_type, initialize_gpsref_type, da_advance_cymdh , &
      initialize_t_tab
   use da_verif_tools, only : map_info, proj_merc, proj_ps,proj_lc,proj_latlon, &
      da_llxy_wrf,da_xyll,da_map_set

   implicit none

   integer      :: num_obs
   character*20 :: obs_type, dummy_c

   character*5  :: stn_id
   integer      :: n, k, kk, l, levels, dummy_i
   real         :: lat, lon, press, height, dummy
   real         :: u_obs, u_inv, u_error, u_inc, &
                   v_obs, v_inv, v_error, v_inc, &
                   t_obs, t_inv, t_error, t_inc, &
                   p_obs, p_inv, p_error, p_inc, &
                   q_obs, q_inv, q_error, q_inc, &
                   spd_obs, spd_inv, spd_err, spd_inc
   real         :: tpw_obs, tpw_inv, tpw_err, tpw_inc
   real         :: ref_obs, ref_inv, ref_err, ref_inc
   integer      :: u_qc, v_qc, t_qc, p_qc, q_qc, tpw_qc, spd_qc, ref_qc
   integer      :: npr, nht, ier, iexp
   character*10 :: date, new_date             ! Current date (ccyymmddhh).
   integer      :: sdate, cdate, edate        ! Starting, current ending dates.
   logical      :: if_write, is_file
   logical, allocatable :: l_skip(:)

   character(len=512)     :: out_dir,filename
   type (surface_type)    :: surface
   type (upr_type)        :: upr, gupr
   type (gpspw_type)      :: gpspw
   type (gpsref_type)     :: gpsref, ggpsref

   integer :: nx, ny, nz, num_date, idate
   real    :: dx, cen_lat, cen_lon, truelat1, truelat2, stand_lon
   integer :: map_proj_wrf
   logical :: l_got_info, inside

   inside = .true.

   nml_unit      = 10
   diag_unit_in  = 50
   diag_unit_out = 20
   info_unit     = 30

   exp_num   = 0
   exp_dirs = ''
   out_dirs = ''

   if_plot_rmse  = .false.
   if_plot_bias  = .false.
   if_plot_abias = .false.

   if_plot_synop     = .false.
   if_plot_sonde_sfc = .false.
   if_plot_metar     = .false.
   if_plot_ships     = .false.
   if_plot_qscat     = .false.
   if_plot_buoy      = .false.

   if_plot_sound     = .false.
   if_plot_geoamv    = .false.
   if_plot_polaramv  = .false.
   if_plot_profiler  = .false.
   if_plot_airep     = .false.
   if_plot_pilot     = .false.

   if_plot_gpspw     = .false.
   if_plot_gpsref    = .false.
   if_plot_airsret   = .false.

   if_plot_tamdar    = .false.

   file_path_string = 'wrfvar/gts_omb_oma_01'
   wrf_file = 'foo.nc'

   istart = 1
   iend   = 10000
   jstart = 1
   jend   = 10000

   ! Read in namelist information defined in module define_cons_types

   open ( unit=nml_unit, file='namelist.plot_diag', STATUS='OLD',  &
         form='formatted' )
   read ( unit=nml_unit, nml=record1, IOSTAT=ier )
   write ( unit=*, nml = record1 )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record1'
      stop
   end if

   read ( unit=nml_unit, nml=record2, iostat=ier )
   write ( unit=*, nml = record2 )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record2'
      stop
   end if
   read ( unit=nml_unit, nml=record3, iostat=ier )
   write ( unit=*, nml = record3 )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record3'
      stop
   end if
   read ( unit=nml_unit, nml=record4, iostat=ier )
   write ( unit=*, nml = record4 )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record4'
      stop
   end if
   read ( unit=nml_unit, nml=record5, iostat=ier )
   write ( unit=*, nml = record5 )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record5'
      stop
   end if
   read ( unit=nml_unit, nml=record6, iostat=ier )
   if ( ier /= 0 ) then
      write (*,*) 'error in reading namelist record6'
      !stop
   else
      write ( unit=*, nml = record6 )
   end if
   close(nml_unit)
   call initialize_t_tab

   call get_fileinfo
   if ( l_got_info ) then
      call set_mapinfo
      istart = max(1, istart)
      iend   = min(nx, iend)
      jstart = max(1, jstart)
      jend   = min(ny, jend)
   end if

   if_plot_surface = .false.
   if (if_plot_synop .or. if_plot_metar .or. if_plot_ships .or. if_plot_buoy .or. &
        if_plot_sonde_sfc .or. if_plot_qscat   ) if_plot_surface = .true.

   if_plot_upr = .false.
   if (if_plot_sound .or. if_plot_pilot .or. if_plot_profiler   .or.    &
       if_plot_geoamv .or. if_plot_polaramv  .or. if_plot_airep .or.    &
       if_plot_airsret .or. if_plot_tamdar ) if_plot_upr= .true.

   read(start_date(1:10), fmt='(i10)')sdate
   read(end_date(1:10), fmt='(i10)')edate
   write(6,'(4a)')' Diag Starting date ', start_date, ' Ending date ', end_date
   write(6,'(a,i8,a)')' Interval between dates = ', interval, ' hours.'

   num_date = 0
   date = start_date
   do while ( date <= end_date )
      num_date = num_date + 1
      call da_advance_cymdh(date, interval, date)
   end do
   allocate(l_skip(num_date))
   l_skip(:) = .false.

   ! check for missing dates
   idate = 0  ! index of date
   date = start_date
   do while ( date <= end_date )
      idate = idate + 1
      do iexp = 1, exp_num
         filename = TRIM(exp_dirs(iexp))//'/'//date//'/'//trim(file_path_string)
         inquire ( file=trim(filename), exist=is_file)
         if ( .not. is_file ) then
            l_skip(idate) = .true.
         end if
         if ( l_skip(idate) ) exit
      end do
      call da_advance_cymdh(date, interval, date)
   end do

   !---------------------------------------------------------------------------
   ! Loop over experiments

   do iexp  =1,exp_num

      idate = 0
      date = start_date
      cdate = sdate
      call initialize_upr_type(gupr)
      call initialize_gpsref_type(ggpsref)

      do while ( cdate <= edate )
         ! Initialize various types
         call initialize_surface_type(surface)
         call initialize_upr_type(upr)
         call initialize_gpspw_type(gpspw)
         call initialize_gpsref_type(gpsref)

         idate = idate + 1

         ! construct file name
         filename = TRIM(exp_dirs(iexp))//'/'//date//'/'//trim(file_path_string)

         inquire ( file=trim(filename), exist=is_file)
         if ( l_skip(idate) .or. .not. is_file ) then
            print*, 'skipping file ', trim(filename)
            !stop
            !  Write output on outunit
            out_dir=trim(out_dirs(iexp))
            if (if_plot_surface  )  then
               call write_diag_single_level(out_dir,diag_unit_out,date,'surface_u',surface%uomb,surface%uoma)
               call write_diag_single_level(out_dir,diag_unit_out,date,'surface_v',surface%vomb,surface%voma)
               call write_diag_single_level(out_dir,diag_unit_out,date,'surface_t',surface%tomb,surface%toma)
               call write_diag_single_level(out_dir,diag_unit_out,date,'surface_p',surface%pomb,surface%poma)
               call write_diag_single_level(out_dir,diag_unit_out,date,'surface_q',surface%qomb,surface%qoma)
            end if

            if (if_plot_gpspw )  then
               call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
            end if

            if (if_plot_gpsref  )  then
             call write_diag_multi_level_h(out_dir,diag_unit_out,date,'gps_ref',gpsref%refomb,gpsref%refoma)
            end if

            if (if_plot_upr ) then
               call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_u',upr%uomb,upr%uoma)
               call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_v',upr%vomb,upr%voma)
               call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_t',upr%tomb,upr%toma)
               call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_q',upr%qomb,upr%qoma)
            end if
            !     Calculate next date:
            call da_advance_cymdh( date, interval, new_date )
            date = new_date
            read(date(1:10), fmt='(i10)')cdate
            cycle
         end if

         open (unit=diag_unit_in, file=trim(filename), form='formatted',   &
                  status='old', iostat=ier)
1           continue

         if_write = .false.
         read(diag_unit_in,'(a20,i8)', end=2000, err = 1000)obs_type,num_obs
         if (index( obs_type,'synop') > 0 ) then
            if (if_plot_synop ) if_write = .true.
            goto 10

         elseif (index( obs_type,'metar') > 0 ) then
            if (if_plot_metar ) if_write = .true.
            goto 10

         elseif (index( obs_type,'ships') > 0 )  then
            if (if_plot_ships ) if_write = .true.
            goto 10

         elseif (index( obs_type,'buoy' ) > 0 )  then
            if (if_plot_buoy ) if_write = .true.
            goto 10

         elseif (index( obs_type,'sonde_sfc') > 0 )  then
            if (if_plot_sonde_sfc ) if_write = .true.
            goto 10

         elseif (index( obs_type,'polaramv') > 0)  then
            if (if_plot_polaramv ) if_write = .true.
            goto 20

         elseif (index( obs_type,'geoamv'  ) > 0)  then
            if (if_plot_geoamv ) if_write = .true.
            goto 20

         elseif (index( obs_type,'gpspw') > 0)  then
            if ( if_plot_gpspw ) if_write = .true.
            goto 30

         elseif (index( obs_type,'sound') > 0)  then
            if (if_plot_sound ) if_write = .true.
            goto 40

         elseif (index( obs_type,'airep') > 0)  then
            if (if_plot_airep ) if_write = .true.
            goto 50

         elseif (index( obs_type,'pilot')    > 0)  then
            if (if_plot_pilot ) if_write = .true.
            goto 60

         elseif (index( obs_type,'profiler') > 0)  then
            if (if_plot_profiler ) if_write = .true.
            goto 60

         elseif (index( obs_type,'ssmir') > 0)  then
            goto 70

         elseif (index( obs_type,'ssmiT') > 0)  then
            goto 80

         elseif (index( obs_type,'satem') > 0)  then
            goto 90

         elseif (index( obs_type,'ssmt1') > 0)  then
            goto 100

         elseif (index( obs_type,'ssmt2') > 0)  then
            goto 100

         elseif (index( obs_type,'qscat') > 0)  then
            if (if_plot_qscat ) if_write = .true.
            goto 110
         elseif (index( obs_type,'gpsref' ) > 0) then
            if (if_plot_gpsref ) if_write = .true.
               goto 120
         elseif (index( obs_type,'airsr') > 0)  then
               if (if_plot_airsret ) if_write = .true.
               goto 130

         elseif (index( obs_type,'tamdar') > 0)  then
            if (if_plot_tamdar ) if_write = .true.
            goto 140

         else
         print*,' Got unknown OBS_TYPE ',obs_type(1:20),' on unit ',diag_unit_in
         stop
         end if

10       continue      !   Synop, Metar, Ships, Buoy , Sonde_sfc

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                                  kk, l, stn_id, &          ! Station
                                  lat, lon, press, &       ! Lat/lon, pressure
                                  u_obs, u_inv, u_qc, u_error, u_inc, &
                                  v_obs, v_inv, v_qc, v_error, v_inc, &
                                  t_obs, t_inv, t_qc, t_error, t_inc, &
                                  p_obs, p_inv, p_qc, p_error, p_inc, &
                                  q_obs, q_inv, q_qc, q_error, q_inc
                  if (if_write) then
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if (u_qc >=  0) call update_stats(surface%uomb, surface%uoma, u_inv, u_inc)
                        if (v_qc >=  0) call update_stats(surface%vomb, surface%voma, v_inv, v_inc)
                        if (t_qc >=  0) call update_stats(surface%tomb, surface%toma, t_inv, t_inc)
                        if (p_qc >=  0) call update_stats(surface%pomb, surface%poma, p_inv, p_inc)
                        if (q_qc >=  0) call update_stats(surface%qomb, surface%qoma, q_inv, q_inc)
                     end if
                  end if
               end do      !  loop over levels
            end do      !  loop over Obs
         end if
         goto 1

20       continue      !    Polar or Geo AMV's

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk, l, stn_id, &          ! Station
                     lat, lon, press, &        ! Lat/lon, pressure
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc

                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if( u_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
                           call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
                        endif
                        if( v_qc >= 0 .and. npr > 0 ) then
                          call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
                          call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
                        endif
                     end if
                  end if
               end do      !  loop over levels
            end do      !  loop over Obs
         end if

         goto 1

30       continue      !    Gpspw

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk, l, stn_id, &          ! Station
                     lat, lon, dummy, &       ! Lat/lon, dummy
                     tpw_obs, tpw_inv, tpw_qc, tpw_err, tpw_inc
                  if (if_write) then
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if (tpw_qc >=  0) call update_stats(gpspw%tpwomb,gpspw%tpwoma,tpw_inv,tpw_inc)
                     end if
                  end if
               end do      !  loop over levels
            end do      !  loop over Obs
         end if

         goto 1

40       continue      !    Sound

         !  [6] Transfer sound obs:

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, press, &       ! Lat/lon, dummy
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc, &
                     t_obs, t_inv, t_qc, t_error, t_inc, &
                     q_obs, q_inv, q_qc, q_error, q_inc
                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if( u_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
                           call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
                        endif
                        if( v_qc >= 0 .and. npr > 0 ) then
                          call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
                          call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
                        endif
                        if( t_qc >= 0 .and. npr > 0 )  then
                          call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
                          call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
                        endif
                        if( q_qc >= 0 .and. npr > 0 )  then
                          call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
                          call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
                        endif
                     end if
!
                  end if
                end do      !  loop over levels
             end do      !  loop over Obs
         end if
         goto 1

50       continue      !    Airep

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)') levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, press, &       ! Lat/lon, dummy
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc, &
                     t_obs, t_inv, t_qc, t_error, t_inc
                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if( u_qc >= 0 .and. npr > 0 ) then
                          call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
                          call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
                        endif
                        if( v_qc >= 0 .and. npr > 0 ) then
                          call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
                          call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
                        endif
                        if( t_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
                           call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
                        endif
                     end if

                  end if
               end do      !  loop over levels
            end do     !  loop over Obs
         end if

         goto 1

60       continue      !    Pilot & Profiler

         !  [8] Transfer pilot obs:
         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, press, &       ! Lat/lon, dummy
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc
                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if( u_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
                           call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
                        endif
                        if( v_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
                           call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
                        endif
                     end if

                  end if
               end do      !  loop over levels
            end do      !  loop over Obs
         end if
         goto 1

70       continue      !  SSMI retrievals

         if ( num_obs > 0 ) then
          do n = 1, num_obs
            read(diag_unit_in,'(i8)')levels
            do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, dummy, &       ! Lat/lon, dummy
                     spd_obs, spd_inv, spd_qc, spd_err, spd_inc
            end do
          end do
         end if

         goto 1

80       continue      !  SSMI radiance

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,*)dummy_c
               read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,7(2f17.7,i8,2f17.7))', err= 1000)&
                  k, l, stn_id, &          ! Station
                  lat, lon, dummy, &       ! Lat/lon, dummy
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy, &
                  dummy, dummy, dummy_i, dummy, dummy
            end do
         end if
         goto 1

90       continue      !  SATEM

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)') levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, dummy, &       ! Lat/lon, dummy
                     dummy,dummy, dummy_i, dummy, dummy
               end do      !  loop over levels
            end do     !  loop over Obs
         end if

         goto 1

100      continue      !  SSMT1 & 2

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)') levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, dummy, &       ! Lat/lon, dummy
                     dummy,dummy, dummy_i, dummy, dummy
               end do      !  loop over levels
            end do      !  loop over obs
         end if

         goto 1

110      continue      !  Scatrometer winds

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)') levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk, l, stn_id, &          ! Station
                     lat, lon, press, &       ! Lat/lon, dummy
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc
                  if (if_write) then
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if (u_qc >=  0) call update_stats(surface%uomb,surface%uoma,u_inv,u_inc)
                        if (v_qc >=  0) call update_stats(surface%vomb,surface%voma,v_inv,v_inc)
                     end if
                  end if
               end do      !  loop over levels
            end do      !  loop over obs
         end if
         goto 1

120      continue      !  Gpsref

   IF ( num_obs > 0 ) THEN
      DO n = 1, num_obs
         read(diag_unit_in,'(i8)') levels
         DO k = 1, levels
            read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                         kk, l, stn_id, &          ! Station
                         lat, lon, height, &       ! Lat/lon, dummy
                         ref_obs, ref_inv, ref_qc, ref_err, ref_inc
            if (if_write .and. height > 0.0) then
               call get_std_ht_level(height, nht, stdh, nstdh)
               if ( l_got_info ) call check_domain(lat, lon, inside)
               if ( inside ) then
                  if ( ref_qc >=  0) then
                     call update_stats(gpsref%refomb(nht),gpsref%refoma(nht),ref_inv,ref_inc)
                     call update_stats(ggpsref%refomb(nht),ggpsref%refoma(nht),ref_inv,ref_inc)
                  end if
               end if
            end if
         END DO      !  loop over levels
      END DO      !  loop over Obs
   ENDIF
   go to 1
!---------------------------------------------------------------------

130      continue      !  AIRSRET

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &          ! Station
                     lat, lon, press, &       ! Lat/lon, dummy
                     t_obs, t_inv, t_qc, t_error, t_inc, &
                     q_obs, q_inv, q_qc, q_error, q_inc
                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)
                     if ( l_got_info ) call check_domain(lat, lon, inside)
                     if ( inside ) then
                        if( t_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
                           call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
                        endif
                        if( q_qc >= 0 .and. npr > 0 ) then
                           call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
                           call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
                        endif
                     end if
                  end if
               end do      !  loop over levels
            end do      !  loop over obs
         end if
         goto 1

140      continue

         if ( num_obs > 0 ) then
            do n = 1, num_obs
               read(diag_unit_in,'(i8)')levels
               do k = 1, levels
                  read(diag_unit_in,'(2i8,a5,2f9.2,f17.7,5(2f17.7,i8,2f17.7))', err= 1000)&
                     kk,l, stn_id, &
                     lat, lon, press, &
                     u_obs, u_inv, u_qc, u_error, u_inc, &
                     v_obs, v_inv, v_qc, v_error, v_inc, &
                     t_obs, t_inv, t_qc, t_error, t_inc, &
                     q_obs, q_inv, q_qc, q_error, q_inc
                  if (if_write .and. press > 0 ) then
                     call get_std_pr_level(press, npr, stdp, nstd)

                   if( u_qc >=  0) then
                     call update_stats(upr%uomb(npr),upr%uoma(npr),u_inv,u_inc)
                     call update_stats(gupr%uomb(npr),gupr%uoma(npr),u_inv,u_inc)
                   endif
                   if( v_qc >=  0) then
                    call update_stats(upr%vomb(npr),upr%voma(npr),v_inv,v_inc)
                    call update_stats(gupr%vomb(npr),gupr%voma(npr),v_inv,v_inc)
                   endif
                   if( t_qc >=  0)  then
                    call update_stats(upr%tomb(npr),upr%toma(npr),t_inv,t_inc)
                    call update_stats(gupr%tomb(npr),gupr%toma(npr),t_inv,t_inc)
                   endif
                   if( q_qc >=  0)  then
                    call update_stats(upr%qomb(npr),upr%qoma(npr),q_inv,q_inc)
                    call update_stats(gupr%qomb(npr),gupr%qoma(npr),q_inv,q_inc)
                   endif

                  end if
                end do
             end do
         end if
         goto 1

         ! Now process the diagnostics
2000     continue

         close (diag_unit_in)
         !  Write output on outunit
         out_dir=trim(out_dirs(iexp))
         if (if_plot_surface  )  then
            call write_diag_single_level(out_dir,diag_unit_out,date,'surface_u',surface%uomb,surface%uoma)
            call write_diag_single_level(out_dir,diag_unit_out,date,'surface_v',surface%vomb,surface%voma)
            call write_diag_single_level(out_dir,diag_unit_out,date,'surface_t',surface%tomb,surface%toma)
            call write_diag_single_level(out_dir,diag_unit_out,date,'surface_p',surface%pomb,surface%poma)
            call write_diag_single_level(out_dir,diag_unit_out,date,'surface_q',surface%qomb,surface%qoma)
         end if

         if (if_plot_gpspw )  then
            call write_diag_single_level(out_dir,diag_unit_out,date,'gpspw_tpw',gpspw%tpwomb,gpspw%tpwoma)
         end if

         if (if_plot_gpsref  )  then
          call write_diag_multi_level_h(out_dir,diag_unit_out,date,'gps_ref',gpsref%refomb,gpsref%refoma)
!rizvi          call write_diag_single_level(out_dir,diag_unit_out,date,'avgh_ref',avgh%refomb,avgh%refoma)
         end if

         if (if_plot_upr ) then
            call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_u',upr%uomb,upr%uoma)
            call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_v',upr%vomb,upr%voma)
            call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_t',upr%tomb,upr%toma)
            call write_diag_multi_level(out_dir,diag_unit_out,date,'upr_q',upr%qomb,upr%qoma)
         end if
         !     Calculate next date:
         call da_advance_cymdh( date, interval, new_date )
         date = new_date
         read(date(1:10), fmt='(i10)')cdate
      end do     ! End loop over date
       if( if_plot_upr ) then
        call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_u',gupr%uomb,gupr%uoma)
        call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_v',gupr%vomb,gupr%voma)
        call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_t',gupr%tomb,gupr%toma)
        call write_diag_multi_level(out_dir,diag_unit_out,date,'gupr_q',gupr%qomb,gupr%qoma)
       endif
       if (if_plot_gpsref  )  then
        call write_diag_multi_level_h(out_dir,diag_unit_out,date,'ggps_ref',ggpsref%refomb,ggpsref%refoma)
       endif

   end do   ! Loop over experiments
   stop

1000 print*,' Error while reading unit ',diag_unit_in,' for experiment ',exp_dirs(iexp)
   stop

contains

subroutine get_std_pr_level(prs, npr, stdp, nstd)

   implicit none

   integer, intent(in )      :: nstd
   real,    intent(in)       :: stdp(nstd)
   integer, intent(out)      :: npr
   real,    intent(in)       :: prs

   real             :: pr
   integer               :: k

   npr = num_miss  ! initialize as a missing value
   pr = prs/100.0
   if        ( pr >= stdp(1)    ) then
       npr = 1
       return
   else if ( pr < stdp(nstd-1) .and. pr >= stdp(nstd) ) then
      npr = nstd
      return
   else
      do k = 2,nstd - 1
         if (pr   >= stdp(k) ) then
             npr = k
            return
         end if
      end do
   end if

end subroutine get_std_pr_level

subroutine get_std_ht_level(height, nht, stdh, nstdh)

   implicit none

   integer, intent(in )      :: nstdh
   real,    intent(in)       :: stdh(nstdh)
   integer, intent(out)      :: nht
   real,    intent(in)       :: height

   real    :: ht
   integer :: k

   ht = height*0.001  ! m to km
   if ( ht <= stdh(1)    ) then
      nht = 1
      return
   else if ( ht > stdh(nstdh-1) ) then
      nht = nstdh
      return
   else
      do k = 2,nstdh - 1
         if ( ht <= stdh(k) ) then
            nht = k
            return
         end if
      end do
   end if

end subroutine get_std_ht_level

subroutine update_stats(stats_omb, stats_oma, omb, oma)

   implicit none
   type(stats_value),   intent(inout)   :: stats_omb, stats_oma
   real, intent (in)                    :: omb, oma

   real      :: x1, x2

   stats_omb%num  = stats_omb%num + 1
   stats_oma%num  = stats_omb%num

   x1 = 1.0/ stats_omb%num
   x2 =  (stats_omb%num-1)*x1

   stats_omb%bias  = x2*stats_omb%bias  + omb  * x1
   stats_oma%bias  = x2*stats_oma%bias  + oma  * x1

   stats_omb%abias = x2*stats_omb%abias + abs(omb) * x1
   stats_oma%abias = x2*stats_oma%abias + abs(oma) * x1

   stats_omb%rmse  = x2*stats_omb%rmse  + omb*omb  * x1
   stats_oma%rmse  = x2*stats_oma%rmse  + oma*oma  * x1

end subroutine update_stats

subroutine write_diag_single_level(out_dir,ounit,ldate,obs_type,omb,oma)

   implicit none

   integer, intent(in)            :: ounit
   character*512,intent(in)       :: out_dir
   character*10,intent(in)        :: ldate
   character*(*),intent(in)       :: obs_type
   type (stats_value),intent(in)  :: omb
   type (stats_value),intent(in)  :: oma

   character*512                  :: filename
   integer                        :: ounit1, ounit2
   real                           :: sigt,bar

   ounit1 = ounit
   ounit2 = ounit + 1

   filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
   open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
   filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
   open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')
   if ( omb%num <= 1 ) then
     sigt=1.       ; bar =  rmiss
     write(ounit1,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
     write(ounit2,'(1x,a10,1x,6(f6.2,1x))') ldate,rmiss, rmiss, rmiss, rmiss,bar,sigt
   else
      ! write(ounit1,'(5x,a10,4(2x,a9))') trim(obs_type),' Number','BIAS','ABIAS','RMSE '
     if (index(obs_type,'_q') > 0 ) then
     call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
     bar=bar*1000.0
     write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,omb%num, 1000.0*omb%bias, 1000.0*omb%abias, 1000.0*sqrt(omb%rmse),bar,sigt
     call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
     bar=bar*1000.0
     write(ounit2,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,oma%num, 1000.0*oma%bias, 1000.0*oma%abias, 1000.0*sqrt(oma%rmse),bar,sigt
     else if( index(obs_type,'_p') > 0 ) then
     call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
     bar=bar/100.0
     write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))')ldate,omb%num, omb%bias/100.0, omb%abias/100.0, sqrt(omb%rmse)/100.0,bar,sigt
     call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
     bar=bar/100.0
     write(ounit2,'(1x,a10,1x,i9,5(1x,f6.2))') ldate,oma%num, oma%bias/100.0, oma%abias/100.0, sqrt(oma%rmse)/100.0,bar,sigt
     else
     call sig_test(omb%num, omb%bias, omb%rmse, sigt,bar)
     write(ounit1,'(1x,a10,1x,i9,1x,5(f6.2,1x))') ldate,omb%num, omb%bias, omb%abias, sqrt(omb%rmse),bar,sigt
     call sig_test(oma%num, oma%bias, oma%rmse, sigt,bar)
     write(ounit2,'(1x,a10,1x,i9,5(1x,f6.2))') ldate,oma%num, oma%bias, oma%abias, sqrt(oma%rmse),bar,sigt
     endif
!
   end if
   close(ounit1)
   close(ounit2)

end subroutine write_diag_single_level

subroutine write_diag_multi_level(out_dir,ounit,ldate,obs_type,omb,oma)

   implicit none

   integer, intent(in)            :: ounit
   character*512,intent(in)       :: out_dir
   character*10,intent(in)        :: ldate
   character*(*),intent(in)       :: obs_type
   type (stats_value),intent(in)  :: omb(nstd)
   type (stats_value),intent(in)  :: oma(nstd)

   character*512                  :: filename
   integer                        :: k
   real                           :: xnum(nstd)
   integer                        :: num(nstd)
   real, dimension(nstd)          :: rmse, bias, abias,sigt,bar
   integer                        :: ounit1, ounit2

   ounit1 = ounit
   ounit2 = ounit + 1

   filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
   open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
   filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
   open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')

   do k = 1, nstd
      num(k) = omb(k)%num
      if (num(k) <= 1 ) then
         xnum(k)  = rmiss
         rmse(k)  = rmiss
         bias(k)  = rmiss
         abias(k) = rmiss
         bar(k)   = rmiss
         sigt(k)  = 1.0
      else
         if (index(obs_type,'_q') > 0 ) then
            rmse(k) = sqrt(omb(k)%rmse) * 1000
            bias(k) = omb(k)%bias * 1000
            abias(k) = omb(k)%abias * 1000
            call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
            bar(k) = bar(k)*1000.
         else
            rmse(k) = sqrt(omb(k)%rmse)
            bias(k) = omb(k)%bias
            abias(k) = omb(k)%abias
            call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
         end if
      xnum(k) = num(k)
      end if
   end do

    write(ounit1,'(1x,a10,1x,16(6(1x,f12.2)))')ldate, (xnum(k), bias(k), abias(k),&
         rmse(k),bar(k),sigt(k),k=1,nstd)

   do k = 1, nstd
      num(k) = oma(k)%num
      if( num(k) <= 1 ) then
         xnum(k)  = rmiss
         rmse(k)  = rmiss
         bias(k)  = rmiss
         abias(k) = rmiss
      else
         if (index(obs_type,'_q') > 0 ) then
            rmse(k) = sqrt(oma(k)%rmse) * 1000
            bias(k) = oma(k)%bias * 1000
            abias(k) = oma(k)%abias * 1000
            call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
            bar(k) = bar(k)*1000.
         else
            rmse(k) = sqrt(oma(k)%rmse)
            bias(k) = oma(k)%bias
            abias(k) = oma(k)%abias
            call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
         end if
       xnum(k) = num(k)
      end if
   end do
   write(ounit2,'(1x,a10,1x,16(6(1x,f12.2)))')ldate, (xnum(k), bias(k), abias(k),&
             rmse(k),bar(k),sigt(k),k=1,nstd)

   close(ounit1)
   close(ounit2)

end subroutine write_diag_multi_level
     subroutine write_diag_multi_level_h(out_dir,ounit,date,obs_type,omb,oma)
     implicit none
     integer, intent(in)            :: ounit
     character*512,intent(in)       :: out_dir
     character*10,intent(in)        :: date
     character*(*),intent(in)       :: obs_type
     type (stats_value),intent(in)  :: omb(nstdh)
     type (stats_value),intent(in)  :: oma(nstdh)
!
     character*512                  :: filename
     integer                        :: k
     real                           :: xnum(nstdh)
     integer                        :: num(nstdh)
     real, dimension(nstdh)         :: rmse, bias, abias, sigt, bar

     integer                        :: ounit1, ounit2
!
     ounit1 = ounit
     ounit2 = ounit + 1
!

     filename = trim(out_dir)//'/'//trim(obs_type)//'_omb.diag'
     open (ounit1, file = trim(filename), form='formatted',status='unknown',position='append')
     filename = trim(out_dir)//'/'//trim(obs_type)//'_oma.diag'
     open (ounit2, file = trim(filename), form='formatted',status='unknown',position='append')

     do k = 1, nstdh
     num(k) = omb(k)%num
     if( num(k) <= 1 ) then
     xnum(k)  = rmiss
     rmse(k)  = rmiss
     bias(k)  = rmiss
     abias(k) = rmiss
     bar(k)   = rmiss
     sigt(k)  = 1.0
     else
        if( index(obs_type,'_q') > 0 ) then

         rmse(k) = sqrt(omb(k)%rmse) * 1000
         bias(k) = omb(k)%bias * 1000
         abias(k) = omb(k)%abias * 1000
         call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
         bar(k) = bar(k)*1000.
       else
        rmse(k) = sqrt(omb(k)%rmse)
        bias(k) = omb(k)%bias
        abias(k) = omb(k)%abias
        call sig_test(num(k), omb(k)%bias, omb(k)%rmse, sigt(k),bar(k))
       endif
      xnum(k) = num(k)
     endif
     enddo
     write(ounit1,'(1x,a10,1x,150(6(1x,f12.2)))')date, (xnum(k), bias(k), abias(k), &
           rmse(k),bar(k),sigt(k), k=1,nstdh)

     do k = 1, nstdh
     num(k) = oma(k)%num
     if( num(k) <= 1 ) then
     xnum(k)  = rmiss
     rmse(k)  = rmiss
     bias(k)  = rmiss
     abias(k) = rmiss
     bar(k)   = rmiss
     sigt(k)  = 1.0
     else
        if( index(obs_type,'_q') > 0 ) then

         rmse(k) = sqrt(oma(k)%rmse) * 1000
         bias(k) = oma(k)%bias * 1000
         abias(k) = oma(k)%abias * 1000
         call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
         bar(k) = bar(k)*1000.
       else
        rmse(k) = sqrt(oma(k)%rmse)
        bias(k) = oma(k)%bias
        abias(k) = oma(k)%abias
        call sig_test(num(k), oma(k)%bias, oma(k)%rmse, sigt(k),bar(k))
       endif
      xnum(k) = num(k)
     endif
     enddo
     write(ounit2,'(1x,a10,1x,150(6(1x,f12.2)))')date, (xnum(k), bias(k), abias(k), &
           rmse(k),bar(k),sigt(k), k=1,nstdh)

!
     close(ounit1)
     close(ounit2)
!
     end subroutine write_diag_multi_level_h
!
     subroutine sig_test(num, bias, rmse, sigt,bar)
     implicit none
     integer, intent(in)      :: num
     real,    intent(in)      :: bias, rmse
     real,    intent(out)     :: sigt, bar

     real                     :: t_val, sd, tmp
!
     sigt=0.
     tmp = num/real(num-1)
!     sd = sqrt ( tmp*rmse - bias*bias/real((n-1)) )
     sd = sqrt ( tmp*( rmse - bias*bias )  )
     do k=2,34
       if (real(num-1) < alpha(k,1)) exit
     end do

     t_val = bias*sqrt( real(num) ) /sd
     bar = alpha(k-1,2) * sd /sqrt( real(num) )

     if (abs(t_val) >= alpha(k-1,2)) sigt=1.

    end subroutine sig_test

   subroutine get_fileinfo
      implicit none
#include "netcdf.inc"
      integer :: iost(12)
      integer :: ncid
      l_got_info = .false.
      iost(1) = nf_open(trim(wrf_file), NF_NOWRITE, ncid)
      if ( iost(1) /= NF_NOERR ) then
         print*, 'INFO: wrf_file: '//trim(wrf_file)//' does not exist for retrieving mapping info'
         return
      else
         print*, 'Retrieving mapping info from wrf_file: ',trim(wrf_file)
      end if
      iost(2)  = nf_get_att_int(ncid, NF_GLOBAL, 'WEST-EAST_GRID_DIMENSION', nx)
      iost(3)  = nf_get_att_int(ncid, NF_GLOBAL, 'SOUTH-NORTH_GRID_DIMENSION', ny)
      iost(4)  = nf_get_att_int(ncid, NF_GLOBAL, 'BOTTOM-TOP_GRID_DIMENSION', nz)
      iost(5)  = nf_get_att_double(ncid, NF_GLOBAL, 'DX', dx)
      iost(6)  = nf_get_att_double(ncid, NF_GLOBAL, 'CEN_LAT', cen_lat)
      iost(7)  = nf_get_att_double(ncid, NF_GLOBAL, 'CEN_LON', cen_lon)
      iost(8)  = nf_get_att_double(ncid, NF_GLOBAL, 'TRUELAT1', truelat1)
      iost(9)  = nf_get_att_double(ncid, NF_GLOBAL, 'TRUELAT2', truelat2)
      iost(10) = nf_get_att_double(ncid, NF_GLOBAL, 'STAND_LON', stand_lon)
      iost(11) = nf_get_att_int(ncid, NF_GLOBAL, 'MAP_PROJ', map_proj_wrf)
      iost(12) = nf_close(ncid)
      if ( .not. any(iost/=NF_NOERR) ) then
         l_got_info = .true.
      end if
      print*, 'nx, ny, nz, dx, map_proj, cen_lat, cen_lon, truelat1, truelat2, stand_lon = ', &
         nx, ny, nz, dx, map_proj_wrf, cen_lat, cen_lon, truelat1, truelat2, stand_lon
   end subroutine get_fileinfo

   subroutine set_mapinfo
      implicit none
      integer :: map_proj_util
      real    :: xref, yref
      real :: start_x, start_y
      xref = nx/2.0
      yref = ny/2.0
      if ( map_proj_wrf == 0 .or. map_proj_wrf == 6 ) then
         map_proj_util = proj_latlon
      else if ( map_proj_wrf == 3 ) then
         map_proj_util = proj_merc
      else if ( map_proj_wrf == 1 ) then
         map_proj_util = proj_lc
      else if ( map_proj_wrf == 2 ) then
         map_proj_util = proj_ps
      end if
      call da_map_set(map_proj_util, cen_lat,cen_lon, xref, yref, dx, &
         stand_lon, truelat1, truelat2, truelat1, stand_lon, map_info)
      !call da_llxy_wrf(map_info, cen_lat, cen_lon, start_x, start_y)
   end subroutine set_mapinfo

   subroutine check_domain(lat, lon, inside)
      implicit none
      real,    intent(in)  :: lat, lon
      logical, intent(out) :: inside
      real                 :: xx, yy
      call da_llxy_wrf(map_info, lat, lon, xx, yy)
      inside = .false.
      if ( xx >= istart .and. xx <= iend .and.  &
           yy >= jstart .and. yy <= jend ) then
         inside = .true.
      end if
   end subroutine check_domain

end program da_verif_obs

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb