text:
!-----------------------------------------------------------------------
      subroutine zndf(y,n,norig,ndf)
!-----------------------------------------------------------------------
!
!      subroutine to calculate number of degrees of freedom (ndf) of a time
!      series y, with length n, and physical dimension np. If y has been
!      interpolated then norig is the original number of points in the shear
!      wave window
!
!      calculated according to appendix in Silver and Chan 1991
!
!      For a gaussian white noise ndf should be equal to n.
!      The true ndf should be less than n
!
!      variables
!    in:
!      y(np)                    real      time series
!      n                        int      number of points
!      norig                    int      original (uninterpolated) number of points
!                               in the shear wave window
!    out:
!      ndf                        int      number of degrees of freedom
!    local:
!      np                        int      array dimension
!                                    (read from SIZE_np.h at compile time)
!      n2                        int      next power of 2 up from n
!      yint(np)                  real      y interpolated onto n2 points
!      yfft(np*2)                   real      fft of interpolated (to n2) y in NR format
!      yfft_amp(np)             real      amplitude spectum
!      nf                        int      number of freq points for n data points
!      yfft_amp                  real      interpolated amp spec on nf points
!      f2,f4                        real      f2 and f4 from Silver and Chan 1991
!                                    used to calc ndf
!
!-----------------------------------------------------------------------
!      N. Teanby      5-7-02   Original code, based on ndf_fun2.f
!                              and ndf_spect.f in DTM
!      J. Wookey      7-4-14   Modified ndf calculation following
!                              Walsh, Arnold and Savage (JGR, 2013)
!-----------------------------------------------------------------------
      use array_sizes ! use the array_sizes modules
!-----------------------------------------------------------------------
      implicit none
      integer ndf
      integer n,n2,nf,norig
      real y(np),yint(np),yfft(np*2),yfft_amp(np)
      real sce2,sce4
      integer scndf
      real e2,e4
      integer i

!  -- calc the fft of y --
!  ** calc n2, the next power of 2 from n (or n if n=power of 2) **
      n2 = 2**(int(log(real(n)-0.1)/log(2.0)) + 1)
!  ** interpolate to have n2 points **
      call zlinint(y,n,np,n2,yint)
!  ** perform fft (1 means do forward fft)**
      call zfft(yint,n2,np,1,yfft)
!  ** determine number of frequency points for original time series.
      nf=(norig+2)/2
!  ** calculate amplitude spectrum **
      do i=1,nf
         yfft_amp(i) = sqrt(yfft(2*i-1)**2 + yfft(2*i)**2)
      enddo

      e2=0.
      e4=0.
      do i=1,nf
   !!          yfft_amp(i)**2 is |f_i|^2|g_i|^2
   !!          note that have used i as the index rather than n
         e2   = e2 + yfft_amp(i)**2
   !!          yfft_amp(i)**4 is |f_i|^4|g_i|^4
   !!          note that have used i as the index rather than n
         e4   = e4 + yfft_amp(i)**4

         if(i.eq.1.or.i.eq.nf)then
   !!            coefficients from SC equation A3
   !!            note adding 1/2 is the same as adding 1 then subtracting 1/2
            sce2 = e2 - 0.5*yfft_amp(i)**2
   !!            coefficients from SC original code (see Walsh eq (33))
   !!            note adding 1/2 is the same as adding 1 then subtracting 1/2
            sce4 = e4 - 0.5*yfft_amp(i)**4
         endif
      enddo
      e2 = sce2

!!    Walsh et al eq 33 shows that the coefficients are
!!    1/3 for n=0,N and 4/3 otherwise
!!    e4 has coefficients 1 for all n so we multiply
!!    e4 by 4/3. The end points (n=0,N) now have coefficients
!!    4/3 so we must subtract off 1 to get 1/3
      e4 = 4.*e4/3. - yfft_amp(1)**4 - yfft_amp(nf)**4
      ndf = nint( 2.0 * (2.0*e2**2/e4 - 1.) )

      ! original version
      scndf = nint( 2.0 * (2.0*sce2**2/sce4 - 1.) )

!      print *, "e2,e4,sce2,sce4,e2,e4,ndf,scndf are "
!      print *, e2,e4,sce2,sce4,e2,e4,ndf,scndf

      return
      end

! MFAST VERSION
! c-----------------------------------------------------------------------
!    subroutine zndf(y,n,norig,ndf)
! c-----------------------------------------------------------------------
! c
! c   subroutine to calculate number of degrees of freedom (ndf) of a time
! c   series y, with length n, and physical dimension np. If y has been
! c   interpolated then norig is the original number of points in the shear
! c   wave window
! c
! c   calculated according to appendix in Silver and Chan 1991
! c    Modifed by Walsh, Arnold and Savage 2013. (Silver and Chan Revisited,
! c     submitted to J. Geophysical Research June 18, 2013.
! c
! c   For a gaussian white noise ndf should be equal to n.
! c   The true ndf should be less than n
! c
! c   variables
! c    in:
! c   y(np)            real   time series
! c   n            int   number of points
! c   norig            int   original (uninterpolated) number of points
! c                  in the shear wave window
! c    out:
! c   ndf            int   number of degrees of freedom
! c    local:
! c   np            int   array dimension
! c                  (read from SIZE_np.h at compile time)
! c   n2            int   next power of 2 up from n
! c   yint(np)         real   y interpolated onto n2 points
! c   yfft(np*2)          real   fft of interpolated (to n2) y in NR format
! c   yfft_amp(np)       real   amplitude spectum
! c   nf            int   number of freq points for n data points
! c   yfft_amp         real   interpolated amp spec on nf points
! c   sce2,sce4            real   f2 and f4 from Silver and Chan 1991
! c       scndf,ndf            int   number off degrees of freedom from Silver and Chan--i.e., previous codes
! c                  used to calc ndf
! c       e2, e4            real   e2 and e4 from Walsh et al. 2013
!
! c
! c-----------------------------------------------------------------------
! c   N. Teanby   5-7-02   Original code, based on ndf_fun2.f
! c               and ndf_spect.f in DTM
! c
! c   M. Savage 15 June 2013  Fixing to use modified f2 and f4 from Walsh et al.
! c
! c       Walsh, E., R. Arnold and M. K. Savage (2013) "Silver and Chan
! c      Revisited), submitted to JGR June 2013.
! c-----------------------------------------------------------------------
!
!    implicit none
!    integer ndf
!    integer n,np,n2,nf,norig
!    include "SIZE_np.h"
!    real y(np),yint(np),yfft(np*2),yfft_amp(np)
! !   real f2,f4
!    real sce2,sce4
!         integer scndf
!    real e2,e4
!    integer i
!
! c  -- calc the fft of y --
! c  ** calc n2, the next power of 2 from n (or n if n=power of 2) **
!    n2 = 2**(int(log(real(n)-0.1)/log(2.0)) + 1)
! c  ** interpolate to have n2 points **
!    call zlinint(y,n,np,n2,yint)
! c  ** perform fft (1 means do forward fft)**
!    call zfft(yint,n2,np,1,yfft)
! c  ** determine number of frequency points for original time series.
!    nf=(norig+2)/2
! c  ** calculate amplitude spectrum **
!    do 2 i=1,nf
!       yfft_amp(i) = sqrt(yfft(2*i-1)**2 + yfft(2*i)**2)
! 2   continue
!
! !! below is the original code based on the SC codes
! !! they refer to f2 and f4 but these are actually
! !! e2 and e4
! !! SC note that the noise is a convolution E \approx \sum_{n=0}^N a_n x f_n*g_n
! !! where a is 1/2 for n=0,N and 1 otherwise(SC equation A3)
! !! we cannot access f but we can estimate it by calculating E
! !! from the data
! !!c  -- calc f2 and f4 from Silver and Chan 1991 --
! !! which are actually e2 and e4
! !!   f2=0.
! !!   f4=0.
! !!   do 7 i=1,nf
! !!        yfft_amp(i)**2 is |f_i|^2|g_i|^2
! !!        note that have used i as the index rather than n
! !!     f2   = f2 + yfft_amp(i)**2
! !!        yfft_amp(i)**4 is |f_i|^4|g_i|^4
! !!        note that have used i as the index rather than n
! !!     f4   = f4 + yfft_amp(i)**4
! !!cc     for zero frequency and for Nyquist:
! !!c     if(i.eq.1.or.i.eq.nf)then
! !!          coefficients from SC equation A3
! !!          note adding 1/2 is the same as adding 1 then subtracting 1/2
! !!c       f2 = f2 - 0.5*yfft_amp(i)**2
! !!          coefficients from SC original code (see Walsh eq (33))
! !!          note adding 1/2 is the same as adding 1 then subtracting 1/2
! !!c       f4 = f4 - 0.5*yfft_amp(i)**4
! !!c     endif
!
!
! !!  testing SC coefficients to compare  (equation 33 Walsh et al.)
! !!  relabel coefficients to make it less confusing
!
!    e2=0.
!    e4=0.
!    do 7 i=1,nf
! !!        yfft_amp(i)**2 is |f_i|^2|g_i|^2
! !!        note that have used i as the index rather than n
!      e2   = e2 + yfft_amp(i)**2
! !!        yfft_amp(i)**4 is |f_i|^4|g_i|^4
! !!        note that have used i as the index rather than n
!      e4   = e4 + yfft_amp(i)**4
!
!      if(i.eq.1.or.i.eq.nf)then
! !!          coefficients from SC equation A3
! !!          note adding 1/2 is the same as adding 1 then subtracting 1/2
!        sce2 = e2 - 0.5*yfft_amp(i)**2
! !!          coefficients from SC original code (see Walsh eq (33))
! !!          note adding 1/2 is the same as adding 1 then subtracting 1/2
!        sce4 = e4 - 0.5*yfft_amp(i)**4
!      endif
! 7    continue
! c
! !!  put in Walsh' coefficients
! !!  Replace sce2 and sce4 with e2 and e4 from Walsh et al. (eqn 25 and 27)
! !!  Below eqn 26 it notes that the bn coefficients in eq 25 are equal to the an
! !!  coefficients
! !!  Note that the first equation in section 2.5 (which happens to be un-numbered)
! !!  mentions that the an coefficients are 1/2 for n=0,N and 1 otherwise
! !!  so the sce2 and e2 are the same
! !!  Walsh et al eq 25
!         e2 = sce2
! !!  below is the main point of different between the two degrees of freedom
! !!  compare Walsh et al (eq 27,28,33)
! !!  note that eq 27 = e4 and eq 33=sce4
! !!  eq 28 which is what is in the SC article has only ever been implemented
! !!  in the original SC codes and later scrapped because eq 33 gave superior
! !!  performance (see Walsh et al section 2.5 and Figure 3)
! !!  you should find that ndf \approx scndf * 3/4 roughly
!
! !!  Walsh et al eq 33 shows that the coefficients are
! !!  1/3 for n=0,N and 4/3 otherwise
! !!  e4 has coefficients 1 for all n so we multiply
! !!  e4 by 4/3. The end points (n=0,N) now have coefficients
! !!  4/3 so we must subtract off 1 to get 1/3
!         e4 = 4.*e4/3. - yfft_amp(1)**4 - yfft_amp(nf)**4
!
! !! estimate the degrees of freedom using (Walsh et al 31 aka SC A12)
! !! based on our defintion of e4 (Walsh et al eq 27)
!    ndf = nint( 2.0 * (2.0*e2**2/e4 - 1.) )
! !! estimate the degrees of freedom using (Walsh et al 31 aka SC A12)
! !! based on our the original SC code definition of e4 (Walsh et al eq 33)
!    scndf = nint( 2.0 * (2.0*sce2**2/sce4 - 1.) )
! c  MKS 15 June 2013 replacing with Walsh method (Equation 31)
!         print *, "e2,e4,sce2,sce4,e2,e4,ndf,scndf are "
!         print *, e2,e4,sce2,sce4,e2,e4,ndf,scndf
!
!         open (unit=59,file="df.txt",action="write"
!      x  ,form="formatted")
! c        open (unit=59,file="df.txt",action="write",status="new"
! c     x  ,form="formatted")
!         write (59,*) e2,e4,sce2,sce4,e2,e4,ndf,scndf
!         close (59)
! c
!    return
!    end

