module read_l2bufr_mod
!$$$   module documentation block
!                .      .    .                                       .
! module:    read_l2bufr_mod module for processing level2 radar bufr files
!   prgmmr: parrish          org: np22                date: 2005-06-06
!
! abstract: Process level 2 radar bufr files, converting radial wind observations
!             to superobs, which are read by subroutine read_superwinds.
!             Because these files can be very big (5-10 Gb), some additions were
!             made to bufrlib to allow use of mpi-io routines when reading the
!             bufr file.  This reduced processing time from 5-8 minutes to
!             60-90 seconds for a 5Gb test file.
!
! program history log:
!   2005-08-01  parrish
!   2006-04-21  parrish, complete rewrite.
!   2006-05-22  parrish, fix bug which causes infinite loop when no data pass initial checks
!   2007-10-24  parrish  add l2superob_only option
!   2009-04-18  woollen  improve mpi_io interface with bufrlib routines
!   2009-11-24  parrish  change time variable from regional_time (passed from gridmod) to
!                          iadate (passed from obsmod), to prevent all radar data being tossed.
!   2011-07-04  todling  - fixes to run either single or double precision
!
! subroutines included:
!   sub initialize_superob_radar - initialize superob parameters to defaults
!   sub radar_bufr_read_all      - process input level 2 bufr file and write out superobs
!
! Variable Definitions:
!   def del_azimuth     - azimuth range for superob box  (default 5 degrees)
!   def del_elev        - elevation angle range for superob box  (default .05 degrees)
!   def del_range       - radial range for superob box  (default 5 km)
!   def del_time        - 1/2 time range for superob box  (default .5 hours)
!   def elev_angle_max  - max elevation angle (default of 5 deg recommended by S. Liu)
!   def minnum                  - minimum number of samples needed to make a superob
!   def range_max       - max radial range to use in constructing superobs  (default 100km)
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

  use kinds, only: r_kind,i_kind
  implicit none

! set default to private
  private
! set subroutines to public
  public :: initialize_superob_radar
  public :: radar_bufr_read_all
! set passed variables to public
  public :: range_max,del_time,l2superob_only,elev_angle_max,del_azimuth
  public :: minnum,del_range,del_elev

  public :: invtllv

  integer(i_kind) minnum
  real(r_kind) del_azimuth,del_elev,del_range,del_time,elev_angle_max,range_max
  logical l2superob_only

contains

  subroutine initialize_superob_radar
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    initialize_superob_radar
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-21  safford -- add subprogram doc block, rm unused uses
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:
!
!$$$
    use constants, only: quarter,half,five
    implicit none

    del_azimuth=five            !   (5 degrees)
    del_elev=quarter            !   (.25 degrees for elevation angle bin)
    del_range=5000._r_kind      !  (5 km)
    del_time=half               !   (hours)
    elev_angle_max=five         !  recommended by S. Liu to avoid heavy convection problems
    minnum=50
    range_max=100000._r_kind    !  (100km)
    l2superob_only=.false.

  end subroutine initialize_superob_radar

  subroutine radar_bufr_read_all(npe,mype)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    radar_bufr_read_all
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-21  safford -- add subprogram doc block, rm unused uses
!   2010-08-30  treadon - changes for reproducibile results
!   2016-08-12  lippi   - Logic for running single radar test.
!   2016-12-14  lippi   - Read 'l2rwbufr' from the OBS_INPUT table rather than
!                         hardcoding that name into the code. Logic to exit code
!                         if the infile (l2rwbufr) is not present. Also, added
!                         new namelist parameter vadwnd_l2rw_qc.
!
!   input argument list:
!     mype     - mpi task id
!     npe
!
!   output argument list:
!
! attributes:
!   language:  f90
!   machine:
!
!$$$  end documentation block

    use kinds, only: r_double,r_quad
    use constants, only: zero,half,one,two,rearth,deg2rad,rad2deg,zero_quad,one_quad
    use mpimod, only: mpi_comm_world,mpi_min,mpi_sum,mpi_real8,ierror,mpi_real16
    use mpimod, only: mpi_max,mpi_integer4
    use obsmod,only: iadate,ndat,dsis,dfile,dtype
    use qcmod, only: vadwnd_l2rw_qc
    use oneobmod, only: lsingleradar,singleradar
    use mpeu_util, only: IndexSet, IndexSort
    use gsi_io, only: verbose
    implicit none

    integer(i_kind),intent(in):: npe,mype

    integer(i_kind),parameter:: max_num_radars=150
    integer(i_kind),parameter:: n_gates_max=4000
    real(r_kind),parameter:: four_thirds = 4.0_r_kind / 3.0_r_kind
    real(r_kind),parameter:: r8     = 8.0_r_kind
    real(r_kind),parameter:: r89_5  = 89.5_r_kind
    real(r_kind),parameter:: r90    = 90.0_r_kind
    real(r_kind),parameter:: r125   = 125.0_r_kind
    real(r_kind),parameter:: r250   = 250.0_r_kind
    real(r_kind),parameter:: r360   = 360.0_r_kind
    real(r_kind),parameter:: r720   = 720.0_r_kind
    real(r_kind),parameter:: r99999 = 99999._r_kind
    real(r_kind),parameter:: rinv60 = 1.0_r_kind/60.0_r_kind
    real(r_double),parameter:: r1e5_double = 1.0e5_r_double

    integer(i_kind) nazbin,nrbin,nelbin
    integer(i_kind) i,ibyte,idate,inbufr,iret,isubset,krad,levs,lundx,n_gates
    integer(i_kind) k,iii,j,jj,jjall,jjj,numzzzz,num_radars_0,iloc
    integer(i_kind) iyref,imref,idref,ihref
    integer(i_kind) iazbin,irbin,ielbin
    integer(i_kind) nminref,nminthis
    integer(i_kind) num_radars,ireadsb,ireadmg,next
    integer(i_kind) nsuper,nsuperall
    integer(i_kind) nthisrad,nthisbins
    integer(i_kind) idups,idups0
    integer(i_kind) nradials_in,nradials_fail_angmax,nradials_fail_time,nradials_fail_elb
    integer(i_kind) nradials_in1,nradials_fail_angmax1,nradials_fail_time1,nradials_fail_elb1
    integer(i_kind) nobs_in,nobs_badvr,nobs_badsr,nobs_lrbin,nobs_hrbin,nrange_max,irad
    integer(i_kind) nobs_in1,nobs_badvr1,nobs_badsr1,nobs_lrbin1,nobs_hrbin1,nrange_max1
    integer(i_kind) num_radars_max,num_radars_min
    integer(i_kind):: histo_el
    integer(i_kind) idate5(5)
    integer(i_kind),allocatable,dimension(:) :: indx,icount
    integer(i_kind),allocatable,dimension(:,:) :: ibins,ibins2

    real(r_double) hdr(10),hdr2(12),rwnd(3,n_gates_max),rdisttest(n_gates_max)

    real(r_kind) delaz,delel,delr,t
    real(r_kind) ddiffmin,ddiffmin0,distfact,range
    real(r_kind) stn_lat,stn_lon,stn_hgt,stn_az,stn_el
    real(r_kind) timemax,timemin
    real(r_kind) timemax1,timemin1
    real(r_kind) this_stalat,this_stalon,this_stahgt
    real(r_kind) rlon0,clat0,slat0,rlonglob,rlatglob,clat1,caz0,saz0,cdlon,sdlon,caz1,saz1
    real(r_kind) this_stalatr,thisazimuthr,thistiltr
    real(r_kind) thisrange,thisazimuth,thistilt,thisvr,thisvr2
    real(r_kind) corrected_tilt
    real(r_kind) thiserr,thistime,thislat,thislon,corrected_azimuth
    real(r_kind) rad_per_meter,thishgt
    real(r_kind) rlonloc,rlatloc
    real(r_kind) a43,aactual,b,c,selev0,celev0,epsh,erad,h,ha
    real(r_kind) celev,selev,gamma
    real(r_kind) vrmax,vrmin,errmax,errmin
    real(r_kind) vrmaxall,vrminall,errmaxall,errminall
    real(r_kind) delazmmax,rdelaz,rdelr,rdelel
    real(r_kind) delazmmaxall
    real(r_kind) deltiltmaxall,deltiltmax
    real(r_kind) deltiltminall,deltiltmin
    real(r_kind) deldistmaxall,deldistmax
    real(r_kind) deldistminall,deldistmin
    real(r_kind),dimension(max_num_radars) :: stn_lat_table,stn_lon_table, &
                 stn_hgt_table,master_lat_table,master_lon_table,master_hgt_table
    real(r_kind),dimension(max_num_radars*npe) :: stn_lat_table_all,stn_lon_table_all, &
                 stn_hgt_table_all

    real(r_quad) thiscount
    real(r_quad),dimension(6) :: binsx
    real(r_quad),allocatable,dimension(:,:,:) :: bins,bins_work
!            bins(1,...) radial distance
!            bins(2,...) azimuth
!            bins(3,...) elev angle
!            bins(4,...) vr
!            bins(5,...) vr**2
!            bins(6,...) relative time
!            ibins(...) count

!       bins(  :,   :,   :)
!            var  rbin*azbin*elbin rad#
    character(8) chdr,chdr2,subset
    character(10) date
    character(20) infile
    character(4),dimension(max_num_radars):: stn_id_table,master_stn_table
    character(4) stn_id_table_all(max_num_radars*npe)
    character(4) stn_id,this_staid
    character(4*max_num_radars) cstn_id_table,cmaster_stn_table
    equivalence(stn_id_table(1),cstn_id_table)
    equivalence(master_stn_table(1),cmaster_stn_table)
    equivalence (chdr,hdr(1))
    equivalence (chdr2,hdr2(1))

    logical rite,print_verbose
    logical lradar

    print_verbose=.false.
    if(verbose) print_verbose=.true.

    ! define infile if using either option for radial winds.
    do i=1,ndat
       if(trim(dtype(i))=='rw'.and.trim(dsis(i))=='l2rw'.and.vadwnd_l2rw_qc)then
          infile=trim(dfile(i)) !l2rwbufr
          exit
       else if(trim(dtype(i))=='rw'.and.trim(dsis(i))=='l2rw'.and..not.vadwnd_l2rw_qc)then
          infile=trim(dfile(i)) !l2rwbufr with no VAD QC
          exit
       else
          infile=''
       end if
    end do
    if(mype==0) write(6,*)'RADAR_BUFR_READ_ALL:  radial wind infile for superobbing = ',infile
    if (infile=='') then
       if(mype==0) write(6,*) '***WARNING*** NOT USING ANY RADIAL WIND OBSERVATIONS!!!'
       return ! Go back to gsisub if nothing to do.
    end if

    rad_per_meter= one/rearth
    erad = rearth
    nazbin=nint(r360/del_azimuth)
    nrbin=nint(range_max/del_range)
    nelbin=nint(elev_angle_max/del_elev)
    delaz=r360/nazbin
    delr=range_max/nrbin
    delel=elev_angle_max/nelbin
    rdelaz=one/delaz
    rdelr =one/delr
    rdelel=one/delel

    num_radars=0
    do i=1,max_num_radars
       stn_id_table(i)='ZZZZ'
    end do
    stn_lat_table=r99999
    stn_lon_table=r99999
    stn_hgt_table=r99999

    rite = .false.
    if (mype==0) rite=.true.

!   Open bufr file with openbf to initialize bufr table, etc in bufrlib
    inbufr=10
    open(inbufr,file=infile,form='unformatted')
    read(inbufr,iostat=iret)subset
    if(iret/=0) then
       if(rite) write(6,*)'RADAR_BUFR_READ_ALL:  problem opening level 2 bufr file "l2rwbufr"'
       close(inbufr)
       return
    end if
    rewind inbufr
    lundx=inbufr
    call openbf(inbufr,'IN',lundx)
    call datelen(10)
    if(l2superob_only) then
       write(date,'( i10)') idate
       read (date,'(i4,3i2)') iyref,imref,idref,ihref
       if(rite) write(6,*)' create superobs only, radar file date = ',iyref,imref,idref,ihref
    else
       iyref=iadate(1)              !????  add mods so can be used in global mode
       imref=iadate(2)
       idref=iadate(3)
       ihref=iadate(4)
       if(rite) write(6,*)' using restart file date = ',iyref,imref,idref,ihref
    end if
    if(rite) write(6,*)'RADAR_BUFR_READ_ALL: analysis time is ',iyref,imref,idref,ihref
    idate5(1)=iyref
    idate5(2)=imref
    idate5(3)=idref
    idate5(4)=ihref
    idate5(5)=0          ! minutes
    call w3fs21(idate5,nminref)

!    Do an initial read of a bit of data to infer what multiplying factor is for
!    radial distance.  There is a possible ambiguity, where the scaling is either
!    125 or 250. If the minimum difference between gate distances is 2, then factor
!    is 125, if = 1 then factor is 250.

    idups=0
    nobs_in=0
    ddiffmin=huge(ddiffmin)
    next=0
    do while(ireadmg(inbufr,subset,idate)>=0)
       next=next+1
       if(next == npe)next=zero
       if(next /= mype)cycle
       read(subset,'(2x,i6)')isubset
       if(isubset>6033) then
          iret=6034
          exit
       end if
       do while (ireadsb(inbufr)==0)
          call ufbint(inbufr,rdisttest,1,n_gates_max,n_gates,'DIST125M')
          if(n_gates>1) then
             do i=1,n_gates-1
                if(nint(abs(rdisttest(i+1)-rdisttest(i)))==0) then
                   idups=idups+1
                else
                   ddiffmin=min(abs(rdisttest(i+1)-rdisttest(i)),ddiffmin)
                end if
             end do
          end if
          call ufbint(inbufr,hdr2,12,1,levs,'SSTN CLAT CLON HSMSL HSALG ANEL YEAR MNTH DAYS HOUR MINU SECO')
          if(hdr2(6)>elev_angle_max) cycle
          idate5(1)=nint(hdr2(7)) ; idate5(2)=nint(hdr2(8)) ; idate5(3)=nint(hdr2(9))
          idate5(4)=nint(hdr2(10)) ; idate5(5)=nint(hdr2(11))
          call w3fs21(idate5,nminthis)
          t=(real(nminthis-nminref,r_kind)+real(nint(hdr2(12)),r_kind)*rinv60)*rinv60
          if(abs(t)>del_time) cycle
          nobs_in=nobs_in+n_gates
          stn_id=chdr2
          ibyte=index(cstn_id_table,stn_id)
          if(ibyte==0) then
             num_radars=num_radars+1
             if(num_radars>max_num_radars) then
                write(6,*)'RADAR_BUFR_READ_ALL:  stop processing level 2 radar ',&
                     'bufr file--increase parameter max_num_radars'
                call stop2(99)
             end if
             stn_lat=hdr2(2)
             stn_lon=hdr2(3)
             stn_hgt=hdr2(4)+hdr2(5)
             stn_id_table(num_radars)=stn_id
             stn_lon_table(num_radars)=stn_lon
             stn_lat_table(num_radars)=stn_lat
             stn_hgt_table(num_radars)=stn_hgt
          end if
       end do
    end do

    call mpi_allreduce(num_radars,num_radars_max,1,mpi_integer4,mpi_max,mpi_comm_world,ierror)
    if(num_radars_max<=0) then
       if(rite) write(6,*)'RADAR_BUFR_READ_ALL:  NO RADARS KEPT IN radar_bufr_read_all, ',&
            'continue without level 2 data'
       call closbf(inbufr)
       return
    end if
    call mpi_reduce(num_radars,num_radars_min,1,mpi_integer4,mpi_min,0,mpi_comm_world,ierror)
    if(rite) write(6,*)' min,max num_radars=',num_radars_min,num_radars_max