max_stars_repo_name: jwookey/sheba

max_issues_repo_name: jwookey/sheba

max_forks_repo_name: jwookey/sheba


text:
##||Address()||

Rebecca Wendlandt Custom Clothing Designs is a custom fashion designer located in Davis, California.  She creates oneofakind garments including custom wedding dresses, evening wear, theatrical and dance costumes, and http://www.rebeccawendlandt.com/wearableart.htm Wearable Art. She also provides professional alterations and teaches private sewing lessons ranging from beginning to advanced sewing.

Services:
  Individually made dresses including: http://www.rebeccawendlandt.com/weddingdresses.htm custom wedding dresses, http://www.rebeccawendlandt.com/eveningwear.htm evening wear, and prom dresses.
  http://www.rebeccawendlandt.com/theatricalcostumedesigns.htm Theatrical Costume design for plays, ballet, and dance performances
  http://www.rebeccawendlandt.com/alterations.htm Clothing alterations are available, but her schedule typically books out months in advance.
  Oneonone http://www.rebeccawendlandt.com/sewinglessons.htm sewing lessons for all skill levels including hand sewing, machine sewing, draping and pattern making.

20101213 13:14:35 nbsp I commissioned Rebecca to custom make my wedding garment. Rebecca made the process collaborative and enjoyable. The product was STELLAR. I wanted a garment I could wear again and again, not just on the special day, and she succeeded with flying colors. She was professional, prompt, and affordable. She had a great eye for detail. Id recommend her to anyone anywhere! Thank you Rebecca!  Users/themichelle

20120819 13:24:05 nbsp Rebecca is a great teacher! I needed help with my tailoring and sewing skills and she was very skilled! Willingess to work with your schedule is a plus!  Users/hokusman

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c    This code is not protected by the DART copyright agreement.
c    DART $Id$

c    Yet another revision of the BUFR preparation program.  It has options
c    to output specific humidity, relative humidity, or dewpoint obs, under
c    namelist control.  New code added by Ryan Torn.  14 Jan 2011
c    (Requires the matching update to the create_real_obs program.)
c
c    This is the 5th version of the BUFR preparation program for the
c    newly revised DART version. Ps and moisture obs are outputed as well.
c
c    12/2009 added additional time window flavors, specified all of the variable
c    types.   (this is planned to be removed and observation windowing done
c    as a post process.  14jan2011)
c
c    A bug with the radiosonde T events selection j2 is fixed at 04/06/2005.
c    and surface data is included in this version.
c    For radiosonde T: T of Pc=1 and 6 is dry T; T of pc=8 is virtual T.
c    For Aircraft T: No moistuerobs and all T are dry T
c    For Surface obs, the T of pc=1 is already virtual T.
c
c    the READPB() routine currently has a select for obs types based on
c    name; this should be removed and completely under namelist control.
c    the previously encountered problem was fixed and not related to this.
c    See the prepdecode/docs directory for the key to all the bufr codes.
c
      REAL*8     R8BFMS
      PARAMETER ( R8BFMS = 10.0E10 )
C                                      "Missing" value for BUFR data
      PARAMETER ( MXR8PM = 10 )
C                                      Maximum number of BUFR parameters
      PARAMETER ( MXR8LV = 400 )
C                                      Maximum number of BUFR levels
      PARAMETER ( MXR8VN = 10 )
C                                      Maximum number of BUFR event sequences
      PARAMETER ( MXR8VT = 6 )
C                                      Maximum number of BUFR variable types
      PARAMETER ( MXSTRL = 80 )
C                                      Maximum size of a string
C
      REAL*8 hdr (MXR8PM), evns(MXR8PM, MXR8LV, MXR8VN, MXR8VT)
C
      COMMON /PREPBC/ hdr, evns, nlev

      PARAMETER (STRLN = 180)

      INTEGER, PARAMETER :: max_otype=300, max_qctype=15
      REAL, PARAMETER    :: MISSING = -9999.
      INTEGER, PARAMETER :: I_MISSING = -9999
      CHARACTER  outstg*(200), subset*8, inf*19 , outf*20
      CHARACTER  var(MXR8VT)/'P','Q','T','Z','U','V'/
C
      character (len=9),  parameter :: infile  = 'prepqm.in'
      character (len=10), parameter :: outfile = 'prepqm.out'

      PARAMETER  ( NFILO = 15 )
      INTEGER  iunso ( NFILO )
     +          /   51,   52,   53,   54,   55,
     +              56,   57,   58,   59,   60,
     +              61,   62,   63,   64,   65  /
      CHARACTER*6     filo ( NFILO )
C     original list:
     +          / 'ADPUPA', 'AIRCAR', 'AIRCFT', 'SATWND', 'PROFLR',
     +            'VADWND', 'SATBOG', 'SATEMP', 'ADPSFC', 'SFCSHP',
     +            'SFCBOG', 'SPSSMI', 'SYNDAT', 'ERS1DA', 'GOESND'  /
C    updated list if you want to handle quikscat winds.
C     +          / 'ADPUPA', 'AIRCAR', 'AIRCFT', 'SATWND', 'PROFLR',
C     +            'VADWND', 'SATEMP', 'ADPSFC', 'SFCSHP', 'SFCBOG',
C     +            'SPSSMI', 'SYNDAT', 'ERS1DA', 'GOESND', 'QKSWND'/

      dimension tdata(8), udata(8), vdata(8), qdata(8), pdata(8)
      dimension zdata(8)
      integer :: wtype, ptype, qtype, ttype, ztype
c    The pc values are the 'program codes' that tell you what processing
c    was done on this observation.  As of now, these are unused, but could
c    be used for selection or diagnosis.
      integer :: pc_t, pc_q, pc_u, pc_v, pc_p, pc_z
      integer :: tqm, pqm, qqm, uqm, vqm, zqm
      logical :: found, uotype, uqcflag, use_this_data_real,
     +           use_this_data_int, processed
      logical :: debug = .false.

c    Namelist Parameters
c----------------------------------------------------------------------

      integer :: qctype_use(max_otype)  ! qc values to accept (default all)
      real :: otype_use(max_otype),     ! report types to use (default all)
     +        obs_window       = 1.5,   ! observation time window (+/-hours)
     +        obs_window_upa   = 1.5,   ! sonde time window (+/-hours)
     +        obs_window_air   = 1.5,   ! aircraft obs time window (+/-hours)
     +        obs_window_sfc   = 0.8,   ! surface obs time window (+/-hours)
     +        obs_window_cw    = 1.5,   ! cloud wind obs time window (+/-hours)
     +        land_temp_error  = 2.5,   ! assumed err surface temp. obs (K)
     +        land_wind_error  = 3.5,   ! assumed err surface wind obs (m/s)
     +        land_moist_error = 0.2    ! assumed err surface moist. obs (%)

      namelist /prep_bufr_nml/ obs_window,
     +                         obs_window_upa,
     +                         obs_window_air,
     +                         obs_window_sfc,
     +                         obs_window_cw,
     +                         otype_use,
     +                         qctype_use,
     +                         land_temp_error,
     +                         land_wind_error,
     +                         land_moist_error

c    constants needed to compute saturation water vapor for qobs
c----------------------------------------------------------------------
      real*8  l0

      eps= .622
      es0= 6.11          ! in mb, qsat pressure at 273.16K
      cl = 4187.
      cpv= 1876.5
      rv = 461.5         ! m2/s2/K
      l0 = 2.501e6
      t0 = 273.16
      fact1 = (cpv - cl) / rv
      fact2 = (l0 + (cl - cpv) * t0) / rv
      fact3 = 1. / t0
      omeps = 1.-eps

c    read the prep_bufr_nml namelist from input.nml
c----------------------------------------------------------------------

      otype_use(:) = MISSING;  qctype_use(:) = I_MISSING
      inml_unit = 15
      open(inml_unit,file='input.nml', status='old',
     &        access='sequential', form='formatted')
      read(inml_unit, nml = prep_bufr_nml)
      close(inml_unit)

      inum_otype=0
      do i = 1, max_otype
        if ( otype_use(i) .eq. MISSING ) exit
        inum_otype = inum_otype + 1
      enddo

      inum_qctype=0
      do i = 1, max_qctype
        if ( qctype_use(i) .eq. MISSING ) exit
        inum_qctype = inum_qctype + 1
      enddo

c     overwrite observation windows if obs_window is defined
c----------------------------------------------------------------------

      if ( obs_window > 0.0 ) then
        print*, "Using the same observation window for all data: ",
     &     obs_window
        obs_window_upa = obs_window
        obs_window_air = obs_window
        obs_window_sfc = obs_window
        obs_window_cw  = obs_window
      endif

c    open the input bufr file and the output text file
c----------------------------------------------------------------------

      iprof = 0
      lunobs=210
      open(unit=lunobs, file = outfile, form='formatted')

      ibufr_unit = 11
      OPEN  ( UNIT = ibufr_unit, FILE = infile, FORM = 'UNFORMATTED' )
      CALL OPENBF  ( ibufr_unit, 'IN', ibufr_unit )
      CALL DATELEN  ( 10 )

c    read the next station report from the input bufr file
c----------------------------------------------------------------------

10    CALL READPB  ( ibufr_unit, subset, idate, ierrpb )

      if ( debug ) print *, 'next obs type: ', subset(1:6)

      idate00 = idate/100
      hour01 = idate - idate00*100
      if( hour01 .eq. 0.0 ) hour01 = 24
      if ( debug ) print*, 'idate= ', idate, idate00, hour01

      IF ( ierrpb .eq. -1 )  STOP

      iprof = iprof + 1    !! station number

c    set the appropriate output file unit number
c----------------------------------------------------------------------

      ii = 1
      found = .false.
      DO WHILE  ( ( .not. found ) .and. ( ii .le. NFILO ) )
        IF( subset(1:6) .eq. filo ( ii ) )  THEN
          found = .true.
          iuno = iunso ( ii )

          time0 = hdr(4)
          stat_elv = hdr(5)

          if ( debug ) print*, 'idate= ', idate, idate00, hour01, time0
          hour01 =  hour01 + time0

          if ( USE_THIS_DATA_REAL(hdr(6),otype_use,inum_otype) .and.
     &                    debug ) then  !  print data read in

            WRITE  ( UNIT = iuno, FMT = '("#", 132("-") )' )
            WRITE  ( UNIT = iuno, FMT = '("#", A3,6x,a6,1x,a6,3x,
     &             a5,5x,a4,3x,a5,4x,a4,3x,a4,1x,a3,1x,a3,2x,a6,
     &             6x,a3,6x,a3,4x,a5,1x,a8,2x,a7,6x,a3 )' )
     &             'SID', 'XOB', 'YOB','DHR','ELV','TYP','T29','ITP',
     &             'lev','var','OB','qm', 'pc', 'rc', 'fc', 'an', 'err'

            WRITE  ( UNIT = iuno, FMT = '("#", 132("-") )' )

          endif

        ELSE

          ii = ii + 1

        END IF

      END DO

c    check the observation time, skip if outside obs. window
c    the prepbufr_03Z.f version tests hour01, this one tests abs(time0)
c----------------------------------------------------------------------
      IF ( ( .not. found ) .and. ( ierrpb .eq. 0 ) )  THEN
        if ( debug ) print*, 'record found w/no label match, val was: ',
     &            subset(1:6)
        GO TO 200
      ENDIF
      IF ( subset(1:6).eq.'SATWND' ) then
        IF ( abs(time0) .gt. obs_window_cw ) then
          if (debug) print*, 'satwnd outside time window, diff was: ',
     &              abs(time0)
          GO TO 200
        ENDIF
      ELSE IF (subset(1:6).eq.'ADPUPA'.or.
     &         subset(1:6).eq.'SATEMP' ) THEN
        IF ( abs(time0) .gt. obs_window_upa ) THEN
         if ( debug ) print*,
     &              'upper-air outside time window, diff was: ',
     &              abs(time0)
          GO TO 200
        ENDIF
      ELSE IF (subset(1:6).eq.'AIRCFT'.or.subset(1:6).eq.'AIRCAR') THEN
        IF ( abs(time0) .gt. obs_window_air ) THEN
          if ( debug ) print*,
     &              'aircraft outside time window, diff was: ',
     &              abs(time0)
          GO TO 200
	ENDIF
      ELSE IF (subset(1:6).eq.'ADPSFC'.or.subset(1:6).eq.'SFCSHP') THEN
        IF ( abs(time0) .gt. obs_window_sfc ) THEN
          if ( debug ) print*,
     &              'surface outside time window, diff was: ',
     &              abs(time0)
          GO TO 200
	ENDIF
      ENDIF

c    place the location data in the appropriate array
c----------------------------------------------------------------------

      if( hdr(2) .ge. 360.0) hdr(2) = hdr(2) - 360.0
      if( hdr(2) .lt.   0.0) hdr(2) = hdr(2) + 360.0

      time = hdr(4)

      tdata(2) = hdr(2)
      tdata(3) = hdr(3)
      tdata(6) = 0
      tdata(7) = iprof + 1000000
      tdata(8) = hour01
         ttype = hdr(6)

      pdata(2) = hdr(2)
      pdata(3) = hdr(3)
      pdata(6) = 0
      pdata(7) = iprof + 3000000
      pdata(8) = hour01
         ptype = hdr(6)

      qdata(2) = hdr(2)
      qdata(3) = hdr(3)
      qdata(6) = 0
      qdata(7) = iprof + 5000000
      qdata(8) = hour01
         qtype = hdr(6)

      zdata(2) = hdr(2)
      zdata(3) = hdr(3)
      zdata(6) = 0
      zdata(7) = iprof + 7000000
      zdata(8) = hour01
         ztype = hdr(6)

      udata(2) = hdr(2)
      udata(3) = hdr(3)
      udata(7) = iprof + 2000000
      udata(8) = hour01           !! time

      vdata(2) = hdr(2)
      vdata(3) = hdr(3)
      vdata(7) = iprof + 9000000
      vdata(8) = hour01

c    u and v wind data are merged with T array for some obs.
c    types.  Account for this and place in appropriate reports
c----------------------------------------------------------------------

      if ( hdr(6) .eq. 120.0 .or. hdr(6) .eq. 130.0 .or.
     &     hdr(6) .eq. 131.0 .or. hdr(6) .eq. 132.0 .or.
     &     hdr(6) .eq. 133.0 ) then
        wtype=hdr(6)+100
      else
        wtype=hdr(6)
      endif

c    check to see if this observation type is desired
      if ( .NOT. USE_THIS_DATA_REAL(real(hdr(6)),otype_use,inum_otype)
     &         ) then
        if (debug) print *, 'this obs type not in use-list, num was: ',
     &       hdr(6)
        GO TO 200
      endif

      DO lv = 1, nlev  !  loop over all levels in the report