!  Create master station list

!  First gather all stn id and lat,lon,hgt lists
    call mpi_allgather(stn_id_table,max_num_radars,mpi_integer4, &
         stn_id_table_all,max_num_radars,mpi_integer4,mpi_comm_world,ierror)
    call mpi_allgather(stn_lat_table,max_num_radars,mpi_real8, &
         stn_lat_table_all,max_num_radars,mpi_real8,mpi_comm_world,ierror)
    call mpi_allgather(stn_lon_table,max_num_radars,mpi_real8, &
         stn_lon_table_all,max_num_radars,mpi_real8,mpi_comm_world,ierror)
    call mpi_allgather(stn_hgt_table,max_num_radars,mpi_real8, &
         stn_hgt_table_all,max_num_radars,mpi_real8,mpi_comm_world,ierror)

!   Create unique master list of all radar names,lats,lons
    jj=0
    do j=1,max_num_radars*npe
       if(stn_id_table_all(j) /= 'ZZZZ')then
          jj=jj+1
          stn_id_table_all(jj)=stn_id_table_all(j)
          stn_lat_table_all(jj)=stn_lat_table_all(j)
          stn_lon_table_all(jj)=stn_lon_table_all(j)
          stn_hgt_table_all(jj)=stn_hgt_table_all(j)
       end if
    end do
    jjall=jj
    num_radars_0=0
    jj=1
    outer: do j=1,jjall
        num_radars_0=num_radars_0+1
        if(num_radars_0 > max_num_radars) then
	  write(6,*)'RADAR_BUFR_READ_ALL:  stop processing level 2 radar ',&
	     'bufr file--increase parameter max_num_radars'
          call stop2(99)
        end if
        master_stn_table(num_radars_0)=stn_id_table_all(jj)
        master_lat_table(num_radars_0)=stn_lat_table_all(jj)
        master_lon_table(num_radars_0)=stn_lon_table_all(jj)
        master_hgt_table(num_radars_0)=stn_hgt_table_all(jj)
        numzzzz=0
        do jjj=jj+1,jjall
          if(stn_id_table_all(jjj) /= 'ZZZZ')then
	     if(stn_id_table_all(jjj)==master_stn_table(num_radars_0)) then
	        stn_id_table_all(jjj)='ZZZZ'
	     else
	        if(numzzzz == 0)numzzzz=jjj
	     end if
	  end if
        end do
        if(numzzzz==0) exit outer
        jj=numzzzz
    end do outer
    allocate(indx(num_radars_0))
    call indexset(indx)
    call indexsort(indx,master_stn_table)

    call mpi_allreduce(ddiffmin,ddiffmin0,1,mpi_real8,mpi_min,mpi_comm_world,ierror)
    call mpi_allreduce(idups,idups0,1,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
    distfact=zero
    if(nint(ddiffmin0)==1)     distfact=r250
    if(nint(ddiffmin0)==2)     distfact=r125
    if(distfact==zero) then
       write(6,*)'RADAR_BUFR_READ_ALL:  problem with level 2 bufr file, ',&
	    'gate distance scale factor undetermined, going with 125'
       distfact=r125
    end if

    timemax=-huge(timemax)
    timemin=huge(timemin)
    nradials_in=0
    nradials_fail_angmax=0
    nradials_fail_time=0
    nradials_fail_elb=0
    nobs_badvr=0
    nobs_badsr=0
    nobs_lrbin=0
    nobs_hrbin=0
    nrange_max=0
    nthisrad=nrbin*nazbin*nelbin
    nthisbins=6*nthisrad

! reopen and reread the file for data this time

    call closbf(inbufr)
    open(inbufr,file=infile,form='unformatted')
    call openbf(inbufr,'IN',inbufr)

    allocate(bins(6,nthisrad,num_radars_0),ibins(nthisrad,num_radars_0))
    bins=zero_quad
    ibins=0

    next=0
    do while(ireadmg(inbufr,subset,idate)>=0)
       next=next+1
       if(next == npe)next=zero
       if(next /= mype)cycle
       read(subset,'(2x,i6)')isubset
       if(isubset>6033) then
	  iret=6034
	  exit
       end if
       do while (ireadsb(inbufr)==0)
	  call ufbint(inbufr,hdr,10,1,levs, &
	       'SSTN YEAR MNTH DAYS HOUR MINU SECO ANAZ ANEL QCRW')
	  nradials_in=nradials_in+1
	  if(hdr(9)>elev_angle_max) then
	     nradials_fail_angmax=nradials_fail_angmax+1
	     cycle
	  end if
	  idate5(1)=nint(hdr(2)) ; idate5(2)=nint(hdr(3)) ; idate5(3)=nint(hdr(4))
	  idate5(4)=nint(hdr(5)) ; idate5(5)=nint(hdr(6))
	  call w3fs21(idate5,nminthis)
	  t=(real(nminthis-nminref,r_kind)+real(nint(hdr(7)),r_kind)*rinv60)*rinv60
	  timemax=max(t,timemax)
	  timemin=min(t,timemin)
	  if(abs(t)>del_time) then
	     nradials_fail_time=nradials_fail_time+1
	     cycle
	  end if
	  stn_az=r90-hdr(8)
	  stn_el=hdr(9)
	  iazbin=stn_az*rdelaz
	  iazbin=mod(iazbin,nazbin)
	  if(iazbin<0) iazbin=iazbin+nazbin
	  iazbin=iazbin+1
	  if(iazbin<=0.or.iazbin>nazbin) then
	     write(6,*)'RADAR_BUFR_READ_ALL:  error in getting iazbin, program stops'
	     call stop2(99)
	  end if
	  ielbin=ceiling(stn_el*rdelel)
	  if(ielbin<1.or.ielbin>nelbin) then
	     nradials_fail_elb=nradials_fail_elb+1
	     cycle
	  end if
	  stn_id=chdr
	  ibyte=index(cmaster_stn_table,stn_id)
	  if(ibyte==0) then
	     write(6,*) ' index error in radar_bufr_read_all -- program stops -- ',ibyte,stn_id
	     call stop2(99)
	  else
	     krad=1+(ibyte-1)/4
	  end if

	  call ufbint(inbufr,rwnd,3,n_gates_max,n_gates,'DIST125M DMVR DVSW')
	  do i=1,n_gates
	     range=distfact*rwnd(1,i)
	     if(range>range_max) then
		nrange_max=nrange_max+1
		cycle
	     end if
	     if(rwnd(2,i)>r1e5_double) then
		nobs_badvr=nobs_badvr+1
		cycle
	     end if
	     if(rwnd(3,i)>r1e5_double) then
		nobs_badsr=nobs_badsr+1
		cycle
	     end if
	     irbin=ceiling(range*rdelr)
	     if(irbin<1) then
		nobs_lrbin=nobs_lrbin+1
		cycle
	     end if
	     if(irbin>nrbin) then
		nobs_hrbin=nobs_hrbin+1
		cycle
	     end if
	     iloc=nrbin*(nazbin*(ielbin-1)+(iazbin-1))+irbin
	     bins(1,iloc,krad)=bins(1,iloc,krad)+range
	     bins(2,iloc,krad)=bins(2,iloc,krad)+stn_az
	     bins(3,iloc,krad)=bins(3,iloc,krad)+stn_el
	     bins(4,iloc,krad)=bins(4,iloc,krad)+rwnd(2,i)
	     bins(5,iloc,krad)=bins(5,iloc,krad)+rwnd(2,i)**2
	     bins(6,iloc,krad)=bins(6,iloc,krad)+t
	     ibins(iloc,krad)=ibins(iloc,krad)+1
	  end do

       end do          !  end do while
    end do             !  loop over blocks
    call closbf(inbufr)

    allocate(ibins2(nthisrad,num_radars_0))
    call mpi_allreduce(ibins,ibins2,nthisrad*num_radars_0,mpi_integer4,mpi_sum,mpi_comm_world,ierror)
    deallocate(ibins)

    call mpi_reduce(nradials_in,nradials_in1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nradials_fail_angmax,nradials_fail_angmax1,1,&
	 mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nradials_fail_time,nradials_fail_time1,1,&
	 mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nradials_fail_elb,nradials_fail_elb1,1,&
	 mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nobs_in,nobs_in1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nobs_badvr,nobs_badvr1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nobs_badsr,nobs_badsr1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nobs_lrbin,nobs_lrbin1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nobs_hrbin,nobs_hrbin1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(nrange_max,nrange_max1,1,mpi_integer4,mpi_sum,0,mpi_comm_world,ierror)
    call mpi_reduce(timemax,timemax1,1,mpi_real8,mpi_max,0,mpi_comm_world,ierror)
    call mpi_reduce(timemin,timemin1,1,mpi_real8,mpi_min,0,mpi_comm_world,ierror)

    if(mype==0 ) then

       allocate(icount(num_radars_0))
       if(rite)write(6,*)'RADAR_BUFR_READ_ALL:  num_radars_0 = ',num_radars_0
       do irad=1,num_radars_0
          krad=indx(irad)
          icount(krad)=0
          do i=1,nthisrad
            if(ibins2(i,krad) >= minnum)icount(krad)=icount(krad)+1
          end do
	  if(rite)write(6,'(" master list radar ",i3," stn id,lat,lon,hgt,num = ",a4,2f10.2,f8.1,i10)') &
	       irad,master_stn_table(krad),master_lat_table(krad),master_lon_table(krad), &
               master_hgt_table(krad),icount(krad)
       end do
       if(rite)then
          write(6,*)'RADAR_BUFR_READ_ALL:  ddiffmin,distfact,idups=',ddiffmin0,distfact,idups0
          write(6,*)' nthisrad=',nthisrad
          write(6,*)' nthisbins=',nthisbins
          write(6,*)' timemin,max=',timemin1,timemax1
          write(6,*)' nradials_in=',nradials_in1
          write(6,*)' nradials_fail_angmax=',nradials_fail_angmax1
          write(6,*)' nradials_fail_time=',nradials_fail_time1
          write(6,*)' nradials_fail_elb=',nradials_fail_elb1
          write(6,*)' nobs_in=',nobs_in1
          write(6,*)' nobs_badvr=',nobs_badvr1
          write(6,*)' nobs_badsr=',nobs_badsr1
          write(6,*)' nobs_lrbin=',nobs_lrbin1
          write(6,*)' nobs_hrbin=',nobs_hrbin1
          write(6,*)' nrange_max=',nrange_max1
       end if

!    Print out histogram of counts by ielbin to see where angles are
       do ielbin=1,nelbin
          histo_el=0
          do krad=1,num_radars_0
	     do iazbin=1,nazbin
		do irbin=1,nrbin
		   iloc=nrbin*(nazbin*(ielbin-1)+(iazbin-1))+irbin
		   histo_el=histo_el+ibins2(iloc,krad)
		end do
	     end do
	  end do
	  if(rite)write(6,'(" ielbin,histo_el=",i6,i20)')ielbin,histo_el
       end do

!   Prepare to create superobs and write out.
       open(inbufr,file='radar_supobs_from_level2',form='unformatted',iostat=iret)
       rewind inbufr
       nsuperall=0
       vrmaxall=-huge(vrmaxall)
       vrminall=huge(vrminall)
       errmaxall=-huge(errmaxall)
       errminall=huge(errminall)
       delazmmaxall=-huge(delazmmaxall)
       deltiltmaxall=-huge(deltiltmaxall)
       deldistmaxall=-huge(deldistmaxall)
       deltiltminall=huge(deltiltminall)
       deldistminall=huge(deldistminall)
    end if
    allocate(bins_work(6,nthisrad,npe))
    do irad=1,num_radars_0
       krad=indx(irad)
       if(r_double==r_quad) then
          call mpi_gather(bins(1,1,krad),nthisbins,mpi_real8 ,bins_work,nthisbins, &
                       mpi_real8 ,0,mpi_comm_world,ierror)
       else
          call mpi_gather(bins(1,1,krad),nthisbins,mpi_real16,bins_work,nthisbins, &
                       mpi_real16,0,mpi_comm_world,ierror)
       endif
       if(mype == 0)then

!   Create superobs and write out.
	  nsuper=0
	  vrmax=-huge(vrmax)
	  vrmin=huge(vrmin)
	  errmax=-huge(errmax)
	  errmin=huge(errmin)
	  delazmmax=-huge(delazmmax)
	  deltiltmax=-huge(deltiltmax)
	  deldistmax=-huge(deldistmax)
	  deltiltmin=huge(deltiltmin)
	  deldistmin=huge(deldistmin)
	  this_stalat=master_lat_table(krad)
	  if(abs(this_stalat)>r89_5) cycle
	  this_stalon=master_lon_table(krad)
	  rlon0=deg2rad*this_stalon
	  this_stalatr=this_stalat*deg2rad
	  clat0=cos(this_stalatr) ; slat0=sin(this_stalatr)
	  this_staid=master_stn_table(krad)
	  this_stahgt=master_hgt_table(krad)
	  do iii=1,nthisrad
! If single radar exp, set lradar false if not the specified radar.
           lradar=.true.
           if(lsingleradar) then
             if(this_staid /= singleradar) lradar=.false.
           end if

           if(lradar) then ! Logical for when running single radar exp.
	     if(ibins2(iii,krad) < minnum) cycle

	     thiscount=one_quad/real(ibins2(iii,krad),r_quad)
             do i=1,6
               binsx(i)=bins_work(i,iii,1)
             end do
             do k=2,npe
               do i=1,6
	         binsx(i)=binsx(i)+bins_work(i,iii,k)
               end do
             end do
             do i=1,6
               binsx(i)=binsx(i)*thiscount
             end do
	     thisrange=  binsx(1)
	     thisazimuth=binsx(2)
	     thistilt=binsx(3)
	     thisvr=binsx(4)
	     vrmax=max(vrmax,thisvr)
	     vrmin=min(vrmin,thisvr)
	     thisvr2=binsx(5)
	     thiserr=sqrt(abs(thisvr2-thisvr**2))
	     errmax=max(errmax,thiserr)
	     errmin=min(errmin,thiserr)
	     thistime=binsx(6)

!            Compute obs height here
!            Use 4/3rds rule to get elevation of radar beam
!            (if local temperature, moisture available, then vertical position
!             might be estimated with greater accuracy by ray tracing )

	     aactual=erad+this_stahgt
	     a43=four_thirds*aactual
	     thistiltr=thistilt*deg2rad
	     selev0=sin(thistiltr)
	     celev0=cos(thistiltr)
	     b=thisrange*(thisrange+two*aactual*selev0)
	     c=sqrt(aactual*aactual+b)
	     ha=b/(aactual+c)
	     epsh=(thisrange*thisrange-ha*ha)/(r8*aactual)
	     h=ha-epsh
	     thishgt=this_stahgt+h

!            Get corrected tilt angle
	     celev=celev0
	     selev=selev0
	     if(thisrange>=one) then
		celev=a43*celev0/(a43+h)
		selev=(thisrange*thisrange+h*h+two*a43*h)/(two*thisrange*(a43+h))
	     end if
	     corrected_tilt=atan2(selev,celev)*rad2deg
	     deltiltmax=max(corrected_tilt-thistilt,deltiltmax)
	     deltiltmin=min(corrected_tilt-thistilt,deltiltmin)
	     gamma=half*thisrange*(celev0+celev)
	     deldistmax=max(gamma-thisrange,deldistmax)
	     deldistmin=min(gamma-thisrange,deldistmin)

!            Get earth lat lon of superob
	     thisazimuthr=thisazimuth*deg2rad
	     rlonloc=rad_per_meter*gamma*cos(thisazimuthr)
             rlatloc=rad_per_meter*gamma*sin(thisazimuthr)
             call invtllv(rlonloc,rlatloc,rlon0,clat0,slat0,rlonglob,rlatglob)
             thislat=rlatglob*rad2deg
             thislon=rlonglob*rad2deg

!            Keep away from poles, rather than properly deal with polar singularity
	     if(abs(thislat)>r89_5) cycle

!            Get corrected azimuth
             clat1=cos(rlatglob)
             caz0=cos(thisazimuthr)
             saz0=sin(thisazimuthr)
             cdlon=cos(rlonglob-rlon0)
             sdlon=sin(rlonglob-rlon0)
             caz1=clat0*caz0/clat1
             saz1=saz0*cdlon-caz0*sdlon*slat0
             corrected_azimuth=atan2(saz1,caz1)*rad2deg
             delazmmax=max(min(abs(corrected_azimuth-thisazimuth-r720),&
                  abs(corrected_azimuth-thisazimuth-r360),&
                  abs(corrected_azimuth-thisazimuth     ),&
                  abs(corrected_azimuth-thisazimuth+r360),&
                  abs(corrected_azimuth-thisazimuth+r720)),delazmmax)

             write(inbufr) this_staid,this_stalat,this_stalon,this_stahgt, &
                  thistime,thislat,thislon,thishgt,thisvr,corrected_azimuth,&
                  thiserr,corrected_tilt,gamma
             nsuper=nsuper+1
           end if
          end do
          if(nsuper > 0)then
            write(6,*)' for radar ',this_staid,' nsuper=',nsuper,' delazmmax=',delazmmax
            write(6,*)' vrmin,max=',vrmin,vrmax,' errmin,max=',errmin,errmax
            write(6,*)' deltiltmin,max=',deltiltmin,deltiltmax,' deldistmin,max=',deldistmin,deldistmax
            vrminall=min(vrminall,vrmin)
            vrmaxall=max(vrmaxall,vrmax)
            errminall=min(errminall,errmin)
            errmaxall=max(errmaxall,errmax)
            delazmmaxall=max(delazmmaxall,delazmmax)
            deltiltmaxall=max(deltiltmaxall,deltiltmax)
            deldistmaxall=max(deldistmaxall,deldistmax)
            deltiltminall=min(deltiltminall,deltiltmin)
            deldistminall=min(deldistminall,deldistmin)
            nsuperall=nsuperall+nsuper
          end if
       end if
    end do
    if(mype == 0)then
       if(rite)then
          write(6,*)' total number of superobs written=',nsuperall
          write(6,*)'  vrmin,maxall=',vrminall,vrmaxall
          write(6,*)' errmin,maxall=',errminall,errmaxall
          write(6,*)' delazmmaxall=',delazmmaxall
          write(6,*)' deltiltmin,maxall=',deltiltminall,deltiltmaxall
          write(6,*)' deldistmin,maxall=',deldistminall,deldistmaxall
       end if
       close(inbufr)
    end if
    deallocate(bins_work,bins,ibins2)
    if(l2superob_only) then
       call mpi_finalize(ierror)
       stop
    end if

end subroutine radar_bufr_read_all

SUBROUTINE tllv(ALM,APH,TLMO,CTPH0,STPH0,TLM,TPH)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    tllv
!
!   prgrmmr:
!
! abstract:
!
! program history log:
!   2008-04-21  safford -- add subprogram doc block, rm unused uses
!
!   input argument list:
!     alm   -- input earth longitude
!     aph   -- input earth latitude
!     tlmo  -- input earth longitude of rotated grid origin (radrees)
!     ctph0 -- cos(earth lat of rotated grid origin)
!     stph0 -- sin(earth lat of rotated grid origin)
!
!   output argument list:
!     tlm   -- rotated grid longitude
!     tph   -- rotated grid latitude
!
! attributes:
!   language:  f90
!   machine:
!
!$$$ end documentation block

  use kinds, only:  r_kind
  implicit none

  real(r_kind),intent(in   ) :: alm,aph,tlmo,ctph0,stph0
  real(r_kind),intent(  out) :: tlm,tph

  real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom

  RELM=ALM-TLMO
  SRLM=SIN(RELM)
  CRLM=COS(RELM)
  SPH=SIN(APH)
  CPH=COS(APH)
  CC=CPH*CRLM
  ANUM=CPH*SRLM
  DENOM=CTPH0*CC+STPH0*SPH
  TLM=ATAN2(ANUM,DENOM)
  TPH=ASIN(CTPH0*SPH-STPH0*CC)

END SUBROUTINE tllv

SUBROUTINE invtllv(ALM,APH,TLMO,CTPH0,STPH0,TLM,TPH)
!$$$  subprogram documentation block
!                .      .    .
! subprogram:    invtllv
!
!   prgrmmr:
!
! abstract:  inverse of tllv:  input ALM,APH is rotated lon,lat
!                   output is earth lon,lat, TLM,TPH
!
! program history log:
!   2008-03-25  safford -- add subprogram doc block, rm unused uses
!
!   input argument list:
!     alm   -- input earth longitude
!     aph   -- input earth latitude
!     tlmo  -- input earth longitude of rotated grid origin (radrees)
!     ctph0 -- cos(earth lat of rotated grid origin)
!     stph0 -- sin(earth lat of rotated grid origin)
!
!   output argument list:
!     tlm   -- rotated grid longitude
!     tph   -- rotated grid latitude
!
! attributes:
!   language:  f90
!   machine:
!
!$$$ end documentation block

  use kinds, only:  r_kind
  implicit none

  real(r_kind),intent(in   ) :: alm,aph,tlmo,ctph0,stph0
  real(r_kind),intent(  out) :: tlm,tph

  real(r_kind):: relm,srlm,crlm,sph,cph,cc,anum,denom

  RELM=ALM
  SRLM=SIN(RELM)
  CRLM=COS(RELM)
  SPH=SIN(APH)
  CPH=COS(APH)
  CC=CPH*CRLM
  ANUM=CPH*SRLM
  DENOM=CTPH0*CC-STPH0*SPH
  TLM=tlmo+ATAN2(ANUM,DENOM)
  TPH=ASIN(CTPH0*SPH+STPH0*CC)

END SUBROUTINE invtllv
end module read_l2bufr_mod


text:
cxx      SUBROUTINE BSUBSTF( ZS,N,IMODEL,LAG,K,IL,LG1,LG2,F,CNST,ZMEAN,SUM,
      SUBROUTINE BSUBSTF( ZS,N,IMODEL,LAG,K,IL,LG1,LG2,ZMEAN,SUM,
     *M,AICM,SDM,A1,SD,AIC,DIC,AICB,SDB,EK,A2,IND,C,C1,C2,B,OEIC,ESUM,
     *OMEAN,OM,E,EMEAN,VARI,SKEW,PEAK,COV,SXX )
C

C
cc      PROGRAM BSUBST
C.......................................................................
C.....PLANNED BY H.AKAIKE...............................................
C.....DESIGNED BY H.AKAIKE AND G.KITAGAWA...............................
C.....PROGRAMMED BY G.KITAGAWA AND F.TADA...............................
C.....ADDRESS: THE INSTITUTE OF STATISTICAL MATHEMATICS, 4-6-7 MINAMI-AZ
C..............MINATO-KU, TOKYO 106, JAPAN..............................
C.....DATE OF THE LATEST REVISION:  MAY  14, 1979.......................
C.......................................................................
C.....THIS PROGRAM WAS ORIGINALLY PUBLISHED IN "TIMSAC-78", BY H.AKAIKE,
C.....G.KITAGAWA, E.ARAHATA AND F.TADA, COMPUTER SCIENCE MONOGRAPHS, NO.
C.....THE INSTITUTE OF STATISTICAL MATHEMATICS, TOKYO, 1979.............
C.......................................................................
C     TIMSAC 78.1.3.
C     _                 ____ _
C     BAYESIAN TYPE ALL SUBSET ANALYSIS OF TIME SERIES BY A MODEL LINEAR
C     IN PARAMETERS
C
C     THIS PROGRAM PRODUCES BAYESIAN ESTIMATES OF TIME SERIES MODELS SUC
C     PURE AR MODELS, AR-MODELS WITH NON-LINEAR TERMS, AR-MODELS WITH PO
C     TYPE MEAN VALUE FUNCTIONS, ETC.  THE GOODNESS OF FIT OF A MODEL IS
C     CHECKED BY THE ANALYSIS OF SEVERAL STEPS AHEAD PREDICTION ERRORS.
C     BY PREPARING AN EXTERNAL SUBROUTINE SETX PROPERLY, ANY TIME SERIES
C     WHICH IS LINEAR IN PARAMETERS CAN BE TREATED.
C       ----------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS PROGRAM:
C             REDATA
C             REDLAG
C             SETLAG
C             REDREG
C             REDUCT
C             ARMFIT
C             SBBAYS
C             CHECK
C       ----------------------------------------------------------------
C
C       INPUTS REQUIRED:
C          MT:    ORIGINAL DATA INPUT DEVICE SPECIFICATION
C          IMODEL:=1  AUTOREGRESSIVE MODEL
C                 =2  POLYNOMIAL TYPE NON-LINEAR MODEL (LAG'S READ IN )
C                 =3  POLYNOMIAL TYPE NON-LINEAR MODEL (LAG'S AUTOMATICA
C                 =4  AR-MODEL WITH POLYNOMIAL MEAN VALUE FUNCTION
C                 =5  ANY NON-LINEAR MODEL
C                 =6  POLYNOMIAL TYPE EXPONENTIALLY DAMPED NON-LINEAR MO
C                 =7  THIS MODEL IS RESERVED FOR THE USER'S OPTIONAL USE
C          LAG:   MAXIMUM TIME LAG USED IN THE MODEL
C          K:     NUMBER OF REGRESSORS
C          IL:    PREDICTION ERRORS CHECKING (UP TO IL-STEPS AHEAD) IS R
C                 N*IL SHOULD BE LESS THAN OR EQUAL TO 20000
C
C       --   THE FOLLOWING INPUTS ARE REQUIRED AT SUBROUTINE REDATA   --
C
C          TITLE:   ORIGINAL DATA SPECIFICATION
C          N:       DATA LENGTH
C          DFORM:   INPUT DATA FORMAT SPECIFICATION STATEMENT
C                   -- EXAMPLE --  (8F10.5 )
C          X(I) (I=1,N):   ORIGINAL DATA
C       ----------------------------------------------------------------
C
cc      !DEC$ ATTRIBUTES DLLEXPORT :: BSUBSTF
C
      PARAMETER  ( MJ2 = 101 )
C
cxx      IMPLICIT  REAL * 8  ( A-H , O-Z )
CC      REAL*4     Z(10000) , TITLE(20) , TTL(5) , E(20000)
cc      REAL*4     TITLE(20) , TTL(5)
cc      DIMENSION  Z(10000), E(20000)
cc      DIMENSION  X(200,101), D(200) , A(100) , B(100)
cc      DATA  TTL  / 4H   B,4HAYES,4HIAN ,4HMODE,4HL    /
cxx      COMMON     / BBB / L1(50) , L2(50) , L3(50) , SD0 , CONST
cc      REAL * 4   CONST , SD0
cxx      DIMENSION  ZS(N), Z(N), E(N,IL)
cxx      DIMENSION  X(N,K+1), A1(K), A2(K), B(K)
cxx      DIMENSION  LG1(3,K), LG2(5)
cxx      DIMENSION  SD(K+1), AIC(K+1), DIC(K+1)
cxx      DIMENSION  IND(K), C(K), C1(K+1), C2(K), ESUM(K+1)
cxx      DIMENSION  EMEAN(IL), VARI(IL), SKEW(IL), PEAK(IL), COV(MJ2)
cxx      DIMENSION  SXX(121), FA(N-LAG,IL)
cc      CHARACTER*4 F(20,K+1)
cxx      INTEGER*1  F(80,K+1)
cxx      CHARACTER(4) G
      INTEGER :: N, IMODEL, LAG, K, IL, LG1(3,K), LG2(5), M, IND(K)
      REAL(8) :: ZS(N), ZMEAN, SUM, AICM, SDM, A1(K), SD(K+1), AIC(K+1),
     1           DIC(K+1), AICB, SDB, EK, A2(K), C(K), C1(K+1), C2(K),
     2           B(K), OEIC, ESUM(K+1), OMEAN, OM, E(N,IL),
     3           EMEAN(IL), VARI(IL), SKEW(IL), PEAK(IL), COV(MJ2),
     4           SXX(121)
      REAL(8) :: Z(N), X(N,K+1), FA(N-LAG,IL), SD0, CONST
      CHARACTER(4) G
C
      COMMON     / BBB / L1(50) , L2(50) , L3(50) , SD0 , CONST
      COMMON     / EEE /  G(20,31)
      COMMON     / AAA /  NN
C
C        EXTERNAL SUBROUTINE DECLARATION:
C
      EXTERNAL  SETX1
      EXTERNAL  SETX2
      EXTERNAL  SETX4
cx      EXTERNAL  SETX5
cx      EXTERNAL  SETX6
cc      EXTERNAL  SETX7
      EXTERNAL  PRDCT1
      EXTERNAL  PRDCT2
cx      EXTERNAL  PRDCT3
cx      EXTERNAL  PRDCT6
C
cc	CHARACTER(100) IFLNAM,OFLNAM
cc	CALL FLNAM2( IFLNAM,OFLNAM,NFL )
cc	IF ( NFL.EQ.0 ) GO TO 999
cc	IF ( NFL.EQ.2 ) THEN
cc	   OPEN( 6,FILE=OFLNAM,ERR=900,IOSTAT=IVAR)
cc	ELSE
cc	   CALL SETWND
cc	END IF
C
C
C        PARAMETERS:
C             MJ1:  ABSOLUTE DIMENSION FOR SUBROUTINE CALL
C
cc      IF ((IMODEL.LE.0) .OR. (IMODEL.GE.8)) GO TO 150
      IF ((IMODEL.LE.0) .OR. (IMODEL.GE.7)) GO TO 150
cc      MJ = 1000
cc      MJ1 = 200
      NN = N
      MJ = N
      MJ1 = N
      ISW = 1
      IPR = 2
C
CC      READ( 5,1 )     MT
cc      MT = 5
cc      OPEN( MT,FILE=IFLNAM,ERR=910,IOSTAT=IVAR,STATUS='OLD' )
cc      READ( 5,1 )     IMODEL , LAG , K , IL
cc      WRITE( 6,3 )
cc      IF( IMODEL.EQ.1 )  WRITE( 6,4 )   IMODEL
cc      IF( IMODEL.EQ.2 )  WRITE( 6,5 )   IMODEL
cc      IF( IMODEL.EQ.3 )  WRITE( 6,5 )   IMODEL
cc      IF( IMODEL.EQ.5 )  WRITE( 6,5 )   IMODEL
cc      IF( IMODEL.EQ.6 )  WRITE( 6,11 )   IMODEL
cc      WRITE( 6,6 )
cc      IF( IMODEL.EQ.1 )  WRITE( 6,7 )  K
cc      IF( IMODEL.EQ.2 )  WRITE( 6,8 )   LAG , K
cc      WRITE( 6,2 )     MT
C
C
C          ---------------------------------------
C          ORIGINAL DATA LOADING AND MEAN DELETION
C          ---------------------------------------
C
cc      CALL  REDATA( Z,N,MT,TITLE )
      CALL  REDATA( ZS,Z,N,ZMEAN,SUM )
      NMK = N - LAG
      LAG1 = LAG + 1
C
C          ---------------------
C          HOUSEHOLDER REDUCTION
C          ---------------------
C
cc      GO TO ( 10,20,30,40,50,60,70 ), IMODEL
cx      GO TO ( 10,20,30,40,50,60 ), IMODEL
cxx      GO TO ( 10,20,30,40 ), IMODEL
      IF ( IMODEL .EQ. 1 ) GO TO 10
      IF ( IMODEL .EQ. 2 ) GO TO 20
      IF ( IMODEL .EQ. 3 ) GO TO 30
      IF ( IMODEL .EQ. 4 ) GO TO 40
C
   10 K = LAG
cc      CALL  REDUCT( SETX1,Z,D,NMK,0,K,MJ1,LAG,X )
      CALL  REDUCT( SETX1,Z,NMK,0,K,MJ1,LAG,X )
      GO TO 100
C
cc   20 CALL  REDLAG( K )
   20 CONTINUE
      DO 21 I=1,K
         L1(I) = LG1(1,I)
         L2(I) = LG1(2,I)
         L3(I) = LG1(3,I)
   21 CONTINUE
cc      CALL  REDUCT( SETX2,Z,D,NMK,0,K,MJ1,LAG,X )
      CALL  REDUCT( SETX2,Z,NMK,0,K,MJ1,LAG,X )
      GO TO 100
C
cc   30 CALL  SETLAG( K )
cc      CALL  REDUCT( SETX2,Z,D,NMK,0,K,MJ1,LAG,X )
cxx   30 CALL  SETLAG( K,LG2(1),LG2(2),LG2(3),LG2(4),LG2(5) )
   30 CALL  SETLAG( KK,LG2(1),LG2(2),LG2(3),LG2(4),LG2(5) )
      DO 31 I=1,K
         LG1(1,I) = L1(I)
         LG1(2,I) = L2(I)
         LG1(3,I) = L3(I)
   31 CONTINUE
      CALL  REDUCT( SETX2,Z,NMK,0,K,MJ1,LAG,X )
      GO TO 100
C
cc   40 CALL  REDUCT( SETX4,Z,D,NMK,0,K,MJ1,LAG,X )
   40 CALL  REDUCT( SETX4,Z,NMK,0,K,MJ1,LAG,X )
      GO TO 100
C
cc   50 CALL  REDREG( K )
cx   50 CONTINUE
cx      DO 51 J=1,K+1
cx         DO 51 I=1,20
cx            II = (I-1)*4+1
cx            G(I,J) = CHAR(F(II,J)) // CHAR(F(II+1,J))
cx     *                             // CHAR(F(II+2,J)) //CHAR(F(II+3,J))
cx   51 CONTINUE
cc      CALL  REDUCT( SETX5,Z,D,NMK,0,K,MJ1,LAG,X )
cx      CALL  REDUCT( SETX5,Z,NMK,0,K,MJ1,LAG,X )
cx      GO TO 100
C
cx   60 CONTINUE
cx	SD0 = SUM
cx      CONST = CNST
cc      SD0 = 0.D0
cc      DO 65  I=1,N
cc   65 SD0 = SD0 + Z(I)*Z(I)
cc      SD0 = SD0 / N
cc      READ( 5,9 )     CONST
cc      WRITE( 6,12 )   SD0 , CONST
cc      CALL  REDLAG( K )
cx      DO 66 I=1,K
cx         L1(I) = LG1(1,I)
cx         L2(I) = LG1(2,I)
cx         L3(I) = LG1(3,I)
cx   66 CONTINUE
cc      CALL  REDUCT( SETX6,Z,D,NMK,0,K,MJ1,LAG,X )
cx      CALL  REDUCT( SETX6,Z,NMK,0,K,MJ1,LAG,X )
cx      GO TO 100
C
cc   70 CALL  REDUCT( SETX7,Z,D,NMK,0,K,MJ1,LAG,X )
cc   70 CALL  REDUCT( SETX7,Z,NMK,0,K,MJ1,LAG,X )
C
  100 CONTINUE
cc	CLOSE( MT )
C
C          ---------------
C          MAICE PROCEDURE
C          ---------------
C
cc      CALL  ARMFIT( X,K,LAG,NMK,ISW,TITLE,MJ1,A,SD,M )
cx      IFG=0
cx      CALL ARMFIT( X,K,LAG,NMK,ISW,MJ1,A1,M,SD,AIC,DIC,SDM,AICM,
cx     *             IFG,LU )
      CALL ARMFIT( X,K,LAG,NMK,ISW,MJ1,A1,M,SD,AIC,DIC,SDM,AICM)
C
C          ------------------
C          BAYESIAN PROCEDURE
C          ------------------
C
cc      CALL  SBBAYS( X,D,K,NMK,IPR,MJ1,A,SD )
      CALL  SBBAYS( X,K,NMK,IPR,MJ1,A2,SDB,EK,AICB,IND,C,C1,C2,B,
     *OEIC,ESUM,OMEAN,OM  )
              if( k.eq.32 ) return
C
cc      IF( IMODEL .EQ. 1 )  CALL  NRASPE( SD,A,B,K,0,121,TITLE )
      IF( IMODEL .EQ. 1 )  CALL  NRASPE( SDB,A2,B,K,0,120,SXX )
C
      NPS = LAG+1
C
C          -------------------------
C          PREDICTION ERROR CHECKING
C          -------------------------
C
cx      GO TO ( 110,120,120,150,130,140 ), IMODEL
cxx      GO TO ( 110,120,120,150 ), IMODEL
      IF ( IMODEL .EQ. 1 ) GO TO 110
      IF ( IMODEL .EQ. 2 ) GO TO 120
      IF ( IMODEL .EQ. 3 ) GO TO 120
      GO TO 150
C
cc  110 CALL CHECK( PRDCT1,Z,A,K,0,IL,NPS,N,0,MJ,E )
cxx  110 CALL CHECK( PRDCT1,Z,A2,K,0,IL,NPS,N,0,MJ,E,FA,EMEAN,VARI,SKEW,
  110 CALL CHECK( PRDCT1,Z,A2,K,0,IL,NPS,N,MJ,E,FA,EMEAN,VARI,SKEW,
     *            PEAK,COV,MJ2 )
      GO TO 150
C
cc  120 CALL  CHECK( PRDCT2,Z,A,K,0,IL,NPS,N,0,MJ,E )
cxx  120 CALL  CHECK( PRDCT2,Z,A2,K,0,IL,NPS,N,0,MJ,E,FA,EMEAN,VARI,SKEW,
  120 CALL  CHECK( PRDCT2,Z,A2,K,0,IL,NPS,N,MJ,E,FA,EMEAN,VARI,SKEW,
     *             PEAK,COV,MJ2 )
      GO TO 150
C
cc  130 CALL CHECK( PRDCT3,Z,A,K,0,IL,NPS,N,0,MJ,E )
cx  130 CALL CHECK( PRDCT3,Z,A2,K,0,IL,NPS,N,0,MJ,E,FA,EMEAN,VARI,SKEW,
cx     *            PEAK,COV,MJ2 )
C
cc  140 CALL  CHECK( PRDCT6,Z,A,K,0,IL,NPS,N,0,MJ,E )
cx  140 CALL  CHECK( PRDCT6,Z,A2,K,0,IL,NPS,N,0,MJ,E,FA,EMEAN,VARI,SKEW,
cx     *             PEAK,COV,MJ2 )
C
  150 CONTINUE
      RETURN
C
cxx    1 FORMAT( 16I5 )
cxx    2 FORMAT( /1H ,'ORIGINAL DATA INPUT DEVICE  MT =',I4 )
cxx    3 FORMAT( ' PROGRAM TIMSAC 78.1.3',/,'   SCALAR TIME SERIES MODEL FI
cxx     *TTING;     BAYESIAN PROCEDURE ( ALL SUBSET REGRESSIN TYPE )' )
cxx    4 FORMAT( //1H ,'MODEL TYPE',I2,/,'   < AUTOREGRESSIVE MODEL >',/,
cxx     11H ,10X,'Z(I) = A(1)*Z(I-1) + A(2)*Z(I-2) + ... + A(M)*Z(I-M) + E(
cxx     2I)' )
cxx    5 FORMAT( //1H ,'MODEL TYPE',I2,/,'   < NON-LINEAR MODEL >',/,1H ,10
cxx     1X,'Z(I) = A(1)*Y(I,1) + A(2)*Y(I,2) + ... + A(K)*Y(I,K) + E(I)' )
cxx    6 FORMAT( 1H ,2X,'WHERE',/,11X,'M:     ORDER OF THE MODEL',/,11X,'E(
cxx     1I):  GAUSSIAN WHITE NOISE WITH MEAN 0  AND  VARIANCE SD(M).' )
cxx    7 FORMAT( 1H ,'FITTING UP TO THE ORDER',I3,2X,'IS TRIED' )
cxx    8 FORMAT( 1H ,'MAXIMUM LAG =',I4,/,' NUMBER OF REGRESSORS =',I4 )
cxx    9 FORMAT( F10.0 )
cxx   11 FORMAT( //1H ,'MODEL TYPE',I2,/,'   < EXPONENTIALLY DAMPED NON-LIN
cxx     *EAR MODEL >',/,1H ,10X,'Z(I) = A(1)*Y(I,1) + A(2)*Y(I,2) + ... + A
cxx     *(K)*Y(I,K) + E(I)' )
cxx   12 FORMAT( 1H ,'SIGMA2 =',D15.5,5X,'CONST =',D15.5 )
C
      E N D
cc      REAL FUNCTION  BICOEF * 8( K,J )
cxx      REAL*8 FUNCTION  BICOEF( K,J )
      DOUBLE PRECISION FUNCTION BICOEF( K,J )
C
C     THIS FUNCTION RETURNS BINOMIAL COEFFICIENTS
C
C          F(K,J) = K]/(J]*(K-J)])
C
C       INPUTS:
C          K:     NUMBER OF OBJECTS
C          J:     NUMBER OF OBJECTS TAKEN
C
C       OUTPUT:
C          F:     NUMBER OF COMBINATIONS OF SELECTING J OBJECTS FROM
C                 A SET OF K OBJECTS
C
cxx      IMPLICIT REAL * 8 ( A-H , O-Z )
      REAL(8) :: SUM, DI
C
      KMJ = K-J
      SUM = 0.D0
      DO 10   I=1,K
      DI = I
cxx   10 SUM = SUM + DLOG( DI )
      SUM = SUM + DLOG( DI )
   10 CONTINUE
C
      IF( J .EQ. 0 )   GO TO 30
      DO 20   I=1,J
      DI = I
cxx   20 SUM = SUM - DLOG( DI )
      SUM = SUM - DLOG( DI )
   20 CONTINUE
C
   30 IF( KMJ .EQ. 0 )   GO TO 50
      DO 40   I=1,KMJ
      DI = I
cxx   40 SUM = SUM - DLOG( DI )
      SUM = SUM - DLOG( DI )
   40 CONTINUE
C
   50 BICOEF = DEXP( SUM )
      RETURN
C
      END
cc      SUBROUTINE  CHECK( PRDCT,X,A,K,L,IL,NPS,NPE,IPR,MJ,E )
cxx      SUBROUTINE  CHECK( PRDCT,X,A,K,L,IL,NPS,NPE,IPR,MJ,E,F,EMEAN,VARI,
      SUBROUTINE  CHECK( PRDCT,X,A,K,L,IL,NPS,NPE,MJ,E,F,EMEAN,VARI,
     *                   SKEW,PEAK,COV,MJ2 )
C
C     THIS SUBROUTINE DRAWS HISTGRAMS AND AUTOCOVARIANCE FUNCTION OF ORI
C     DATA OR PREDICTION ERRORS.
C
C       ----------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS SUBROUTINE
C             GRAPH1
C             GRAPH2
C             MOMENT
C             (PRDCT)
C       ----------------------------------------------------------------
C       INPUTS:
C          PRDCT:  EXTERNAL SUBROUTINE DESIGNATION
C          X:      ORIGINAL DATA
C          A:      REGRESSION COEFFICIENTS
C          K:      NUMBER OF REGRESSORS
C          L:      MA-ORDER ( THIS ARGUMENT IS ONLY USED FOR THE CHECKIN
C                  OF AR-MA MODEL)
C          IL:     MAXIMUM SPAN OF LONG RANGE PREDICTION
C                  =0     ANALYSIS OF ORIGINAL DATA
C                  >0     ANALYSIS OF MULTI-STEP (UP TO IL) PREDICTION E
C          NPS:    PREDICTION STARTING POSITION
C          NPE:    PREDICTION ENDING POSITION
C          IPR:    =0  MATRIX OF SEVERAL STEP AHEAD PREDICTION ERRORS SU
C                  =1  MATRIX OF SEVERAL STEP AHEAD PREDICTION ERRORS IS
C                      OUT
C          MJ:     ABSOLUTE DIMENSION OF E IN THE MAIN PROGRAM
C
C       OUTPUT:
C          E:      SEVERAL-STEPS PREDICTION ERRORS
C
C
cxx      IMPLICIT  REAL * 8  ( A-H , O-Z )
CC      REAL * 4  X , F , E
cx      DIMENSION  X(1) , E(MJ,1)
cc      DIMENSION  A(1) , COV(120)
cxx      DIMENSION  X(NPE) , E(MJ,IL)
cxx      DIMENSION  A(K) , COV(MJ2)
cxx      DIMENSION  EMEAN(IL), VARI(IL), SKEW(IL), PEAK(IL)
cc      COMMON  /COMXX/ F(2000)
cxx      DIMENSION  F(NPE-NPS+1,IL)
      INTEGER :: K, L, IL, NPS, NPE, MJ, MJ2
      REAL(8) :: X(NPE), A(K), E(MJ,IL), F(NPE-NPS+1,IL), EMEAN(IL),
     1           VARI(IL), SKEW(IL), PEAK(IL), COV(MJ2)
      REAL(8) :: SUM, COV1, SD
C
C
      ISTEP = 1
      ISW = IL
      LAGH = 100
      N = NPE - NPS - 1
      IF( LAGH .GE. N )     LAGH = N - 1
      LAG1 = LAGH + 1
      NMK = N - K
      IF( ISW .GT. 0 )     GO TO 20
C
      DO 10  I=NPS,NPE
cxx   10 E(I,1) = X(I)
      E(I,1) = X(I)
   10 CONTINUE
      IL = 1
      GO TO 36
C
C       ---  SEVERAL STEP AHEAD PREDICITON  ---
C
   20 CONTINUE
      CALL  PRDCT( X,A,K,L,IL,NPS,NPE,MJ,E )
C
C       ---  PREDICTION ERROR  ---
C
cxx      DO 30  II=NPS,NPE
      DO 31  II=NPS,NPE
         I = NPE-II+NPS
         DO 30  J=1,IL
         JJ = I-J+1
cxx   30 E(I,J) = X(I) - E(JJ,J)
      E(I,J) = X(I) - E(JJ,J)
   30 CONTINUE
   31 CONTINUE
      IF( IL .EQ. 1 )     GO TO 34
      DO 35  J=2,IL
      JJ = J-1
cxx      DO 35  I=1,JJ
      DO 33  I=1,JJ
      II = I+NPS-1
cxx   35 E(II,J) = 0.D0
      E(II,J) = 0.D0
   33 CONTINUE
   35 CONTINUE
   34 CONTINUE
cc      WRITE( 6,690 )
cc      IF( IPR .EQ. 0 )   GO TO 36
cc      WRITE( 6,670 )     (I,I=1,IL)
cc      DO 37  I=NPS,NPE
cc   37 WRITE( 6,640 )     I , (E(I,J),J=1,IL)
   36 CONTINUE
C
C       ---  MOMENT COMPUTATION  ---
C
      DO 50  KK=1,IL
C
      II = NPS+KK-1
      DO 40  I=II,NPE
      J = I - II + 1
cc   40 F(J) = E(I,KK)
cxx   40 F(J,KK) = E(I,KK)
      F(J,KK) = E(I,KK)
   40 CONTINUE
      NMK = NPE-NPS-(KK-2)
C
cc      CALL  MOMENT( F,NMK,EMEAN,VARI,SKEW,PEAK )
      CALL  MOMENT( F(1,KK),NMK,EMEAN(KK),VARI(KK),SKEW(KK),PEAK(KK) )
C
cc      IF( ISW .GT. 1 )     WRITE( 6,610 )   KK
cc      IF( ISW .EQ. 0 )     WRITE( 6,650 )
cc      WRITE( 6,620 )     EMEAN , VARI , SKEW , PEAK
C
C       ---  HISTOGRAM OF F(I)  ---
C
cc      SIG = DSQRT(VARI)
cc      CALL  GRAPH1( F,1,NMK,SIG )
   50 CONTINUE
C
C       ---  AUTOCORRELATION FUNCTION COMPUTATION  ---
C
      DO 100  KK=1,IL
C
      DO 70   II=1,LAG1
      JJ = NPS + KK - 1
      IE = NPE - II + 1
      SUM = 0.D0
      DO 60   I=JJ,IE
      J = I + II- 1
cxx   60 SUM = SUM + E(I,KK)*E(J,KK)
      SUM = SUM + E(I,KK)*E(J,KK)
   60 CONTINUE
cxx   70 COV(II) = SUM / (NPE-NPS-KK+2)
      COV(II) = SUM / (NPE-NPS-KK+2)
   70 CONTINUE
C
      COV1 = COV(1)
      DO 80   I=1,LAG1
cxx   80 COV(I) = COV(I) / COV1
      COV(I) = COV(I) / COV1
   80 CONTINUE
C
cc      IF( ISW .GT. 0 )   WRITE( 6,630 )     KK
cc      IF( ISW .EQ. 0 )   WRITE( 6,660 )
cc      WRITE( 6,680 )     (COV(I),I=1,50)
C
C       ---  AUTOCORRELATION FUNCTION DISPLAY  ---
C
      SD = NPE - NPS + 1
      SD = DSQRT( 1.D0/SD )
cc      CALL  GRAPH2( COV,LAG1,SD )
C
      IF( ISTEP .EQ. KK )     GO TO 110
C
  100 CONTINUE
C
C
  110 CONTINUE
cxx      RETURN
cxx  610 FORMAT( //1H ,I3,'-STEP AHEAD PREDICTION ERROR' )
cxx  620 FORMAT( 1H ,'MEAN       =',D15.8,/,' VARIANCE   =',D15.8,/,' SKEWN
cxx     1ESS   =',D15.8,/,' PEAKEDNESS =',D15.8 )
cxx  630 FORMAT( //1H ,'AUTOCORRELATION FUNCTION OF ',I3,'-STEP AHEAD PREDI
cxx     1CTION ERROR' )
cxx  640 FORMAT( 1H ,I5,5X,10D12.4 )
cxx  650 FORMAT( //,' ORIGINAL DATA' )
cxx  660 FORMAT( //1H ,'AUTOCORRELATION FUNCTION OF ORIGINAL DATA' )
cxx  670 FORMAT( 1H ,10X,10(4X,'J =',I2,3X) )
cxx  680 FORMAT( 1H ,10D13.5 )
cxx  690 FORMAT( ///1H ,45(1H-),2X,'<< J-STEP AHEAD PREDICTION ERROR >>',2X
cxx     1,45(1H-) )
      END
      SUBROUTINE  MOMENT( X,N,F1,F2,F3,F4 )
C
C          +--------------------+
C          ! MOMENT COMPUTATION !
C          +--------------------+
C
C     THIS SUBROUTINE COMPUTES MOMENTS.
C
C       INPUTS:
C          X:     ORIGINAL DATA VECTOR
C          N:     DATA LENGTH
C
C       OUTPUTS:
C          F1:    MEAN OF X
C          F2:    VARIANCE OF X
C          F3:    SKEWNESS OF X
C          F4:    PEAKEDNESS OF X
C
cxx      IMPLICIT  REAL * 8  ( F )
CC      DIMENSION  X(1)
cx	REAL*8  X(1)
cxx      REAL*8  X(N)
      INTEGER :: N
      REAL(8) :: X(N), F1, F2, F3, F4
      REAL(8) :: FN, FSUM, FF
C
      FN = N
      FSUM = 0.D0
      DO  10     I=1,N
cxx   10 FSUM = FSUM + X(I)
      FSUM = FSUM + X(I)
   10 CONTINUE
C
      F1 = FSUM / FN
C
      F2 = 0.D0
      F3 = 0.D0
      F4 = 0.D0
      DO  20     I=1,N
      FF = X(I) - F1
      F2 = F2 + FF*FF
      F3 = F3 + FF**3
cxx   20 F4 = F4 + FF**4
      F4 = F4 + FF**4
   20 CONTINUE
C
      F2 = F2 / FN
      F3 = F3 / (FN*F2*DSQRT(F2))
      F4 = F4 / (FN*F2*F2)
C
      RETURN
      END
      SUBROUTINE  PRDCT1( Z,A,M,L,IL,NPS,NPE,MJ,EZ )
C
C     THIS SUBROUTINE COMPUTES SEVARAL STEP AHEAD PREDICTION VALUE OF AN
C     AUTOREGRESSIVE MOVING AVERAGE MODEL.
C
C       INPUTS:
C          Z:      ORIGINAL DATA VECTOR
C          A:      AR-MA COEFFICIENTS
C          M:      AR-ORDER
C          L:      MA-ORDER
C          IL:     MAXIMUM SPAN OF LONG RANGE PREDICTION
C          NPS:    PREDICTION STARTING POSITION
C          NPE:    PREDICTION ENDING POSITION
C          MJ:     ABSOLUTE DIMENSION OF EZ
C
C       OUTPUT:
C          EZ:     PREDICTION VALUE MATRIX
C
cxx      IMPLICIT  REAL * 8  ( A-H , O-Z )
CC      REAL * 4  Z(1) , EZ(MJ,1)
cx      DIMENSION  Z(1) , EZ(MJ,1)
cx	DIMENSION  A(1)
cxx      DIMENSION  Z(NPE) , EZ(MJ,IL)
cxx      DIMENSION  A(M)
      INTEGER :: M, L, IL, NPS, NPE, MJ
      REAL(8) :: Z(NPE), A(M), EZ(MJ,IL)
      REAL(8) :: SUM
C
C
      DO  100     II=NPS,NPE
C
cxx      DO  90     KK=1,IL
      DO  91     KK=1,IL
      KKM1 = KK - 1
      SUM = 0.D0
      IF( KK .EQ. 1 )     GO TO 30
      DO  20     I=1,KKM1
      KI = KK - I
cxx   20 SUM = SUM + A(I)*EZ(II,KI)
      SUM = SUM + A(I)*EZ(II,KI)
   20 CONTINUE
   30 IF( KK .GT. M )     GO TO 50
      DO  40     I=KK,M
      I1 = II + KKM1 - I
cxx   40 SUM = SUM + A(I)*Z(I1)
      SUM = SUM + A(I)*Z(I1)
   40 CONTINUE
C
   50 IF( L .LE. 0 )     GO TO 90
      IF( KK .GT. L )     GO TO 90
      DO  60     I=KK,L
      I1 = M + I
      I2 = II + KKM1 - I
      IF( I2 .GE. II )     GO TO 60
      SUM = SUM + A(I1)*(Z(I2)-EZ(I2,1))
   60 CONTINUE
   90 EZ(II,KK) = SUM
   91 CONTINUE
  100 CONTINUE
C
      RETURN
      END
      SUBROUTINE  PRDCT2( Z,A,K,L,IL,NPS,NPE,MJ1,EZ )
C
C     THIS SUBROUTINE COMPUTES SEVERAL STEPS AHEAD PREDICTION VALUES OF
C     NON-LINEAR REGRESSION MODEL.
C
C       INPUTS:
C          Z:      ORIGINAL DATA VECTOR
C          A:      VECTOR OF AR-COEFFICIENTS
C          K:      ORDER OF THE AR MODEL
C          L:      THIS DUMMY VARIABLE IS NOT REFERENCED IN THIS SUBROUT
C          IL:     MAXIMUM SPAN OF LONG RANGE PREDICTION
C          NPS:    PREDICTION STARTING POSITION
C          NPE:    PREDICTION ENDING POSITION
C          MJ1:    ABSOLUTE DIMENSION OF EZ
C
C       OUTPUT:
C          EZ:     PREDICTION VALUE MATRIX
C
CC      IMPLICIT  REAL*8 ( A-D,O-Y )
cxx      IMPLICIT  REAL*8 ( A-H,O-Z )
cc      DIMENSION  Z(1) , A(1) , EZ(MJ1,1) , Y(20)
cx      DIMENSION  Z(1) , A(1) , EZ(MJ1,1) , Y(IL)
cxx      DIMENSION Z(NPE) , A(K) , EZ(MJ1,IL) , Y(IL)
cc      REAL * 4   CSTDMY, SD0DMY
      INTEGER :: K, IL, NPS, NPE, MJ1
      REAL(8) :: Z(NPE), A(K), EZ(MJ1,IL)
      REAL(8) :: Y(IL), CSTDMY, SD0DMY, SUM, XX, X
      COMMON     / BBB /  LAG1(50) , LAG2(50) , LAG3(50), CSTDMY, SD0DMY
CC      INTEGER  RETURN
C
cxx      DO 100  II=NPS,NPE
      DO 110  II=NPS,NPE
         DO 50  J1=1,IL
            JJ = J1-1
            SUM = 0.D0
            DO 40  J=1,K
               XX = 1.D0
               LAG = LAG1(J)
CC               ASSIGN 10 TO RETURN
CC               GO TO 200
CC   10          XX = XX*X
               X = 1.D0
               IF ( LAG .GT. 0 ) THEN
                  I = II+JJ-LAG
                  IF ( I .GE. II ) THEN
                     I = I - II + 1
                     X = Y(I)
                  ELSE
                     X = Z(I)
                  END IF
               END IF
               XX = XX*X
C
               LAG = LAG2(J)
CC               ASSIGN 20 TO RETURN
CC               GO TO 200
CC   20          XX = XX*X
               X = 1.D0
               IF ( LAG .GT. 0 ) THEN
                  I = II+JJ-LAG
                  IF ( I .GE. II ) THEN
                     I = I - II + 1
                     X = Y(I)
                  ELSE
                     X = Z(I)
                  END IF
               END IF
               XX = XX*X
C
               LAG = LAG3(J)
CC               ASSIGN 30 TO RETURN
CC               GO TO 200
CC   30          XX = XX*X
               X = 1.D0
               IF ( LAG .GT. 0 ) THEN
                  I = II+JJ-LAG
                  IF ( I .GE. II ) THEN
                     I = I - II + 1
                     X = Y(I)
                  ELSE
                     X = Z(I)
                  END IF
               END IF
               XX = XX*X
C
cxx   40       SUM = SUM + A(J)*XX
            SUM = SUM + A(J)*XX
   40       CONTINUE
cxx   50    Y(J1) = SUM
         Y(J1) = SUM
   50    CONTINUE
C
      DO  100  J=1,IL
      EZ(II,J) = Y(J)
  100 CONTINUE
  110 CONTINUE
CC      GO TO 300
C        -----  INTERNAL SUBROUTINE  -----
C
CC  200 X = 1.D0
CC      IF ( LAG .LE. 0 )     GO TO 220
CC      I = II+JJ-LAG
CC      IF ( I .GE. II )      GO TO 210
CC      X = Z(I)
CC      GO TO 220
CC  210 I = I - II + 1
CC      X = Y(I)
CC  220 GO TO RETURN, ( 10,20,30 )
C        ----------------------------------
C
CC  300 RETURN
C    L : DUMMY
      L = L
      RETURN
      END
cc      SUBROUTINE  SETLAG( K )
      SUBROUTINE  SETLAG( K,LAG1,LAG2,LAG3,LAG4,LAG5 )
C
C     THIS SUBROUTINE  PREPARES SPECIFICATION OF REGRESSORS (L1(I),L2(I)
C     (I=1,...,K) FOR THE FITTING OF (POLYNOMIAL TYPE) NON-LINEAR MODEL.
C     THE OUTPUTS ARE USED AS THE INPUTS TO SUBROUTINE SETX2.
C
C       INPUTS:
C          LAG1:    MAXIMUM TIME LAG OF LINEAR TERM
C          LAG2:    MAXIMUM TIME LAG OF SQUARED TERM
C          LAG3:    MAXIMUM TIME LAG OF QUADRATIC CROSS TERM
C          LAG4:    MAXIMUM TIME LAG OF CUBIC TERM
C          LAG5:    MAXIMUM TIME LAG OF CUBIC CROSS TERM
C
C       OUTPUTS:
C          K:       NUMBER OF REGRESSORS
C          (L1(I),L2(I),L3(I))  (I=1,K):     SPECIFICATION OF REGRESSORS
C
C              ......................................................
C              I-TH REGRESSOR IS DEFINED BY
C                   Z(N-L1(I)) * Z(N-L2(I)) * Z(N-L3(I))
C              WHERE  0-LAG TERM Z(N-0) IS REPLACED BY THE CONSTANT 1.
C              ......................................................
C
cc      REAL * 4   CSTDMY,SD0DMY
cxx      REAL * 8   CSTDMY,SD0DMY
      INTEGER :: K, LAG1, LAG2, LAG3, LAG4, LAG5
      REAL(8) :: CSTDMY, SD0DMY
      COMMON     / BBB /  L1(50),L2(50),L3(50),CSTDMY,SD0DMY
C
cc      READ( 5,1 )     LAG1,LAG2,LAG3,LAG4,LAG5
cc      WRITE( 6,2 )    LAG1,LAG2,LAG3,LAG4,LAG5
      IF(LAG1.LE.0)  GO TO 15
      DO 10  I=1,LAG1
      L1(I) = I
      L2(I) = 0
cxx   10 L3(I) = 0
      L3(I) = 0
   10 CONTINUE
   15 K = LAG1
C
      IF(LAG2.LE.0)   GO TO 30
      DO 20  I=1,LAG2
      K = K+1
      L1(K) = I
      L2(K) = I
cxx   20 L3(K) = 0
      L3(K) = 0
   20 CONTINUE
C
   30 IF(LAG3.LE.1)   GO TO 50
      LL = LAG3-1
cxx      DO 40  I=1,LL
      DO 41  I=1,LL
         I1 = I+1
         DO 40  J=I1,LAG3
         K = K+1
         L1(K) = I
         L2(K) = J
cxx   40    L3(K) = 0
         L3(K) = 0
   40    CONTINUE
   41 CONTINUE
   50 M  = K
C
      IF(LAG4.LE.0)   GO TO 65
      DO 60  I=1,LAG4
      K = K+1
      L1(K) = I
      L2(K) = I
cxx   60 L3(K) = I
      L3(K) = I
   60 CONTINUE
   65 CONTINUE
C
      IF(LAG5.LE.1)   GO TO 80
cxx      DO 70  I=1,LAG5
cxx         DO 70  J=I,LAG5
      DO 72  I=1,LAG5
         DO 71  J=I,LAG5
            DO 70  L=J,LAG5
            IF(I.EQ.J .AND. J.EQ.L)  GO TO 70
            K = K+1
            L1(K) = I
            L2(K) = J
            L3(K) = L
   70       CONTINUE
   71    CONTINUE
   72 CONTINUE
C
cc   80 WRITE( 6,3 )
   80 CONTINUE
cxx      IF( LAG1 .EQ. 0 )  GO TO 100
cxx      DO 90  I=1,LAG1
cc   90 WRITE( 6,4 )     I , L1(I)
cxx   90 CONTINUE
cxx  100 J = LAG1+1
cxx      IF( LAG2+LAG3 .EQ. 0 )   GO TO 120
cxx      DO 110  I=J,M
cc  110 WRITE( 6,5 )     I , L1(I) , L2(I)
cxx  110 CONTINUE
cxx  120 J = M+1
cxx      IF( LAG4+LAG5 .EQ. 0 )   GO TO 140
cxx      DO 130  I=J,K
cc  130 WRITE( 6,6 )     I ,L1(I) , L2(I) ,L3(I)
cxx  130 CONTINUE
cxx  140 CONTINUE
C
      RETURN
cxx    1 FORMAT( 16I5 )
cxx    2 FORMAT( /1H ,'LAG1 =',I3,5X,'LAG2 =',I3,5X,'LAG3 =',I3,5X,
cxx     *        'LAG4 =',I3,5X,'LAG5 =',I3 )
cxx    3 FORMAT( 1H ,4X,'M',5X,'REGRESSOR  Y(I,M)' )
cxx    4 FORMAT( 1H ,I5,5X,'Z(I-',I2,')' )
cxx    5 FORMAT( 1H ,I5,5X,'Z(I-',I2,') * Z(I-',I2,')' )
cxx    6 FORMAT( 1H ,I5,5X,'Z(I-',I2,') * Z(I-',I2,') * Z(I-',I2,')' )
C
      END
      SUBROUTINE  SETX2( Z,N0,L,K,MJ1,JSW,LAG,X )
C
C          +----------------------------------------+
C          ! MATRIX X SET UP (FOR NON-LINEAR MODEL) !
C          +----------------------------------------+
C
C     THIS SUBROUTINE PREPARES DATA MATRIX X FROM DATA VECTOR Z(I) (I=N0
C     N0+K+LAG) FOR THE FITTING OF NON-LINEAR AUTOREGRESSIVE MODEL.  X I
C     USED AS THE INPUT TO SUBROUTINE HUSHLD.
C
C       INPUTS:
C          Z:      ORIGINAL DATA VECTOR
C          N0:     INDEX OF THE END POINT OF DISCARDED FORMER OBSERVATIO
C                  (NEW OBSERVATION STARTS AT N0+LAG+1 AND ENDS AT N0+LA
C          L:      DIMENSION OF THE VECTOR OF NEW OBSERVATIONS
C          K:      NUMBER OF REGRESSORS
C          MJ1:    ABSOLUTE DIMENSION OF X
C          JSW:    =0   TO CONSTRUCT INITIAL L*(K+1) DATA MATRIX
C                  =1   TO AUGMENT ORIGINAL (K+1)*(K+1) MATRIX X BY AN
C                       L*(K+1) DATA MATRIX OF ADDITIONAL OBSERVATIONS
C          LAG:    MAXIMUM TIME LAG
C          KSW:    THIS DUMMY VARIABLE IS NOT REFERENCED IN THIS SUBROUT
C
C--  THE FOLLOWING VARIABLE SPECIFICATION IS GIVEN EITHER BY REDLAG OR S
C         (L1(I) , L2(I) , L3(I))  (I=1,K)
C
C               I-TH REGRESSOR IS DEFINED BY
C                    Z(N-L1(I)) * Z(N-L2(I)) * Z(N-L3(I))
C               WHERE 0-LAG TERM Z(N-0) IS AUTOMATICALLY REPLACED BY CON
C
C       OUTPUT:
C          X:      L*(K+1) MATRIX           IF  JSW = 0
C                  (K+1+L)*(K+1) MATRIX     IF  JSW = 1
C
CC      REAL * 8  X(MJ1,1)
CC      DIMENSION  Z(1)
cc      REAL * 4   CSTDMY, SD0DMY
cx      REAL * 8  X(MJ1,1) ,  Z(1) , ZTEM
cxx      REAL * 8  X(MJ1,K+1) ,  Z(N0+LAG+L) , ZTEM
cxx      REAL * 8  CSTDMY, SD0DMY
      INTEGER :: N0, L, K, MJ1, JSW, LAG
      REAL(8) :: Z(N0+LAG+L), X(MJ1,K+1)
      REAL(8) :: CSTDMY, SD0DMY, ZTEM
      COMMON     / BBB /  L1(50) , L2(50) , L3(50), CSTDMY, SD0DMY
C
      K1 = K + 1
      I0 = K1*JSW
      DO  10     I=1,L
      I1 = I + I0
      J1 = N0 + LAG + I
cxx   10 X(I1,K1) = Z(J1)
      X(I1,K1) = Z(J1)
   10 CONTINUE
C
      DO  70     II=1,K
      LL1 = L1(II)
      LL2 = L2(II)
      LL3 = L3(II)
cxx      DO  60     I=1,L
      DO  61     I=1,L
      ZTEM = 1.D0
      I1 = I + I0
      J1 = N0 + LAG + I
      IF( LL1 .EQ. 0 )     GO TO 40
      M1 = J1 - LL1
      ZTEM = ZTEM * Z(M1)
   40 IF( LL2 .EQ. 0 )     GO TO 50
      M2 = J1 - LL2
      ZTEM = ZTEM * Z(M2)
   50 IF( LL3 .EQ. 0 )     GO TO 60
      M3 = J1 - LL3
      ZTEM = ZTEM * Z(M3)
   60 X(I1,II) = ZTEM
   61 CONTINUE
   70 CONTINUE
C
      RETURN
C
      E N D
      SUBROUTINE  SETX4( Z,NO,L,K,MJ1,JSW,LAG,X )
C
C     THIS SUBROUTINE PREPARES DATA MATRIX X FROM DATA VECTOR Z(I) (I=NO
C     NO+K+L) FOR THE FITTING OF AUTOREGRESSIVE MODEL WITH POLYNOMIAL TY
C     VALUE FUNCTION.  X IS THEN USED AS THE INPUT TO SUBROUTINE HUSHLD.
C       INPUTS:
C          Z:      ORIGINAL DATA VECTOR
C          NO:     INDEX OF THE END POINT OF DISCARDED FORMER OBSERVATIO
C          L:      DIMENSION OF THE VECTOR OF NEW OBSERVATIONS
C          K:      NUMBER OF REGRESSORS
C          MJ1:    ABSOLUTE DIMENSION OF X
C          JSW:    =0   TO CONSTRUCT INITIAL L*(K+1) DATA MATRIX
C                  =1   TO AUGMENT ORIGINAL (K+1)*(K+1) MATRIX X BY AN
C                       L*(K+1) DATA MATRIX OF ADDITIONAL OBSERVATIONS
C          LAG:    MAXIMUM TIME LAG OF THE MODELS
C          N:      DATA LENGTH
C
C       OUTPUT:
C          X:      L*(K+1) MATRIX           IF   JSW = 0
C                  (K+1+L)*(K+1) MATRIX     IF   JSW = 1
C
C
cxx      IMPLICIT  REAL  * 8 (A-H,O-Z)
CC      REAL * 4 Z(1)
cx      DIMENSION  Z(1)
cx      DIMENSION  X(MJ1,1)
cxx      DIMENSION  X(MJ1,K+1) ,  Z(NO+LAG+L)
      INTEGER :: NO, L, K, MJ1, JSW, LAG
      REAL(8) :: Z(NO+LAG+L), X(MJ1,K+1)
      REAL(8) :: BN, Y, XX
cc      COMMON     / AAA /  N , M
      COMMON     / AAA /  N
C
C          M:      ORDER OF POLYNOMIAL OF MEAN VALUE FUNCTION
C
      M = K - LAG - 1
      K1 = K + 1
      I0 = JSW*K1
      M1 = M + 1
      LAG=K-M1
      BN = 2.D0/(N-LAG-1.D0)
cxx      DO 10  I=1,L
      DO 11  I=1,L
      Y= BN*(NO+I-1)-1.D0
      XX= 1.D0
      DO 10 J=1,M1
      II= I+I0
      X(II,J) = XX
cxx   10 XX = XX*Y
      XX = XX*Y
   10 CONTINUE
   11 CONTINUE
C
cxx      DO 20  I=1,L
      DO 21  I=1,L
      II = I+I0
      NN = NO+LAG+I
      X(II,K1) = Z(NN)
      DO 20   J=1,LAG
      NN = NN-1
      JJ = J+M1
cxx   20 X(II,JJ) = Z(NN)
      X(II,JJ) = Z(NN)
   20 CONTINUE
   21 CONTINUE
C
      RETURN
C
      END
      SUBROUTINE  SRTMIN( X,N,IX )
C
C       THIS SUBROUTINE ARRANGES X(I) (I=1,N) IN ORDER OF INCREASING
C       MAGNITUDE OF X(I)
C
C       INPUTS:
C          X:   VECTOR
C          N:   DIMENSION OF THE VECTOR
C       OUTPUTS:
C          X:   ARRANGED VECTOR
C          IND: INDEX OF ARRANGED VECTOR
C
cxx      IMPLICIT  REAL * 8  ( A-H , O-Z )
cx      DIMENSION  X(1) , IX(1)
cxx      DIMENSION  X(N) , IX(N)
      INTEGER :: N, IX(N)
      REAL(8) :: X(N), XMIN, XT
C
      NM1 = N - 1
      DO  30     I=1,N
cxx   30 IX(I) = I
      IX(I) = I
   30 CONTINUE
      DO  20     II=1,NM1
      XMIN = X(II)
      MIN = II
      DO  10     I=II,N
      IF( XMIN .LT. X(I) )     GO TO 10
      XMIN = X(I)
      MIN = I
   10 CONTINUE
      IF( XMIN .EQ. X(II) )     GO TO 20
      XT = X(II)
      X(II) = X(MIN)
      X(MIN) = XT
      IT = IX(II)
      IX(II) = IX(MIN)
      IX(MIN) = IT
   20 CONTINUE
C
      RETURN
      E N D
cc      SUBROUTINE  SUBSPC( B,K,N,IPR,EK )
cxx      SUBROUTINE  SUBSPC( B,K,N,IPR,EK,IND,C,C1,C2,OEIC,ESUM1,OMEAN,OM )
      SUBROUTINE  SUBSPC( B,K,N,EK,IND,C,C1,C2,OEIC,ESUM1,OMEAN,OM )
C
C       THIS SUBROUTINE PRODUCES BAYESIAN ESTIMATES OF PARTIAL CORRELATI
C       BY CHECKING ALL SUBSET REGRESSION MODELS.
C       ----------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS SUBROUTINE
C             BICOEF
C             SRTMIN
C       ----------------------------------------------------------------
C
C         INPUTS:
C           B:   LEAST SQUARES ESTIMATES OF PARTIAL CORRELATIONS
C           K:   DIMENSION OF VECTOR A
C           N:   NUMBER OF OBSERVATIONS USED FOR THE ESTIMATION OF A(I)
C           IPR: =0  TO SUPPRESS THE OUTPUTS
C                >0  TO PRINT OUT THE OUTPUTS
C
C         OUTPUTS:
C           B(I) (I=1,K):   BAYESIAN ESTIMATES OF PARTIAL CORRELATIONS
C           EK:   EQUIVALENT NUMBER OF FREE PARAMETERS IN THE BAYESIAN M
C
cxx      IMPLICIT  REAL * 8 ( A-H , O-Z )
cc      DIMENSION  B(1) , C(50) , D(50,50)
cc      DIMENSION  IND(50) , KND(50) , ESUM(50)
cc      DIMENSION  C1(50), C2(50), ESUM1(50)
cx      DIMENSION  B(1) , C(K) , D(K+1,K+1)
cxx      DIMENSION  B(K) , C(K) , D(K+1,K+1)
cxx      DIMENSION  IND(K) , KND(K+1) , ESUM(K+1)
cxx      DIMENSION  C1(K+1), C2(K), ESUM1(K+1)
      INTEGER :: K, N, IND(K)
      REAL(8) :: B(K), EK, C(K), C1(K+1), C2(K), OEIC, ESUM1(K+1),
     1           OMEAN, OM
      INTEGER :: KND(K+1)
      REAL(8) :: BICOEF, D(K+1,K+1), ESUM(K+1), CC, DN, SUM, EIC,
     1           SUMC, EXIC, OSUM
C
      CC = 1.D0 + DLOG(2.D0)
      K1 = K + 1
      DN = N
cxx      DO 10   I=1,K1
cxx      ESUM(I) = 0.D0
cxx      DO 10   J=1,K1
cxx   10 D(I,J) = 0.D0
      ESUM(1:K1) = 0.D0
      D(1:K1,1:K1) = 0.D0
C
C          SQUARE OF PARTIAL CORRELATIONS ( NORMALISED BY MULTIPLYING N
C
      DO 20   I=1,K
cxx   20 C(I) = B(I)*B(I)*DN
      C(I) = B(I)*B(I)*DN
   20 CONTINUE
C
C          ARRANGEMENT OF C(I) IN ORDER OF INCREASING MAGNITUDE
C
      CALL  SRTMIN( C,K,IND )
cc      IF( IPR .LE. 1 )     GO TO 60
cc      WRITE( 6,7 )
cc      DO  50     I=1,K
cc   50 WRITE( 6,6 )     I , IND(I) , C(I)
cc   60 CONTINUE
C
C          FIND THE MINIMUM OF EIC
C
      OEIC = CC*K
      SUM = 0.D0
      DO 30   I=1,K
      SUM = SUM + C(I)
      EIC = SUM + CC*(K-I)
cxx   30 IF( OEIC .GT. EIC )   OEIC = EIC
      IF( OEIC .GT. EIC )   OEIC = EIC
   30 CONTINUE
cc      WRITE( 6,604 )   OEIC
C
C--------  COMPUTATION OF EIC'S OF WHOLE SUBSET REGRESSION MODELS  -----
C
C          INITIAL SETTING
C
      DO 40   I=1,K
cxx   40 KND(I) = 0
      KND(I) = 0
   40 CONTINUE
      KND(K1) = 1
      SUM = 0.D0
      SUMC = 0.D0
      M = K
      IP = 0
      IQ = 0
C
  100 CONTINUE
C
C          -----  SPECIFICATION OF NEXT SUBSET  -----
cxx               DO 110   I=1,K
               DO 111   I=1,K
               IF( KND(I) .EQ. 0 )   GO TO 110
               KND(I) = 0
               GO TO 120
  110          KND(I) = 1
  111          CONTINUE
  120          CONTINUE
C          ------------------------------------------
C
  130 CONTINUE
      IF( IP .GT. K )   GO TO 200
      IF( KND(IP+1) .EQ. 0 )   GO TO 140
C
      IF( IQ .EQ. 0 )   GO TO 165
      IF( KND(IQ) .EQ. 1 )   GO TO 150
      IQ = IQ-1
      SUMC = SUMC + C(IQ+1)
C
      IF( SUMC + CC*(K-IP+IQ) .GT. OEIC + 40.D0 )   GO TO 180
      GO TO 150
C
  140 IP = IP+1
      IQ = IP-1
      SUMC = C(IP)
      IF( SUMC + CC .GT. OEIC + 40.D0 )   GO TO 200
C
  150 M = K-IP+IQ
      SUM = SUMC
      IF( IQ .EQ. 0 )   GO TO 165
      DO 160   I=1,IQ
      IF( KND(I) .EQ. 1 )   GO TO 160
      M = M-1
      SUM = SUM + C(I)
  160 CONTINUE
  165 CONTINUE
      EIC = SUM + CC*M - OEIC
      IF( EIC .GT. 40.D0 )   GO TO 100
      EXIC = DEXP( -0.5D0*EIC )
      ESUM(M+1) = ESUM(M+1) + EXIC
      DO 170   I=1,K
cxx  170 IF( KND(I) .EQ. 1 )   D(I,M+1) = D(I,M+1) + EXIC
      IF( KND(I) .EQ. 1 )   D(I,M+1) = D(I,M+1) + EXIC
  170 CONTINUE
      GO TO 100
C         --------------------------------------------
  180          DO 190   I=1,IP
cxx  190          KND(I) = 1
               KND(I) = 1
  190          CONTINUE
               KND(IP+1) = 0
               IP = IP+1
               IQ = IP-1
               GO TO 130
C         ---------------------------------------------
C
C--------------------------  WHOLE SUBSETS CHECKED  --------------------
C
  200 CONTINUE
cc      IF( IPR .GE. 2 )     WRITE( 6,8 )
cc      IF( IPR .GE. 2 )     WRITE( 6,607 )     (ESUM(I),I=1,K1)
      DO 201 I=1,K1
cxx  201 ESUM1(I) = ESUM(I)
      ESUM1(I) = ESUM(I)
  201 CONTINUE
C
C          MEAN OF NUMBER OF PARAMETERS
C
      OSUM = 0.D0
      SUM = ESUM(1)
      DO 210   I=1,K
      SUM = SUM + ESUM(I+1)
cxx  210 OSUM = OSUM + I*ESUM(I+1)
      OSUM = OSUM + I*ESUM(I+1)
  210 CONTINUE
      OMEAN = OSUM / SUM
      OM = OMEAN / K
cc      IF( IPR .GE. 2 )     WRITE( 6,608 )     OMEAN , OM
C
C       --  BINOMIAL TYPE DAMPER  --
C
      DO 220   I=1,K1
      J = I-1
      KMJ = K-J
cc  220 C(I) = BICOEF(K,J)*(OM**J)*((1.D0-OM)**KMJ)
cc      C(1) = 1.D0 / (1.D0 + K)
cxx  220 C1(I) = BICOEF(K,J)*(OM**J)*((1.D0-OM)**KMJ)
      C1(I) = BICOEF(K,J)*(OM**J)*((1.D0-OM)**KMJ)
  220 CONTINUE
      C1(1) = 1.D0 / (1.D0 + K)
      DO 221  I=1,K
cc  221 C(I+1) = C(I) * I / (1.D0 + K - I)
cxx  221 C1(I+1) = C1(I) * I / (1.D0 + K - I)
      C1(I+1) = C1(I) * I / (1.D0 + K - I)
  221 CONTINUE
cc      IF( IPR .GE. 2 )     WRITE( 6,609 )
cc      IF( IPR .GE. 2 )     WRITE( 6,607 )     (C(I),I=1,K1)
C
      DO 230   I=1,K1
cc  230 ESUM(I) = ESUM(I)*C(I)
cxx  230 ESUM(I) = ESUM(I)*C1(I)
      ESUM(I) = ESUM(I)*C1(I)
  230 CONTINUE
C
      SUM = 0.D0
      DO 240   I=1,K1
cxx  240 SUM = SUM + ESUM(I)
      SUM = SUM + ESUM(I)
  240 CONTINUE
C
cxx      DO 250   J=1,K1
      DO 251   J=1,K1
      DO 250   I=1,K
cc  250 D(I,J) = D(I,J)*C(J) / SUM
cxx  250 D(I,J) = D(I,J)*C1(J) / SUM
      D(I,J) = D(I,J)*C1(J) / SUM
  250 CONTINUE
  251 CONTINUE
C
C          WEIGHTS OF PARTIAL CORRELATIONS
C
cxx      DO 260   I=1,K
cc  260 C(I) = 0.D0
cxx  260 C2(I) = 0.D0
      C2(1:K) = 0.D0
cxx      DO 270   I=1,K
      DO 271   I=1,K
      DO 270   J=1,K1
cc  270 C(I) = C(I) + D(I,J)
cxx  270 C2(I) = C2(I) + D(I,J)
      C2(I) = C2(I) + D(I,J)
  270 CONTINUE
  271 CONTINUE
cc      IF( IPR .GE. 2 )     WRITE( 6,603 )
cc      IF( IPR .GE. 2 )     WRITE( 6,607 )     (C(I),I=1,K)
C
C          AVERAGING AND REARRANGEMENT OF PARTIAL CORRELATIONS
C
      EK = 1.D0
      DO 280   I=1,K
      J = IND(I)
cc      B(J) = B(J)*C(I)
cc  280 EK = EK + C(I)**2
      B(J) = B(J)*C2(I)
cxx  280 EK = EK + C2(I)**2
      EK = EK + C2(I)**2
  280 CONTINUE
cc      IF( IPR .LE. 1 )     RETURN
cc      WRITE( 6,602 )
cc      DO  290     I=1,K
cc  290 WRITE( 6,609 )     I , B(I)
C
      RETURN
C
cxx    3 FORMAT( 1H ,20I3 )
cxx    4 FORMAT( 1H ,3D20.10,2I5 )
cxx    6 FORMAT( 1H ,2I7,F13.5 )
cxx    7 FORMAT( 1H ,6X,'I IND(I)',4X,'N*B(I)**2' )
cxx    8 FORMAT( 1H ,'ESUM(I) (I=1,M+1)' )
cxx    9 FORMAT( 1H ,'***  BINOMIAL TYPE  ***' )
cxx  602 FORMAT( 1H ,'PARTIAL CORRELATIONS OF THE BAYESIAN MODEL',/,1H ,6X,
cxx     1'I',9X,'A(I)' )
cxx  603 FORMAT( 1H ,'FINAL BAYESIAN WEIGHTS OF PARTIAL CORRELATIONS' )
cxx  604 FORMAT( 1H ,'  OAIC =',D13.5 )
cxx  606 FORMAT( 1H ,'OSUM =',D20.10,5X,'SUM =',D20.10,5X,'COD =',D20.10 )
cxx  607 FORMAT( 1H ,10D13.5 )
cxx  608 FORMAT( 1H ,'OMEAN =',D15.8,5X,'OM =',D15.8 )
cxx  609 FORMAT( 1H ,I7,F13.5 )
      END
cc      SUBROUTINE  SBBAYS( X,D,K,N,IPR,MJ1,A,SD )
      SUBROUTINE  SBBAYS( X,K,N,IPR,MJ1,A,SD,EK,AIC,IND,C,C1,C2,B,
     *OEIC,ESUM,OMEAN,OM  )
C
C     THIS SUBROUTINE PRODUCES BAYESIAN MODEL BASED ON ALL SUBSET
C     REGRESSION MODELS USING THE OUTPUT OF SUBROUTINE REDUCT.
C       ----------------------------------------------------------------
C       THE FOLLOWING SUBROUTINES ARE DIRECTLY CALLED BY THIS SUBROUTINE
C             RECOEF
C             SDCOMP
C             SUBSPC
C       ----------------------------------------------------------------
C
C       INPUTS:
C          X:     N*(K+1) TRIANGULAR MATRIX,OUTPUT OF SUBROUTINE REDUCT
C          K:     NUMBER OF REGRESSORS OF THE BAYESIAN MODEL
C          N:     DATA LENGTH
C          IPR:   =0  TO SUPPRESS THE OUTPUTS
C                 =1  TO PRINT OUT FINAL RESULT
C                 =2  TO PRINT OUT INTERIM AND FINAL RESULTS
C          MJ1:   ABSOLUTE DIMENSION OF X
C
C       OUTPUTS:
C          A(I) (I=1,K):   REGRESSION COEFFICIENTS OF BAYESIAN MODEL
C          SD:    RESIDUAL VARIANCE
C
cxx      IMPLICIT  REAL * 8  (A-H , O-Z )
cc      DIMENSION  X(MJ1,1) , A(1) , D(1)
cc      DIMENSION  B(50) , G(50)
cc      DIMENSION  IND(50), C(50), C1(50), C2(50), ESUM(50)
cx      DIMENSION  X(MJ1,1) , A(1) , D(K)
cxx      DIMENSION  X(MJ1,K+1) , A(K) , D(K)
cxx      DIMENSION  B(K) , G(K)
cxx      DIMENSION  IND(K), C(K), C1(K+1), C2(K), ESUM(K+1)
      INTEGER :: K, N, IPR, MJ1, IND(K)
      REAL(8) :: X(MJ1,K+1), A(K), SD, EK, AIC, C(K), C1(K+1), C2(K),
     1           B(K), OEIC, ESUM(K+1), OMEAN, OM
      REAL(8) :: D(K), G(K), FN, SUM, BB
      K1 = K + 1
      FN = N
cc      IF( IPR .GE. 2 )     WRITE( 6,3 )
C
C          PARTIAL CORRELATIONS COMPUTATION
C
      SUM = X(K1,K1)**2
      DO 10   I=1,K
      J = K1-I
      SUM = SUM + X(J,K1)**2
      G(J) = DSQRT( SUM )
cxx   10 B(J) = X(J,K1)*X(J,J) / (G(J)*DABS(X(J,J)))
      B(J) = X(J,K1)*X(J,J) / (G(J)*DABS(X(J,J)))
   10 CONTINUE
C
C          PARTIAL CORRELATIONS OF BAYESIAN MODEL COMPUTATION
C
cc      CALL  SUBSPC( B,K,N,IPR,EK )
cxx      CALL  SUBSPC( B,K,N,IPR,EK,IND,C,C1,C2,OEIC,ESUM,OMEAN,OM )
      CALL  SUBSPC( B,K,N,EK,IND,C,C1,C2,OEIC,ESUM,OMEAN,OM )
C
C          MODIFICATION OF CROSS-PRODUCTS  X(I,K1) (I=1,K)
C
      DO 30   I=1,K
cc      B(I) = B(I)*X(I,I)*G(I) / DABS(X(I,I))
      BB = B(I)*X(I,I)*G(I) / DABS(X(I,I))
      D(I) = X(I,K1)
cc   30 X(I,K1) = B(I)
cxx   30 X(I,K1) = BB
      X(I,K1) = BB
   30 CONTINUE
C
C          REGRESSION COEFFICIENTS OF BAYSIAN MODEL
C
      CALL  RECOEF( X,K,K,MJ1,A )
C
      DO 40   I=1,K
cxx   40 X(I,K1) = D(I)
      X(I,K1) = D(I)
   40 CONTINUE
C
C          RESIDUAL VARIANCE AND AIC
C
cc      CALL  SDCOMP( X,A,D,N,K,MJ1,SD )
      CALL  SDCOMP( X,A,N,K,MJ1,SD )
C
      IF( IPR .EQ. 0 )     RETURN
      AIC = FN*DLOG(SD) + 2.D0*EK
cc      WRITE( 6,6 )     SD , EK , AIC
      RETURN
C
cxx    3 FORMAT( //1H ,18(1H-),/,' BAYESIAN PROCEDURE',/,1H ,18(1H-) )
cxx    5 FORMAT( 1H ,10D13.5 )
cxx    6 FORMAT( 1H ,'RESIDUAL VARIANCE',16X,'SD =',D19.8,/,1H ,'EQUIVALENT
cxx     1 NUMBER OF PARAMETERS  EK =',F10.3,/,1H ,32X,'AIC =',F15.3 )
      E N D

ext: f


text:
!! Copyright (C) 2008 Pavel Sakov
!!
!! This file is part of EnKF-Matlab. EnKF-Matlab is a free software. See
!! LICENSE for details.

! File:           parameters.f90
!
! Created:        31/08/2007
!
! Last modified:  08/02/2008
!
! Author:         Pavel Sakov
!                 CSIRO Marine and Atmospheric Research
!                 NERSC
!
! Purpose:        Fortran code for QG model. Reads and stores the model
!                 parameters.
!
! Description:
!
! Revisions:

module parameters_mod

  use utils_mod
  implicit none

  save

  public parameters_read

  integer, parameter, private :: FID_PRM = 111

  ! This is the place to change the resolution of the model, by
  ! changing the parameter MREFIN:
  ! To make the ocean a rectangular (at your own risk!) change nx1,
  ! ny1 to (2,4), and adjust domain size (Lx, Ly)...
  !
  integer, parameter, public :: MREFIN = 7
  integer, parameter, public :: NX1 = 2
  integer, parameter, public :: NY1 = 2
  integer, parameter, public :: M = NX1 * 2 ** (MREFIN - 1) + 1
  integer, parameter, public :: N = NY1 * 2 ** (MREFIN - 1) + 1

  real(8), parameter, public :: PI = 3.14159265358979323d0

  real(8), parameter, public :: rf_coeff = 0.1 ! Roberts filter coefficient
                                               ! for the leap-frog scheme

  real(8), public :: dt            ! time step
  real(8), public :: rkb           ! bottom friction
  real(8), public :: rkh           ! horizontal friction
  real(8), public :: rkh2          ! biharmonic friction
  real(8), public :: F             ! Froud number, (L / Rd)^2
  real(8), public :: r             ! factor in front of nonlinear
                                   ! advection term
  integer, public :: rstart = 0
  real(8), public :: tend
  real(8), public :: dtout
  integer, public :: verbose = VERBOSE_DEF
  character(STRLEN), public :: scheme
  character(STRLEN), public :: restartfname
  character(STRLEN), public :: outfname

  ! secondary parameters

  real(8), dimension(N, M), public :: CURLT

  real(8), public :: dx
  real(8), public :: dy

contains

  subroutine parameters_read(prmfname)
    character(STRLEN) :: prmfname
    integer :: j
    integer :: m1 = M - 1

    namelist /parameters/ &
         tend, &
         dtout, &
         dt, &
         rkb, &
         rkh, &
         rkh2, &
         F, &
         r, &
         verbose, &
         scheme, &
         rstart, &
         restartfname, &
         outfname

    call wopen(trim(prmfname), FID_PRM, 'old')
    read(FID_PRM, parameters)

    if (verbose > 0) then
       write(stdout, *) 'parameters:'
       write(stdout, *) '  tend =', tend
       write(stdout, *) '  dtout =', dtout
       write(stdout, *) '  dt =', dt
       write(stdout, *) '  RKB =', rkb
       write(stdout, *) '  RKH =', rkh
       write(stdout, *) '  RKH2 =', rkh2
       write(stdout, *) '  F =', F
       write(stdout, *) '  r =', r
       write(stdout, *) '  verbose =', verbose
       write(stdout, *) '  scheme = "', trim(scheme), '"'
       write(stdout, *) '  rstart =', rstart
       write(stdout, *) '  restart file = "', trim(restartfname), '"'
       write(stdout, *) '  output file = "', trim(outfname), '"'
    end if

    dx = 1.0d0 / dble(N - 1)
    dy = 1.0d0 / dble(M - 1)

    do j = 1, M
       CURLT(:, j) = - 2.0d0 * PI * sin(2.0d0 * PI * (j - 1) / m1)
    end do

    close(FID_PRM)
  end subroutine parameters_read

  subroutine parameters_get(N_out, M_out, tend_out, dtout_out, dt_out, &
       rkb_out, rkh_out, rkh2_out, F_out, r_out, verbose_out, &
       scheme_out, rstart_out, restartfname_out, outfname_out)

    integer, intent(out) :: N_out, M_out
    real(8), intent(out) :: dt_out
    real(8), intent(out) :: rkb_out
    real(8), intent(out) :: rkh_out
    real(8), intent(out) :: rkh2_out
    real(8), intent(out) :: F_out
    real(8), intent(out) :: r_out
    real(8), intent(out) :: tend_out
    real(8), intent(out) :: dtout_out
    integer, intent(out) :: verbose_out
    character(STRLEN), intent(out) :: scheme_out
    integer, intent(out) :: rstart_out
    character(STRLEN), intent(out) :: restartfname_out
    character(STRLEN), intent(out) :: outfname_out

    N_out = N
    M_out = M
    tend_out = tend
    dtout_out = dtout
    dt_out = dt
    rkb_out = rkb
    rkh_out = rkh
    rkh2_out = rkh2
    F_out = F
    r_out = r
    verbose_out = verbose
    scheme_out = scheme
    rstart_out = rstart
    restartfname_out = restartfname
    outfname_out = outfname
  end subroutine parameters_get

end module parameters_mod

max_stars_repo_name: brajard/DAPPER

max_issues_repo_name: brajard/DAPPER

max_forks_repo_name: brajard/DAPPER


text:
C    %W% %G%
c
c     This module consists of a set of subroutines for reading from
c     the direct-access history file

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine gethisw (fltvar)
      real fltvar

c     This subroutine reads the next word from the record buffer
c     without advancing the position pointer
c
      include 'tspinc/dahist.inc'
c
      fltvar = dabuff(posptr)
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine gethisc (charr, kount)

c     This subroutine reads in kount 4-byte words into the charr
c     array from the buffer for da file unit l8 and reads one or
c     more following records if the array requires data beyond the
c     current one.

c     By ANSI standard, a real param can't be passed up to a
c     character argument, so equivalencing to char*4 must be done
c     for each incoming real.
c
      character charr(*)*4, ch4*4
      integer kount
      real flt4
      equivalence (ch4,flt4)
c
      do la = 1,kount
        call gethisf (flt4, 1)
        charr(la) = ch4
      enddo
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine gethisi (intarr, kount)
      real intarr(*)
      integer kount

c     This subroutine reads in kount integers from the intarr array
c     from the buffer for da file unit l8 and reads one or more
c     following records if the array requires data beyond the current
c     one.

c     intarr is declared real so as to directly copy from the real
c     buffer.
c
      call gethisf (intarr, kount)
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine gethisf (fltarr, kount)
      real fltarr(*)
      integer kount

c     This subroutine reads in kount integers from the intarr array
c     from the buffer for da file unit l8 and reads one or more
c     following records if the array requires data beyond the current
c     one.

c     The convention is to have record & position pointers point to
c     file record & position in buffer next piece of data will
c     come from.
c     The minimun value for both pointers is 1.
c
      include 'tspinc/blkcom1.inc'
      include 'tspinc/dahist.inc'
c
      lb = posptr
      iend = lb + kount - 1
      if (iend .le. recsiz) then

c       All data we want resides in rest of buffer

        do la = 1,kount
          fltarr(la) = dabuff(lb)
          lb = lb + 1
        enddo
        posptr = lb
      else

c       We want more data than buffer holds.  Must read in one or
c       more records.
c       First get rest of what's in buffer now and advance recptr

        la = 1
        lb = posptr
        do while (lb .le. recsiz)
          fltarr(la) = dabuff(lb)
          la = la + 1
          lb = lb + 1
        enddo

c       If reading several records, can put incoming data directly
c       into caller's data array

        ka = la
        kb = ka + recsiz - 1

c       read full records only

        do while (kb .lt. kount)
          recptr = recptr + 1
          read (l8,rec=recptr) (fltarr(lc),lc = ka,kb)
          ka = ka + recsiz
          kb = kb + recsiz
        enddo

c       Only partial record needed, read into buffer first

        recptr = recptr + 1
        read (l8,rec=recptr) (dabuff(lc),lc = 1,recsiz)
        lb = 1
        do la = ka,kount
          fltarr(la) = dabuff(lb)
          lb = lb + 1
        enddo
      endif
      posptr = lb

c     if buffer "used up", read next one from file

      if (posptr .gt. recsiz) then
        recptr = recptr + 1
        posptr = posptr - recsiz
        read (l8,rec=recptr) (dabuff(lc),lc = 1,recsiz)
      endif
c
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine sethpos (recrd,postn)
      integer recrd, postn
c
c     This subroutine Moves the dir-acc file position to the recrd +
c     position pointer pair.
c
      include 'tspinc/blkcom1.inc'
      include 'tspinc/dahist.inc'
c
      if (recrd .ne. recptr) then
        recptr = recrd
        read (l8,rec=recptr) (dabuff(lb), lb = 1,recsiz)
      endif
      posptr = postn
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine shohpos (recrd,postn)
      integer recrd, postn
c
c     This subroutine returns the current dir-acc file recrd +
c     position pointer pair.
c
      include 'tspinc/dahist.inc'
c
      recrd = recptr
      postn = posptr
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine calchpos (oldrec,oldpos,newrec,newpos,addby)
      integer oldrec, oldpos, newrec, newpos, addby
c
c     This subroutine calculates a new file position from an older
c     position + an adding offset.
c
      include 'tspinc/dahist.inc'
c
      iaddrec = addby / recsiz
      iaddpos = addby - (iaddrec * recsiz)
      newrec = oldrec + iaddrec
      newpos = oldpos + iaddpos
      if (newpos .gt. recsiz) then
        newpos = newpos - recsiz
        newrec = newrec + 1
      endif
      return
      end

c   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-   -=-

      subroutine skiphis (numwrds)
      integer numwrds

c     This subroutine advances the current file pointer by nunmwrds.
c     It is similar to gethisi, except no data is read.
c
      include 'tspinc/dahist.inc'
c
      call calchpos (recptr,posptr,inewrec,inewpos,numwrds)
      call sethpos (inewrec,inewpos)
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
-I../../outdir/nv_small/vmod/nvdla/bdma
-I../../outdir/nv_small/vmod/nvdla/cacc
-I../../outdir/nv_small/vmod/nvdla/car
-I../../outdir/nv_small/vmod/nvdla/cbuf
-I../../outdir/nv_small/vmod/nvdla/cdma
-I../../outdir/nv_small/vmod/nvdla/cdp
-I../../outdir/nv_small/vmod/nvdla/cmac
-I../../outdir/nv_small/vmod/nvdla/csc
-I../../outdir/nv_small/vmod/nvdla/glb
-I../../outdir/nv_small/vmod/nvdla/nocif
-I../../outdir/nv_small/vmod/nvdla/pdp
-I../../outdir/nv_small/vmod/nvdla/retiming
-I../../outdir/nv_small/vmod/nvdla/rubik
-I../../outdir/nv_small/vmod/nvdla/sdp
-I../../outdir/nv_small/vmod/nvdla/top
-I../../outdir/nv_small/vmod/nvdla/csb_master
-I../../outdir/nv_small/vmod/nvdla/cfgrom
-I../../outdir/nv_small/vmod/rams/synth
-I../../outdir/nv_small/vmod/vlibs
-I../../outdir/nv_small/vmod/include
-v ../../outdir/nv_small/vmod/vlibs/RANDFUNC.vlib
-v ../../outdir/nv_small/vmod/vlibs/nv_assert_no_x.vlib
-v ../../outdir/nv_small/vmod/nvdla/nocif/NV_NVDLA_XXIF_libs.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_248X82_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_32X192_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_32X224_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_256X80_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X256_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_160X82_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X226_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_256X144_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_32X288_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_32X256_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_20X288_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_248X144_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_80X15_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_60X168_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_20X80_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_64X226_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_64X288_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_32X32_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_128X6_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_128X11_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_64X10_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_16X256_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X72_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_256X4_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_256X7_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_256X8_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_64X116_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_32X16_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_256X11_GL_M4_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_160X16_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_80X14_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X16_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_160X144_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X288_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_80X66_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_160X65_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_60X22_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_80X9_GL_M2_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_64X66_GL_M1_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_128X18_GL_M2_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_16X272_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMDP_16X64_GL_M1_E2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_512X64_GL_M4_D2.v
-v ../../outdir/nv_small/vmod/rams/model/RAMPDP_256X64_GL_M2_D2.v

-DNO_PLI_OR_EMU
-DNO_PLI
-DDESIGNWARE_NOEXIST
-DSYNTHESIS
-Wno-moddup
-Wno-fatal
--top-module NV_nvdla
../../outdir/nv_small/vmod/nvdla/top/NV_nvdla.v

max_stars_repo_name: minisparrow/freedom

max_issues_repo_name: minisparrow/freedom

max_forks_repo_name: minisparrow/freedom


text:
*DECK ALI
      FUNCTION ALI (X)
C***BEGIN PROLOGUE  ALI
C***PURPOSE  Compute the logarithmic integral.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C5
C***TYPE      SINGLE PRECISION (ALI-S, DLI-D)
C***KEYWORDS  FNLIB, LOGARITHMIC INTEGRAL, SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C ALI(X) computes the logarithmic integral; i.e., the
C integral from 0.0 to X of (1.0/ln(t))dt.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  EI, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900326  Removed duplicate information from DESCRIPTION section.
C           (WRB)
C***END PROLOGUE  ALI
C***FIRST EXECUTABLE STATEMENT  ALI
      IF (X .LE. 0.0) CALL XERMSG ('SLATEC', 'ALI',
     +   'LOG INTEGRAL UNDEFINED FOR X LE 0', 1, 2)
      IF (X .EQ. 1.0) CALL XERMSG ('SLATEC', 'ALI',
     +   'LOG INTEGRAL UNDEFINED FOR X = 1', 2, 2)
C
      ALI = EI (LOG(X) )
C
      RETURN
      END


text:
Volleyball is a sport that requires a handful of people and, a ball and a net.  Get the ball from a sports store, grab some friends and hit one of our local places to play!

 Places to Play

  Walnut Park has a sand lot
  on the south east side of Bainer Hall at lunch
  Community Park north of the solar panels.
  Open gym volleyball is played in Davis Senior High, two or three times a week, depending on the season. The floor fee is $5. To get more details, check the Event Calendar at http://www.sacvolleyball.com.  Update April 2008: The hours have changed.   On Fridays and Tuesdays, they play 9pm to 12am at the Davis Senior High Gym.  On Sundays, they play 5pm to 7pm also at he Davis Senior High Gym.  Dean (dinokitadani@yahoo.com) runs the open gym.
  The ARC and the Pavilion also provides several places to play. If you ask one of the member service assistants, the people in the blue uniforms, to setup a volleyball court in the main gym. If you feel that it is too busy in the main gym then you can ask them to setup a court in the Pavilion as well.
  Davis Athletic Club Up until this summer there was an open gym with Volleyball games Tuesday and Sunday from 6 to 9 PM, and Thursday (still ongoing) from about 6 (6:30) to 8 PM. Due to some repairs needed in the gym (has to do with protecting the lights on the ceiling), the management suspended the volleyball sessions on Tuesday and Sunday (summer 2010). Its supposed to resume in September. To play you either need a membership or you can buy pinch cards.

 Organizations
  Davis Volleyball Club

 Teams
  UC Davis Mens Volleyball

20080425 23:48:47 nbsp Theres about 24 players on friday nights which allows 4 6man teams.   There are people that are starting out to people who are experienced.  All are welcome and more nice that other places where theyre more clickish. Users/atwong

20080425 23:50:07 nbsp Oh yeah...  the median age is around 30 to 50. Users/atwong

20100818 03:13:46 nbsp Davis Athletic Club is the only place in Davis where, to my knowledge, you can drop in for Volleyball three times a week: Tuesday, Thursday and Sunday. To be honest, its the only reason I bought a membership, so I hope they keep it that way. Neat atmosphere, great games.  Users/PassingThrough

20100819 13:41:46 nbsp Hi! I am new to Davis and am looking for a place/time to play some decent/semicompetitive pickup sand volleyball.  Any advice??

Thanks :)
Gail Users/gail.dasilva

20100819 13:43:35 nbsp PassingThrough, I missed your comment before!  I also play indoor... what time do people play??  Thanks! Users/gail.dasilva

20100819 17:41:16 nbsp Hi Gail, right now at DAC, they play on Thursday evening, in the soccer room. I cant remember the time but it starts around 6 PM if I remember well. I would call the club before showing up. :)

As for the Tuesday and Sunday sessions, the games start at 6 PM (include the time needed to set up the net and to warm up). At the moment theyre repairing the lights (what protects them) in the gym, so volleyball on these days will resume in September. Hope this helps! Users/PassingThrough

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search