c    find out the event before virtural T and Q check step (pc=8.0)
c----------------------------------------------------------------------

        j2t = 1
        do ievent =1, MXR8VN
          if(evns(3,lv,ievent, 3) .lt. 8.0) then
            j2t = ievent
            exit
          endif
        enddo

        j2q = 1
        do ievent =1, MXR8VN
          if(evns(3,lv,ievent, 2) .lt. 8.0) then
            j2q = ievent
            exit
          endif
        enddo

        DO kk = 1, MXR8VT

c     for T & Q, the events (pc>=8) are virtual T step. The events (pc=6, 1) are dry T
c     set qm of j2t/j2q event according to 1st event
c
          if((subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'ADPSFC' .or.
     &        subset(1:6).eq.'SATEMP' .or.
     &        subset(1:6).eq.'SFCSHP') .and.   var(kk).eq. 'T') then
            if(j2t.ne.1) evns(7,lv, j2t, kk) = evns(7,lv, 1, kk)
            jj = j2t
          else if((subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'ADPSFC'
     &      .or.   subset(1:6).eq.'SFCSHP') .and.  var(kk).eq.'Q') then
            if(j2q.ne.1) evns(7,lv, j2q, kk) = evns(7,lv, 1, kk)
            jj = j2q
          else
            jj = 1
          endif
c
          WRITE( UNIT = outstg, FMT = '( A8, 1X,
     &          2F7.2, 1X, F7.3, 1X, 2F8.1, 1X, F7.1, 1X,
     &          F6.1, I4, 1X, A2, 7(1X,F8.1) )' )
     &     (hdr(ii), ii =1, 8),lv, var(kk), (evns(ii,lv,jj,kk), ii=1,7)

          DO mm = 1, 200
            IF  ( outstg (mm:mm) .eq. '*' ) outstg (mm:mm) = ' '
          END DO

          if (debug) write(iuno, '(A)') trim(outstg)
        END DO

        processed = .false.

c    set up the temperature observation data, use the j2t event or 1 if T
c----------------------------------------------------------------------
        if(subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'ADPSFC' .or.
     &     subset(1:6).eq.'SATEMP' .or.
     &     subset(1:6).eq.'SFCSHP' ) then

           toe = evns(7, lv, j2t, 3)
           tob = evns(1, lv, j2t, 3) + 273.16
           tqm = evns(2, lv, j2t, 3)
          pc_t = evns(3, lv, j2t, 3)

        else  ! use the first event

           toe = evns(7, lv, 1, 3)
           tob = evns(1, lv, 1, 3) + 273.16
           tqm = evns(2, lv, 1, 3)
          pc_t = evns(3, lv, 1, 3)

        endif

        if (pc_t < 0 .or. pc_t > 99) pc_t = 99
        if(tqm .eq. 0) toe = toe*0.9
        if(tqm .eq. 3) toe = toe*1.2

c    set up the moisture observation data, use the j2t event for RH and
c     compute specific humidity in g/kg
c----------------------------------------------------------------------

        if(subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'ADPSFC' .or.
     &     subset(1:6).eq.'SFCSHP' ) then

          qoe  = evns(7, lv, j2q, 2) * 0.1     ! g/kg
          qob  = evns(1, lv, j2q, 2) * 1.0e-3  ! g/kg
          qqm  = evns(2, lv, j2q, 2)
          pc_q = evns(3, lv, j2q, 2)

        else

          qoe  = evns(7, lv, 1, 2) * 0.1       ! g/kg
          qob  = evns(1, lv, 1, 2) * 1.0e-3    ! g/kg
          qqm  = evns(2, lv, 1, 2)
          pc_q = evns(3, lv, 1, 2)

        endif
        if (pc_q < 0 .or. pc_q > 99) pc_q = 99
        if(qqm .eq. 0) qoe = qoe*0.9
        if(qqm .eq. 3) qoe = qoe*1.2

c    set up the pressure observation data, units are mb
c----------------------------------------------------------------------

        poe  = evns(7, lv, 1, 1)
        pob  = evns(1, lv, 1, 1)
        pqm  = evns(2, lv, 1, 1)
        pc_p = evns(3, lv, 1, 1)

        if (pc_p < 0 .or. pc_p > 99) pc_p = 99
        if(pqm .eq. 0) poe = poe*0.9
        if(pqm .eq. 3) poe = poe*1.2
        ppb = evns(1, lv, 1, 1)            ! mb
        !zob = evns(1, lv, 1, 4)            ! m

c    set up the height observation data, units are m
c----------------------------------------------------------------------

        zoe  = evns(7, lv, 1, 4)
        zob = evns(1, lv, 1, 4)            ! m
        zqm  = evns(2, lv, 1, 4)
        pc_z = evns(3, lv, 1, 4)

        if (pc_z < 0 .or. pc_z > 99) pc_z = 99
        if(zqm .eq. 0) zoe = zoe*0.9
        if(zqm .eq. 3) zoe = zoe*1.2

c    set up the wind observation data, units are m/s
c----------------------------------------------------------------------

        uoe  = evns(7, lv, 1, 5)
        uob  = evns(1, lv, 1, 5)
        uqm  = evns(2, lv, 1, 5)
        pc_u = evns(3, lv, 1, 5)
        if(uqm .eq. 0) uoe = uoe*0.9
        if(uqm .eq. 3) uoe = uoe*1.2

        voe  = evns(7, lv, 1, 6)
        vob  = evns(1, lv, 1, 6)
        vqm  = evns(2, lv, 1, 6)
        pc_v = evns(3, lv, 1, 6)

        if (pc_u < 0 .or. pc_u > 99) pc_u = 99
        if (pc_v < 0 .or. pc_v > 99) pc_v = 99
        if(vqm .eq. 0) voe = voe*0.9
        if(vqm .eq. 3) voe = voe*1.2

c    write out temperature observation from ADPUPA, AIRCAR, AIRCFT
c----------------------------------------------------------------------
        if (subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'AIRCAR' .or.
     &     subset(1:6).eq.'SATEMP' .or.
     &                                   subset(1:6).eq.'AIRCFT') then
          if (use_this_data_int(tqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        toe .lt. 1.e9 .and. tob .lt. 1.e9                  ) then

            tdata(1) = toe
            tdata(4) = ppb
            tdata(5) = tob

c    in some old files this appears to be out of range
c    and it seems to be unused in converting to an obs_seq so
c    i feel ok setting it to something that will fit in an I2 field.
            if (pc_t < 0 .or. pc_t > 99) pc_t = 99
            write(lunobs, 800) tdata, ttype, tqm, subset(1:6), pc_t
            processed = .true.

          else
            if ( debug ) print *, 'skip temp, toe,tob,tqm,pqm = ',
     &           toe, tob, tqm, pqm
          endif   ! if use_this_data
        endif

c    write out temperature observation of SFCSHP, ADPSFC
c----------------------------------------------------------------------

        if (subset(1:6).eq.'SFCSHP' .or. subset(1:6).eq.'ADPSFC') then

          if (use_this_data_int(tqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        tob .lt. 1.0e9 ) then

            if ( ttype .gt. 200 )  ttype = ttype - 100
            if ( ttype .eq. 184 )  ttype = 183

            tdata(1) = toe
            if ( toe .ge. 1.e9 ) tdata(1) = land_temp_error
            tdata(4) = zob
            tdata(5) = tob

            if (pc_t < 0 .or. pc_t > 99) pc_t = 99
            write(lunobs, 800) tdata, ttype, tqm, subset(1:6), pc_t
            processed = .true.

           else
             if ( debug ) print *, 'skip temp; tob,tqm,pqm = ',
     &        tob, tqm, pqm
           endif   ! if use_this_data
        endif

c    write out moisture observation from ADPUPA
c----------------------------------------------------------------------

        if ( subset(1:6) .eq. 'ADPUPA' ) then

          if (use_this_data_int(qqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        qoe .lt. 1.e9 .and. qob .lt. 1.e9             ) then

c           compute the error of specific moisture based on RH obs error

            qoe_rh = qoe

            es   = es0 * (tob / t0) ** fact1 * exp(fact2*(fact3-1./tob))
            qsat = eps * es / (pob - omeps * es)
            qoe  = max(0.1, qoe * qsat * 1000.0) ! to g/kg, set min value

            if ( debug ) print*, 'es= ', es, tob-273.16, qoe

            if( .not. use_this_data_int(tqm,qctype_use,inum_qctype))then
              if ( debug ) print*, 'upa bad = ', qoe  ! the T obs cannot be used for qoe
              qoe = 1.0e10
            endif

            if (qoe .lt. 9.9) then  ! skip large qoe obs.

              if (pc_q < 0 .or. pc_q > 99) pc_q = 99
              processed = .true.

              ! specific humidity
              qdata(1) = qoe
              qdata(4) = ppb
              qdata(5) = qob
              qdata(6) = 0.0
              write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

              !  write RH data to file
              qdata(1) = qoe_rh
              qdata(5) = qob / (qsat * 1000.0)
              qdata(6) = 1.0
              write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

              !  write dew point data to file
              eo = (qsat * (1.0-qoe_rh) * pob) /
     &             (eps + omeps * qsat * (1.0-qoe_rh))
              qdata(1) = tob - 1.0 / (1.0/t0 - (rv/l0) * log(eo/es0))
              eo = (qob * pob * 0.001) / (eps + omeps * qob * 0.001)
              qdata(5) = 1.0 / (1.0/t0 - (rv/l0) * log(eo/es0))
              qdata(6) = 2.0
              write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

            endif  ! qoe not too large

         else
            if ( debug ) print *, 'skip moist, qoe,qob,qqm,pqm = ',
     &          qoe, qpb, qqm, pqm
         endif   ! if use_this_data

        endif

c    write out moisture observation from SFCSHP ADPSFC
c----------------------------------------------------------------------

        if (subset(1:6).eq.'SFCSHP' .or. subset(1:6).eq.'ADPSFC') then

          if (use_this_data_int(qqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        qob .lt. 1.0e9                                ) then

           es   = es0 * (tob / t0) ** fact1 * exp(fact2*(fact3-1./tob))
           qsat = eps * es / (pob - omeps * es)
           if (qoe .gt. 1.e9) qoe = land_moist_error
           qoe_rh = qoe
           qoe  = max(0.1, qoe * qsat * 1000.0) ! to g/kg, set min value

           if( .not. use_this_data_int(tqm,qctype_use,inum_qctype)) then
             if ( debug ) print*, 'surface bad = ', qoe
c             ! the T obs cannot be used for qoe
             qoe = 1.0e10
           endif

           if ( qtype .gt. 200 )  qtype = qtype - 100
           if ( qtype .eq. 184 )  qtype = 183

           if(qoe .lt. 9.9) then   ! skip large qoe obs

             if (pc_q < 0 .or. pc_q > 99) pc_q = 99
             processed = .true.

             ! specific humidity
             qdata(1) = qoe
             qdata(4) = zob
             qdata(5) = qob
             qdata(6) = 0.0
             write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

             !  write RH data to file
             qdata(1) = qoe_rh
             qdata(5) = qob / (qsat * 1000.0)
             qdata(6) = 1.0
             write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

             !  write dew point data to file
             eo = (qsat * (1.0-qoe_rh) * pob) /
     &            (eps + omeps * qsat * (1.0-qoe_rh))
             qdata(1) = tob - 1.0 / (1.0/t0 - (rv/l0) * log(eo/es0))
             eo = (qob * pob * 0.001) / (eps + omeps * qob * 0.001)
             qdata(5) = 1.0 / (1.0/t0 - (rv/l0) * log(eo/es0))
             qdata(6) = 2.0
             write(lunobs, 800) qdata, qtype, qqm, subset(1:6), pc_q

           endif  ! if qoe not too large

          else
               if ( debug ) print *, 'skip moist, qob,qqm,pqm = ',
     &          qob, qqm, pqm
          endif   ! if use_this_data

        endif

c    write out surface pressure observations
c----------------------------------------------------------------------

        if (subset(1:6).eq.'ADPUPA' .or.
     &      subset(1:6).eq.'SFCSHP' .or. subset(1:6).eq.'ADPSFC') then

          if(use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &       poe.lt.1.e9 .and. pob.lt.1.e9 .and. stat_elv.eq.zob) then

            if ( ptype .gt. 200 )  ptype = ptype - 100
            if ( ptype .eq. 184 )  ptype = 183

            pdata(1) = poe
            pdata(4) = zob
            pdata(5) = pob

            if (pc_p < 0 .or. pc_p > 99) pc_p = 99
            write(lunobs, 800) pdata, ptype, pqm, subset(1:6), pc_p
            processed = .true.

           else
             if ( debug ) print *, 'skip press, poe,pob,pqm,z,s = ',
     &          poe, pob, pqm, zob, stat_elv
          endif   ! if use_this_data

        endif

c    write out wind observation of ADPUPA, AIRcraft, SATWND
c----------------------------------------------------------------------

        if (subset(1:6).eq.'ADPUPA' .or. subset(1:6).eq.'AIRCAR' .or.
     &      subset(1:6).eq.'AIRCFT' .or. subset(1:6).eq.'SATWND') then

          if(use_this_data_int(uqm,qctype_use,inum_qctype) .and.
     &       use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &       uoe .lt. 1.e9 .and. uob .lt. 1.e9             .and.
     &       voe .lt. 1.e9 .and. vob .lt. 1.e9             ) then

            udata(1) = uoe
            udata(4) = ppb
            udata(5) = uob
            udata(6) = vob

            vdata(1) = voe
            vdata(4) = ppb
            vdata(5) = vob
            vdata(6) = uob

            if (pc_u < 0 .or. pc_u > 99) pc_u = 99
            if (pc_v < 0 .or. pc_v > 99) pc_v = 99
            write(lunobs, 800) udata, wtype, uqm, subset(1:6), pc_u
            write(lunobs, 800) vdata, wtype, vqm, subset(1:6), pc_v
            processed = .true.

          else
             if ( debug ) print *, 'skip wind, uoe,uob,uqm,pqm = ',
     &          uoe, uob, uqm, pqm
          endif  ! if use_this_data

        endif

c    write out wind observation of SFCSHP ADPSFC
c----------------------------------------------------------------------

        if (subset(1:6).eq.'SFCSHP' .or. subset(1:6).eq.'ADPSFC' ) then

          if (use_this_data_int(uqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        uob .lt. 1.e9                                 ) then

            udata(1) = uoe
            if ( uoe .ge. 1.e9 ) udata(1) = land_wind_error
            udata(4) = zob
            udata(5) = uob
            udata(6) = vob

            vdata(1) = voe
            if ( voe .ge. 1.e9 ) vdata(1) = land_wind_error
            vdata(4) = zob
            vdata(5) = vob
            vdata(6) = uob

            if (pc_u < 0 .or. pc_u > 99) pc_u = 99
            if (pc_v < 0 .or. pc_v > 99) pc_v = 99
            write(lunobs, 800) udata, wtype, uqm, subset(1:6), pc_u
            write(lunobs, 800) vdata, wtype, vqm, subset(1:6), pc_v
            processed = .true.

          else
             if ( debug ) print *, 'skip wind, uob,uqm,pqm = ',
     &          uob, uqm, pqm
          endif   ! if use_this_data

        endif

c    write out geopotential height observation of ADPUPA
c----------------------------------------------------------------------

        if (subset(1:6).eq.'ADPUPA') then

          if (use_this_data_int(zqm,qctype_use,inum_qctype) .and.
     &        use_this_data_int(pqm,qctype_use,inum_qctype) .and.
     &        zoe. lt. 1.e9 .and. zob .lt. 1.e9            ) then

            zdata(1) = zoe
            zdata(4) = ppb
            zdata(5) = zob

            if (pc_z < 0 .or. pc_z > 99) pc_z = 99
            write(lunobs, 800) zdata, ztype, zqm, subset(1:6), pc_z
            processed = .true.

          else
             if ( debug ) print *,
     &           'skip geopotential height, zoe,zob,zqm,pqm = ',
     &            zoe, zob, zqm, pqm
          endif   ! if use_this_data

        endif

c----------------------------------------------------------------------
       if (.not. processed) then
         if (debug)print*, 'bot of loop w/o processing, obs was: ',
     &              subset(1:6), ' lv =', lv
       endif

      enddo
200    continue

      IF ( ierrpb .eq. 0 )  GO TO 10

800   format(f5.2,2f9.4,e12.5,f7.2,f7.2,f9.0,f7.3,i4,i2,1x,a6,i2)
      STOP
      END
C-----------------------------------------------------------------------
        SUBROUTINE READPB  ( lunit, subset, idate, iret )
C
C      This subroutine will read and combine the mass and wind subsets
C      of the next station report in the prepbufr file.  It is styled
C      after entry point READNS, and it only requires the prepbufr file
C      to be opened for reading with OPENBF.  The combined station
C      report is returned to the caller in COMMON /PREPBC/.
C      This common area contains the number of levels in the report,
C      a one dimensional array with the header information, and a four
C      dimensional array containing all events from the variables POB,
C      QOB, TOB, ZOB, UOB, and VOB for the report.
C
C      The header array contains the following list of mnemonics:
C
C      SID XOB YOB DHR ELV TYP T29 ITP
C
C      The 4-D array of data, EVNS ( ii, lv, jj, kk ), is indexed
C      as follows:
C
C      "ii" indexes the event data types; these consist of:
C          1) OBservation
C          2) Quality Mark
C          3) Program Code
C          4) Reason Code
C          5) ForeCast value
C          6) ANalysed value
C          7) observation error
C      "lv" indexes the levels of the report
C      "jj" indexes the event stacks
C      "kk" indexes the variable types (p,q,t,z,u,v)
C
C      Note that the structure of this array is identical to one
C      returned from UFBEVN, with an additional (4th) dimension to
C      include the six variable types into the same array.
C
C      The return codes are as follows:
C      iret =  0 - normal return
C           =  1 - the station report within COMMON /PREPBC/ contains the
C                  last available subset from within the prepbufr file
C           = -1 - there are no more subsets available from within the
C                  prepbufr file
C
        INCLUDE  'prepbufr.prm'
        CHARACTER*(*)   subset
        CHARACTER*(MXSTRL)      head
     +          / 'SID XOB YOB DHR ELV TYP T29 ITP' /
Cliu
        CHARACTER*(MXSTRL)      ostr ( MXR8VT )
     +          / 'POB PQM PPC PRC PFC PAN POE',
     +            'QOB QQM QPC QRC QFC QAN QOE',
     +            'TOB TQM TPC TRC TFC TAN TOE',
     +            'ZOB ZQM ZPC ZRC ZFC ZAN ZOE',
     +            'UOB WQM WPC WRC UFC UAN WOE',
     +            'VOB WQM WPC WRC VFC VAN WOE'  /
C
        REAL*8          hdr2 ( MXR8PM ),
     +                  evns2 ( MXR8PM, MXR8LV, MXR8VN, MXR8VT )
        REAL*8          r8sid, r8sid2, pob1, pob2
        CHARACTER*8     csid, csid2, subst2
        LOGICAL         match / .true. /
        EQUIVALENCE     ( r8sid, csid ), ( r8sid2, csid2 )
C
        SAVE            match, subst2, idate2
C-----------------------------------------------------------------------

Cnsc
 1000   continue
Cnsc

        iret = 0
C
C      If the previous call to this subroutine did not yield matching
C      mass and wind subsets, then READNS is already pointing at an
C      unmatched subset.  Otherwise, call READNS to advance the subset
C      pointer to the next subset.
C
        IF  ( match )  THEN
            CALL READNS  ( lunit, subset, idate, jret )
            IF  ( jret .ne. 0 )  THEN
                iret = -1
                RETURN
            END IF

cliu   select data type
       if(subset .ne. 'ADPUPA' .and. subset .ne. 'AIRCAR' .and.
     &    subset .ne. 'SATWND' .and. subset .ne. 'AIRCFT' .and.
     &    subset .ne. 'SATEMP' .and.
     &    subset .ne. 'ADPSFC' .and. subset .ne. 'SFCSHP') go to 1000
cliu

        ELSE
            subset = subst2
            idate = idate2
        END IF
C
C      Read the HDR and EVNS data for the subset that is currently
C      being pointed to.
C
        CALL UFBINT  ( lunit, hdr, MXR8PM, 1, jret, head )
        DO ii = 1, MXR8VT
          CALL UFBEVN ( lunit, evns ( 1, 1, 1, ii ), MXR8PM, MXR8LV,
     +                     MXR8VN, nlev, ostr (ii) )
        END DO
C
C      Now, advance the subset pointer to the following subset and
C      read its HDR data.
C
Cnsc
 2000   continue
Cnsc
        CALL READNS  ( lunit, subst2, idate2, jret )
        IF  ( jret .ne. 0 )  THEN
            iret = 1
            RETURN
        END IF
C
cliu   select data type
       if(subst2.ne. 'ADPUPA' .and. subst2.ne. 'AIRCAR' .and.
     &    subst2.ne. 'SATWND' .and. subst2.ne. 'AIRCFT' .and.
     &    subst2.ne. 'SATEMP' .and.
     &    subst2.ne. 'ADPSFC' .and. subst2.ne. 'SFCSHP') go to 2000
c        ! careful about 2000
cliu

        CALL UFBINT  ( lunit, hdr2, MXR8PM, 1, jret, head )
C
C      Check whether these two subsets have identical SID, YOB, XOB,
C      ELV, and DHR values.  If so, then they are matching mass and
C      wind subsets for a single station report.
C
        match = .true.
C
        IF  ( subset .ne. subst2 )  THEN
            match = .false.
            RETURN
        END IF
C
        r8sid = hdr (1)
        r8sid2 = hdr2 (1)
        IF  ( csid .ne. csid2 )  THEN
            match = .false.
            RETURN
        END IF
C
        DO ii = 2, 5
            IF  ( hdr (ii) .ne. hdr2 (ii) )  THEN
                match = .false.
                RETURN
            END IF
        END DO
C
C      Read the EVNS data for the second of the two matching subsets.
C
        DO ii = 1, MXR8VT
          CALL UFBEVN ( lunit, evns2 ( 1, 1, 1, ii ), MXR8PM, MXR8LV,
     +                     MXR8VN, nlev2, ostr (ii) )
        ENDDO
C
C      Combine the EVNS data for the two matching subsets into a
C      single 4-D array.  Do this by merging the EVNS2 array into the EVNS array.
C
        DO 10 lv2 = 1, nlev2
            DO lv = 1, nlev
                pob1 = evns ( 1, lv, 1, 1 )
                pob2 = evns2 ( 1, lv2, 1, 1 )
                IF  ( pob1 .eq. pob2 )  THEN
C
                  DO kk = 1, MXR8VT
                    DO jj = 1, MXR8VN
                      DO ii = 1, MXR8PM
                        IF  ( evns ( ii, lv, jj, kk ) .eq. R8BFMS ) THEN
                          evns ( ii, lv, jj, kk ) =
     +                          evns2 ( ii, lv2, jj, kk )
                        END IF
                      END DO
                    END DO
                  END DO
                  GO TO 10
                ELSE IF ( (pob2 .gt. pob1) .or. (lv .eq. nlev) ) THEN
C
                  nlev = nlev + 1
                  DO kk = 1, MXR8VT
                    DO jj = 1, MXR8VN
                      DO ii = 1, MXR8PM
                        evns ( ii, nlev, jj, kk ) =
     +                        evns2 ( ii, lv2, jj, kk )
                      END DO
                    END DO
                  END DO
                  GOTO 10
                END IF
            END DO
   10   END DO
C
        RETURN
        END

C-----------------------------------------------------------------------
        FUNCTION USE_THIS_DATA_INT(ifld, ufld, nfld)
C
C       ifld is the integer value to test, ufld is an array of
C       integers 'nfld' long.  if ifld is found in any of the
C       ufld entries, this function returns true, else false.
C       if nfld is 0, it bypasses the search and always
C       returns true; accepting any value.
C
        integer, intent(in) :: nfld, ifld, ufld(nfld)

        integer             :: nn
        logical             :: use_this_data_int

C       bypass test - always return true
        if (nfld == 0) then
          use_this_data_int = .TRUE.
          return
        endif

C       normal search code, look for match in list
        use_this_data_int = .FALSE.
        DO nn = 1, nfld
          IF ( ifld .eq. ufld(nn) ) THEN
            use_this_data_int = .TRUE.
            RETURN
          END IF
        END DO

        RETURN
        END

C-----------------------------------------------------------------------
        FUNCTION USE_THIS_DATA_REAL(rfld, ufld, nfld)
C
C       ifld is the real value to test, ufld is an array of
C       reals 'nfld' long.  if ifld is found in any of the
C       ufld entries, this function returns true, else false.
C       if nfld is 0, it bypasses the search and always
C       returns true; accepting any value.
C

        integer, intent(in) :: nfld
        real, intent(in)    :: rfld, ufld(nfld)

        integer             :: nn
        logical             :: USE_THIS_DATA_REAL

C       bypass test - always return true
        if (nfld == 0) then
          use_this_data_real = .TRUE.
          return
        endif

C       normal search code, look for match in list
        USE_THIS_DATA_REAL = .FALSE.
        DO nn = 1, nfld
          IF ( rfld .eq. ufld(nn) ) THEN
            USE_THIS_DATA_REAL = .TRUE.
            RETURN
          END IF
        END DO

        RETURN
        END

C <next few lines under version control, do not edit>
C $URL$
C $Id$
C $Revision$
C $Date$

max_stars_repo_name: nancycollins/DART

max_issues_repo_name: nancycollins/DART

max_forks_repo_name: nancycollins/DART


text:
A Gelato and Candy shop, now serving bulk and gourmet candies and hot chocolate. On their opening day their computer system was down  resulting in free gelato for all!  It has a window into the Varsity Theatre as well as a regular storefront.  The cafe is owned by Sinisa Novakovic, who also owns Mishkas Cafe and coowns the Varsity with Users/JonFenske

This cafe serves Caffe Classico Gelato (dairy) and Sorbetto (nondairy).  Both are also available at some local grocery stores.

The smallest size gelato costs $3.15 and the larger size increases the cost by a dollar ($4.15, $5.15),  it used to cost an extra $0.50 if you bought it through the theatre but that fee no longer exists.

The walls are painted concrete, and as a result, voices echo loudly and the acoustics are poor inside the small store space.
I havent noticed this, and the walls have some sort of plaster sprayed on them for the texture  Users/StevenDaubert
Now that we sell candy weve got shelves and candy jars all over the walls that help with the acoustics of the shop. Users/MeganYourFace

   Gotta love ice cream, gelato, or frozen yogurt on a hot Summer day!

20060724 23:55:19 nbsp Excellent gelato! Finally someone opened a gelato shop in Davis. And in time for the heat wave too. The only nagging question is: what new flavors are they planning for the second cooler? P.S. hazlenut chocolate... mmm... Users/LiRic
  Gelato has been available in Davis for more than a decade I believe. Konditorei in East Davis not only has baked goods but a glass case of gelato. It would be interesting to hear how the two locations compare in quality. KJM
     Yes, Konditorei has had gelato for a while now. Unfortunately their paltry number of flavors hardly qualifies them to be Davis sole gelato store. I would soonerand dogo to Berkeley for gelato than eat the gelato at Konditorei. Thank goodness Icekrimski Cafe opened. How much is the gelato at Konditorei? I never bothered to check.
   Zias sells gelato...
     Is the gelato at Zias any good? I forgot that they had it until you mentioned it. I think I remember that they only had a few unspecial flavors, though. I think I shouldve specified that its good to finally have a dedicated gelato store in Davis. BTW, thanks to whoever it was who posted pics. Its good to see that they stocked their second cooler.
      No problemo... I was standing there today, and my camera was in my pocket, and I thought, oh duh! I am visible in the reflection in the upper right pic. KJM
     Zias has been selling Gelato Classico brand for over a year now. Same great taste as here, and they seem to change flavors often. They also sell pints in another freezer. Any idea if you can get pints here?
        Does gelato taste good when brought home? If anyones ever bought a pint, is it hard to keep soft? Do you need to lower your fridge temperature?
           The Gelato is happiest at 5 to 9 degrees. But a freezer will do, just keep it away from the walls. The optimal consistency is of thick peanut butter.  Users/StevenDaubert

20060726 14:55:13 nbsp Gelato is wonderful stuff.  Should this be the Gelatski Cafe? Users/JaimeRaba

20060727 17:14:37 nbsp Computer system down  free gelato?? What, so before computers were invented there were no monetary transactions anywhere? Are people that bad at math? Users/BarnabasTruman
  They also were not told the prices for the Gelato sizes on the first day, so it turned into an impromptu promotional freebie.  KJM
  Oh well, I guess I shouldnt be complaining about free ice cream. :)  Users/BarnabasTruman

20060728 21:06:54 nbsp Score! The management likes Apple products, I guess. On opening night, there was a Mac mini and now theres a Bose Sounddock for the iPod. Props to the owner. Users/LiRic

A few plants strategically placed in the corners would help with the excessive echoing, along with acoustic tile on the ceiling to break up the standing waves. Users/KarlMogel

20060731 23:31:57 nbsp I think they have fancy ice cream there too, which is really tasty. I got the dulce de leche. Yummie, :D. Theyre very nice, its good to see theyre busy! Users/ChristyMarsden

20060823 13:20:10 nbsp As an ice cream guru here is the skinny... the gelato is YUMMY!  wonderful!  I love it!  It is the same brand they sell at Zias but Zias isnt open at night.  They also sell pints of it at the Coop. the ice cream is double rainbow brand and is tasty. Also pints at the Coop. The shop is really well placed and done.  I really like the community it adds to downtown.  Another great reason to walk around and mingle the people on the planet you live with. Thanks for opening a great new hang out for ice cream in downtown!!   Users/HeatherMartin

20070228 23:53:35 nbsp A foodie told me that dots you see in the vanilla bean and the double espresso bean shows that it comes from real beans, not extracts. Users/StevenDaubert

2007314 I just stopped by here for the first time tonight, and it made my evening. I didnt have cash when I walked in (they only take cash), but after a few samples I rode my bike to the bank to get money. Yeah, thats how good it was. The fruit flavors tasted EXACTLY like the fruit theyre named after (tangerine, mango, lemon, coconut). Reminds me of France. And the chocolate flavors (chocolate mint, chocolate hazlenut, chocolate malt) were so rich and perfect. Reminds me of Italy. I love trying new and slightly strange things, so I LOVED the ginger flavor (try it!) and the blackberry cabernet sorbetto. Wow. Red wine flavored sorbet. Why didnt they think of that sooner? It tasted like a luscious grownup popsicle. So I guess what Im trying to say is that you really should try this place. Just remember to bring cash! Users/MareGiudice

20070322 16:16:15 nbsp I love gelato so I decided to try this place. It looks nice and the variety of flavors were okay not that great. Pricing was a little steep for me compared to the gelato stores I went to in the Bay Area. It was decent gelato but not the best. Customer service was terrible, the server kept avoiding me and my friend (guess he didnt like Asians). Basically wanted us out quickly and didnt even offer samples to us. From this I wouldnt come back here as often.  Users/BrandonToo

   Any evidence for that charge of racism?  Users/JimEvans
   I assume the idea is that he and his friend are both Asian. Not much evidence, I agree, but there you are.  Users/NickSchmalenberger
  Just seemed that the server didnt want to serve my friend and I. The group before us were treated very nicely but when it came to us he acted like we didnt belong there. We didnt even do anything but smile and greet him. Then a group of guys come around and immediately the server greeted them in a friendly manner and asked them if they wanted any samples. He actually avoided us until I asked him what I wanted. Simply stated I had a bad experience here. And a sidenote is that both groups before and after us were white people so of course it would feel strange as to why we were not treated the same way.  Users/BrandonToo

20070324 23:21:49 nbsp Stumbled upon Icekrimski the other night, had a very pleasant experience.  The two brunettes working were swamped but still able to greet each customer and serve them with efficiency and smiles.  The Gelato was fantastic, Im definatley going to be a regular! Users/JacobWright

20070423 09:47:13 nbsp Great place. Im glad they got rid of the Ice Cream it wasnt needed.
Hope everybody has a chance to try it I will definately be back. Users/KevinRobinson

20070429 18:50:33 nbsp Icekrimskis has really yummy gelato! I went to a huge gelato place several months ago in my hometown of Los Angeles. This place was in the Santa Monica Promenade, and had like a hundred flavors! Still, no matter how many flavors this place had, their gelato was not awesome, and totally paled in comparison to Icekrimski.
Users/ArielaHaro

20070430 17:51:38 nbsp The Pumpkin flavor is simply wonderful. Try it, folks. Users/KarlMogel
  indeed Karl, did you enjoy the nutmeg and the cinnamon at the end to give just a hint of pumpkin pieness?  Users/StevenDaubert

20070506 14:49:06 nbsp With a ton of unique and wonderful options to chose from and the opportunity to sample the most tempting morsels, this is a must try spot in Davis.  The staff is funny and very attentive.  Stopping here leaves a pleasant smile on customers faces with a hint of nutmeg! Users/RachelRagatz

20070519 18:57:58 nbsp When getting gelato tonight there were two people working but only one was actually helping customers.  The guy just sat there talking to friends, I assume, while the girl had to help everyone.  I kept expecting him to stand up and help us, but we ended having to wait a while for the girl to help everyone before us. I dont know if he was on a break or something, but i just thought it was really rude that he wasnt helping her. When we finally got our gelato, it was delicious, just wish the service was more professional.  Users/JoyChu

20070522 07:12:49 nbsp My friend and I went there yesterday. Pretty friendly place. Daubert gave us a lot of samples. Try the tiramisu , dulce de leche or caramello chocolate crunch gelato. I also like mixing the mango sorbetto with any of the fruit flavors. Great service, great gelato! Users/JonasMari

20070523 11:32:48 nbsp yes, as a frequent patron at icekrimski, the blonde guy who works there is very friendly (when he is the only one there); however, it seems like whenever there are two people working, he lets the other person do all the work (i know because i have seen it on more than one occasion, and my friends who come there can verify it as well). i dont mean to cause controversy, but i just call it like i see it. Users/JoyChu

20070524 21:25:33 nbsp Went here today for the first time.  There was only one girl working there, and it wasnt busy at all, so the service was prompt and attentive.  We tried a couple of samples and then I settled on the cherry gelato that had plenty of cherry chunks in it.  It was very good, and I thought the portion size was decent for the price.  My fiancee had half strawberry, half something with a weird name that was effectively cookies and cream.  She was very pleased as well. Users/AlexPomeranz
    Stracciatella?  Users/StevenDaubert
         It was a strange word like that.  Wikipedias definition for http://en.wikipedia.org/wiki/Stracciatella Stracciatella says it can be Vanilla ice cream with chocolate shavings.  Could have been. Users/AlexPomeranz

20070527 17:06:03 nbsp This place is awesome. Come here next time you have a frozen dairy craving. The employees are very friendly and the gelato is excellent. Users/AmitVainsencher

20070528 17:46:35 nbsp I finally went here for the lemon gelato, it was very good but seemingly not the best Ive had (thought its been a while).  Ill have my dad stop by to check it out, he raves about lemon gelato from his old neighborhood in the Bronx. Users/EdHenn

20070608 00:59:14 nbsp Always great service. Always delicious gelato. :) Users/CAkita

20070611 09:42:24 nbsp The gelato here is really good, and refreshing on a hot day.  (So, refreshing just about every day??).  Try the sorbet  its intensely flavored and will knock your socks off.  (Which you dont need anyway, since its so darned hot...) Users/CovertProfessor

20070712 14:34:49 nbsp Im a relatively frequent customer at Icekrimski (I go there about once a week), and its quite tasty.  The mango and guava sorbettos are amazing, and the child size is sufficient simply because of the amount of flavor that it packs.  Of the gelato flavors, Id probably recommend the creme brulee or the hazelnut.  Yummy.

On another note, the service there varied depending on who was working.  The two younger employees (the blonde haired guy and a younger gal) have always been very courteous, but the older lady that was working there (with the useless guy that others have mentioned) wasnt very nice.  In fact, she seemed kind of pissed to serve me and my friends.  Maybe she was annoyed that we got excited over frozen treats. Users/BrianMac

20070929 19:57:04 nbsp I recommend getting a regular sized half mango, half passion fruit sorbetto. The flavors complement each other nicely.  Users/at86

20080201 17:28:41 nbsp does anyone else find it troubling that the employees here touch the cookie that goes into the ice cream, and touch your money, and everyone elses money, and do not wash their hands??  i believe that is a health code violation.  Users/Chantal

20080323 23:59:26 nbsp I like this place, they always pack you selected cup size full of the gelato. The pumpkin gelato is my favorite but another good combination is half chocolate and half banana fosterstastes just like the chocolate dipped frozen bananas from Disneyland.   Users/NtoWolves

20080523 10:50:28 nbsp Went here last night for the first time.  Pleasantly satisfied with service and delicious!  I had the Caramello Chocolate chrunch and the Strawberry Fresca.  The strawberry fresca was really good and my favorite of the two. Users/Jillep

20081021 14:04:21 nbsp Hella tasty.  I would recommend this place to anyone. Users/ChrisWaterstraat

20090124 22:35:51 nbsp They have excellent limoncello poppy seed gelato. However, their hours are kind of wonky. I went in today to get some of my favorite gelato, and they werent open. Their window had a sign posting that their hours were frisun 3:307:30. Pretty lame. I went to BOTH the coop and nugget market to see what flavors of classico they serve and was bummed to find that they pretty much only carried vanilla, tiramisu, a mango kind and maybe one or two others. This place definitely has the best selection, but make sure you go when they are actually open. I think that at least on the weekends, they should be open for longer. Users/BrookeB

20100716 17:20:33 nbsp This is the best gelato with a good atmosphere and location. A little expensive though, but worth it. Users/JasperD

20120127 13:55:29 nbsp Best gelato joint in town. Users/RoamingDavisite

20130106 18:48:20 nbsp They now accept credit and debit! Hooray! Users/ErinHollander

20130522 15:49:45 nbsp Bring back macadamia nut coconut!!! Users/Aaron.Curtin

20140509 02:46:29 nbsp The gelato has a wonderfully thick consistency and retains its flavour superbly, without the taste diminishing as you work your way through a serving. The mango is tart and fruity (like a chilly Gusher, but better!), and the honey lavender is sweet and tastes subtly like tea. Definitely allow the honey lavender the chance to perform its gentle work upon your tasting buds.  Users/JaniceW


text:
!***********************************************************************
!  Integrated Water Flow Model (IWFM)
!  Copyright (C) 2005-2021
!  State of California, Department of Water Resources
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU General Public License
!  as published by the Free Software Foundation; either version 2
!  of the License, or (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!  (http://www.gnu.org/copyleft/gpl.html)
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
!
!  For tecnical support, e-mail: IWFMtechsupport@water.ca.gov
!***********************************************************************
MODULE Budget_Parameters
  IMPLICIT NONE

  PUBLIC

  ! -------------------------------------------------------------
  ! --- DATA TYPES FOR COLUMNS, MAX TITLE LENGTH FOR ASCII FILE OUTPUTS
  ! -------------------------------------------------------------
  INTEGER,PARAMETER :: VR                     = 1   , &     !Volumetric rate
                       VLB                    = 2   , &     !Volume at the beginning of time step
                       VLE                    = 3   , &     !Volume at the end of time step
                       AR                     = 4   , &     !Area
                       LT                     = 5   , &     !Length
                       VR_lwu_PotCUAW         = 6   , &     !Volumetric rate - specific for Potentail CUAW column in Land & Water Use Budget
                       VR_lwu_AgSupplyReq     = 7   , &     !Volumetric rate - specific for Agricultural Supply Requirement column in Land & Water Use Budget
                       VR_lwu_AgShort         = 8   , &     !Volumetric rate - specific for Agricultural Shortage column in Land & WAter Use Budget
                       VR_lwu_AgPump          = 9   , &     !Volumetric rate - specific for Agricultural Pumping column in Land & WAter Use Budget
                       VR_lwu_AgDiv           = 10  , &     !Volumetric rate - specific for Agricultural Diversions column in Land & WAter Use Budget
                       VR_lwu_AgOthIn         = 11  , &     !Volumetric rate - specific for Agricultural Other Inflows column in Land & WAter Use Budget
                       f_iMaxTitleLen         = 1000, &
                       f_iFormatSpecLen       = 500 , &
                       f_iMaxLocationNameLen  = 100 , &
                       f_iColumnHeaderLen     = 100 , &
                       f_iBudgetDescriptorLen = 100 , &
                       f_iPathNameLen         = 80  , &
                       f_iDSSDataUnitLen      = 8
  INTEGER,PARAMETER :: f_iPER_CUM         = 1   , &
                       f_iPER_AVER        = 2
  CHARACTER(LEN=8),PARAMETER :: f_cDataTypes(2)       = ['PER-CUM ','PER-AVER']
  CHARACTER(LEN=9),PARAMETER :: f_cLocationNameMarker = '@LOCNAME@'
  CHARACTER(LEN=6),PARAMETER :: f_cAreaMarker         = '@AREA@'
  CHARACTER(LEN=8),PARAMETER :: f_cLengthUnitMarker   = '@UNITLT@'
  CHARACTER(LEN=8),PARAMETER :: f_cAreaUnitMarker     = '@UNITAR@'
  CHARACTER(LEN=8),PARAMETER :: f_cVolumeUnitMarker   = '@UNITVL@'

END MODULE

max_stars_repo_name: giorgk/IWFM-SAFE

max_issues_repo_name: giorgk/IWFM-SAFE

max_forks_repo_name: giorgk/IWFM-SAFE


text:
*=======================================================================
* File:    tmtest2.f
* Author:  Alex Stivala
* Created: July 2008
*
* Test program for tmatn.f
*
* $Id: tmtss2.f 2968 2009-11-22 01:21:50Z astivala $
*=======================================================================

      program tmtest2

      implicit none

      external tmatn

      integer n1,n2,info
      parameter(n1 = 7, n2 = 7)
      double precision omega1(n1,n1),omega2(n2,n2)
      double precision dmat1(n1,n1),dmat2(n2,n2)
      double precision score
      double precision soln(n1*n2+n1+n2)

* pytableaucreate.py -t dssp -n /var/tmp/d1hhoa_.ent
      data omega1 / 0.00,  1.75, -2.51,  0.93,  1.46,  1.07, -1.67,
     $     1.75,  0.00,  2.01, -1.06, -0.33,  2.18, -1.33,
     $     -2.51, 2.01,  0.00,  2.66, -2.32, -1.59,  1.72,
     $     0.93, -1.06,  2.66,  0.00,  0.73,  1.98, -0.94,
     $     1.46, -0.33, -2.32,  0.73,  0.00,  2.15, -1.20,
     $     1.07,  2.18, -1.59,  1.98,  2.15,  0.00,  2.60,
     $     -1.67, -1.33,  1.72, -0.94, -1.20, 2.60,  0.00 /

* pytableaucreate.py -t dssp -n /var/tmp/d1hhoa_.ent
      data omega2 / 0.00,  1.75, -2.51,  0.93,  1.46,  1.07, -1.67,
     $     1.75,  0.00,  2.01, -1.06, -0.33,  2.18, -1.33,
     $     -2.51, 2.01,  0.00,  2.66, -2.32, -1.59,  1.72,
     $     0.93, -1.06,  2.66,  0.00,  0.73,  1.98, -0.94,
     $     1.46, -0.33, -2.32,  0.73,  0.00,  2.15, -1.20,
     $     1.07,  2.18, -1.59,  1.98,  2.15,  0.00,  2.60,
     $     -1.67, -1.33,  1.72, -0.94, -1.20, 2.60,  0.00 /

      data dmat1 / 0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0 /

      data dmat2 / 0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0,
     $             0, 0, 0, 0, 0, 0, 0 /

      call tmatn(omega1, n1,n1, omega2, n2 , n2, .false., .false.,
     $           dmat1, dmat2,
     $           score, soln, info)

*     score from MATLAB is -16.4933

      if (info .ne. 0) then
         write(*,*),'TABMATN FAILED'
      else
         write(*,900),score
 900     format('score = ',f6.2)
      endif
      end

max_stars_repo_name: stivalaa/qptabsearch

max_issues_repo_name: stivalaa/qptabsearch

max_forks_repo_name: stivalaa/qptabsearch