c
      subroutine out_pseudoPilot(ilevUD,rAz,rp,rHLOS,rNodata,rStdErr
     +               , kyear,kmonth,kday,khour,kminute
     +			, rLat,rLon,l_open_unit,kunit
     +               , id_dwlchannel, channel_mask,nprof)

      integer ilevUD
      real rAz
      real rp(iLevUD)
      real rHLOS(iLevUD)
      real rStdErr(iLevUD)
      integer kyear
      integer kmonth
      integer khour
      integer kminute
      real rLat
      real rLon
c      integer nprof
      logical l_open_unit
      integer kunit
      integer id_dwlchannel
      integer channel_mask(ilevUD)
C
C       written 29-4-1999 by Siebren de Haan
C       for       
C       input file : LIDAR-bufr
C       output     : pseudo_pilot-bufr
C 
C       modified 11-09-2002 by David Tan
C       for
C       input files: LIPAS ascii output
C
C       modified 01-11-2003 by David Tan
C       for
C       dwl channel ident stored in radiosonde type
C
	parameter	(jwork=80000, jbufl=20000, kelem = 20000, jbyte=80000)
	parameter	(jsup=9, jsec0=3, jsec1=40, jsec2=64, jsec3=4)
	parameter	(jsec4=2, jkey=46, jelem=20000, kvals=80000)

	integer		kbuff(jbufl),kbufr(jbufl),ksup(jsup),koutdata(jsup)
	integer		ksec0(jsec0),ksec1(jsec1),ksec2(jsec2),ksec3(jsec3)
	integer		ksec4(jsec4)
c	integer		isec1(jsec1),isec2(jsec2),isec3(jsec3)
c	integer		isec4(jsec4)
	integer		ktdlst(jelem)

	real		values(kvals),outvals(kvals)

	character*64	cnames(kelem)
	character*24	cunits(kelem)
	character*80	cvals(kvals)
	character*80 inname,outname

	real MAXSIGMA
	integer iGROSS,iSIGMA

        real rAz2
        real rLon2
c
c        write(6,*)'Starting bufr output...'
c        write(6,*)'id_dwlchannel=',id_dwlchannel
c

c
c	iNumObs=0
c

	iret=0
	nargs = iargc()

C       MAXSIGMA
	MAXSIGMA=50.0

C       missing data indicator
	rmdi=1.7E+38
	itot=0
	ngross=0
	nMAXSIGMA=0
c
c	nmess=0
	nmess=1

C     If first profile then 
C     open output unit
      if ( l_open_unit ) then
	l_open_unit = .false.
	if (nargs .lt. 1) then
	   write(*,*)'addnoise(out_pseudoPILOT) <BUFR output-file> ...'
	   stop
	endif
	call getarg(1,outname)
	call PBOPEN(kunit,outname,'w',iret)
c	print*,' kunit = ',kunit
	if (iret .eq. -1) STOP 'open failed on bufr.out'
	if (iret .eq. -2) STOP 'invalid file name'
	if (iret .eq. -3) STOP 'invalide open mode specified'
c        kunit = 10
c	print*,' kunit after pbopen = ',kunit
      endif


c	kval=ksup(5)  ! not used in what follows
c	nmess=ksup(6) ! set to 1 above
c
c	nlev=20			! number of levels
	nlev=ilevUD		! number of levels
c
c	iel=12+1+5*nlev		! khead RF DATA
c	1    +1+1+5*nlev+12	! QC RF QC-DATA
c	2    +1+1+1+5*nlev+12	! QC GC GA QC-DATA
	iel=12+1+5*nlev		! khead RF DATA
     +       +1+1+5*nlev+12    ! QC RF QC-DATA
     +       +1+1+1+5*nlev+12	! QC GC GA QC-DATA

c	do n=1,ksup(6)
	do n=1,nmess
	   itot=(n-1)*iel	   
c	   ival0=(n-1)*kel
C	set the header
	   outvals(itot+1)=99	! WMO block number
	   outvals(itot+2)=187	! WMO station number -> satellite
C	   outvals(itot+3)=values(ival0+2)
C	   outvals(itot+4)=values(ival0+3)
	   outvals(itot+3)=id_dwlchannel ! Sonde type - proxy for dwl_channel
	   outvals(itot+4)=rmdi ! Radiosonde computational method
	   
c	   outvals(itot+5)=values(ival0+4)
c	   outvals(itot+6)=values(ival0+5)
c	   outvals(itot+7)=values(ival0+6)
c	   outvals(itot+8)=values(ival0+7)
c	   outvals(itot+9)=values(ival0+8)
c	   outvals(itot+10)=values(ival0+10)
c	   outvals(itot+11)=values(ival0+11)
	   outvals(itot+5)=kyear	! Year
	   outvals(itot+6)=kmonth	! month
	   outvals(itot+7)=kday	        ! day
	   outvals(itot+8)=khour	! hour
	   outvals(itot+9)=kminute	! minute
	   outvals(itot+10)=rLat	! Latitude

c           if( rLon .gt. 540.).or.( rLon.le.-180.)then
c              write(*,*)'Bad longitude:',rLon
c              stop
c           else
c              rLon2 = rLon
c              if( rLon2 .gt. 180. ) rLon2 = rLon2 - 360.
c           endif
           rLon2 = rLon
           do while ( rLon2 .gt. 180. )
              rLon2 = rLon2 - 360.
           enddo
           do while ( rLon2 .le. -180. )
              rLon2 = rLon2 + 360.
           enddo

	   outvals(itot+11)=rLon2	! longitude
c
	   outvals(itot+12)=rmdi ! elevation is missing
c
c	   outvals(itot+13)=20
	   outvals(itot+13)=nlev

	
C       copy data
c	   ival=12
	   iout=13
	   ilev=1
	   do while (ilev.le.nlev)
	      outvals(itot+iout+1)=rp(ilev) ! pressure
	      
c	      iSIGMA=ival0+215+2*(ilev-1)+1 ! sigma position
c	      iGROSS=ival0+215+2*(ilev-1)+2 ! gross error position
c	      iAZIM=ival0+ival+2 ! azim position
c	      iHLOS=ival0+ival+3 ! HLOS position

C	      write(*,*) iGROSS,iSIGMA,iAZIM,iHLOS
C	      write(*,*) values(iGROSS),values(iSIGMA),values(iAZIM),values(iHLOS)
c
c	      if ( (values(iGROSS).eq.0.0).and.
C GJM
	      if ( rHLOS(ilev) .EQ. rNodata) rHLOS(ilev) = rmdi
	      if (
     +             (rStdErr(ilev).lt.MAXSIGMA) ) then	     
				! no gross error and sigma < MAXSIGMA
		 if (abs(rHLOS(ilev)).lt.1.E30) then
c		    if (rHLOS(ilev).gt.0.0) then
c		       outvals(itot+iout+2)=1.0 ! sign of HLOS
c		    else
c		       outvals(itot+iout+2)=0.0
c		    endif
c		    outvals(itot+iout+2)=0.0
		    outvals(itot+iout+2)=rmdi ! Vert sounding significance RMDI
c		    if (rHLOS(ilev).gt.0.0) then
c		       outvals(itot+iout+4)=rAz ! azim of HLOS
c		    else
c		       outvals(itot+iout+4)=rAz + 180.
c		    endif
c
                    rAz2 = rAz
c		    if (rHLOS(ilev).lt.0.0) rAz2=rAz2+180.
c		    rAz2 = AMOD( rAz2, 360. )
		    outvals(itot+iout+4)=rAz2
c
c		    outvals(itot+iout+5)=abs(rHLOS(ilev)) ! abs(los)
                    if(channel_mask(ilev).eq.id_dwlchannel)then
		       outvals(itot+iout+5)=    rHLOS(ilev)  !     los
                    else
		       outvals(itot+iout+5)=    rmdi  !     los
                    end if
		 else
		    outvals(itot+iout+2)=rmdi
		    outvals(itot+iout+5)=rmdi
		    outvals(itot+iout+4)=rmdi
		 endif
		 outvals(itot+iout+3)=rStdErr(ilev)*1000 ! std dev
c		 outvals(itot+iout+4)=rAz ! azim
	      else
		 goto 200
	      endif

	      goto 201
 200	      ngross=ngross+1	! gross error
	      if (rStdErr(ilev).lt.MAXSIGMA) nMAXSIGMA=nMAXSIGMA+1	     

	      outvals(itot+iout+2)=rmdi
	      outvals(itot+iout+3)=rmdi
	      outvals(itot+iout+4)=rmdi
	      outvals(itot+iout+5)=rmdi
	      
 201	      continue

	      iout=iout+5
c	      ival=ival+3
	      ilev=ilev+1
	   end do
C       data-presence indicator
	   iout=iout+1
	   outvals(itot+iout)=0.0 ! QUALITY INFORMATION FOLLOW
	   iout=iout+1
	   outvals(itot+iout)=12+5*nlev ! EXTENDED DELAYED .... 
	   do itemp=1,12+5*nlev
	      iout=iout+1
	      outvals(itot+iout)=0.0
	   end do
C       data-confidence indicator
	   iout=iout+1
	   outvals(itot+iout)=98.0 ! GENERATING CENTRE 
	   iout=iout+1
	   outvals(itot+iout)=10.0 ! GENERATING APPLICATION
	   iout=iout+1
	   outvals(itot+iout)=12+5*nlev ! EXTENDED DELAYED .... 
	   do itemp=1,12+5*nlev
	      iout=iout+1
	      outvals(itot+iout)=100.0
	   end do
	end do
	itot=itot+iout

	
	do itemp=1,itot
	   if (outvals(itemp).gt.1.E30) then
	      outvals(itemp)=rmdi
	   endif
	end do
C	write(*,*) itot,iout,ngross,nlev,nmess,ksup(6)
c	write(*,*) itot,iout,ngross,nlev,nmess

C	1 Initialization
C       1.1 integers
	khead=12
	ktdlen=22
	nset=3

C       1.2 Unexpanded data descriptors
	ktdlst( 1)=301001
C	ktdlst( 2)=002201	! Simulated satellite instrument
C	ktdlst( 3)=002022	! Satellite data processing technique used
	ktdlst( 2)=002011	! Sonde type - proxy for dwl channel
	ktdlst( 3)=002012	! Radiosonde computational method
	ktdlst( 4)=301011	! Year, month, day

	ktdlst( 5)=301012
	ktdlst( 6)=301022
	ktdlst( 7)=105000	! EXTENDED DELAYED DESCRIPTOR REPLICATION FACTOR
	ktdlst( 8)=031001
	ktdlst( 9)=007004	! Pressure
	ktdlst(10)=008001
	ktdlst(11)=010003
	ktdlst(12)=011001
c	ktdlst(13)=011002
	ktdlst(13)=011003       ! U-comp for signed LOS
	ktdlst(14)=222000
	ktdlst(15)=101000
	ktdlst(16)=031002	! EXTENDED DELAYED DESCRIPTOR REPLICATION FACTOR
	ktdlst(17)=031031
	ktdlst(18)=001031	! Generating centre
	ktdlst(19)=001201	! Generating application
	ktdlst(20)=101000
	ktdlst(21)=031002	! EXTENDED DELAYED DESCRIPTOR REPLICATION FACTOR
	ktdlst(22)=033007	! Confidence
        
C       1.3 Section 0 contents
	ksec0(1)=0		! TOTAL LENGTH OF SECTION 0
	ksec0(2)=0		! TOTAL LENGTH OF BUFR MESSAGE
	ksec0(3)=3		! BUFR EDITION NUMBER

C       1.4 Section 1 contents
	ksec1( 1)=18		! TOTAL LENGTH OF SECTION 1
	ksec1( 2)=3		! BUFR EDITION NUMBER
	ksec1( 3)=98		! ORIGINATING CENTRE
	ksec1( 4)=1		! UPDATE SEQUENCE NUMBER
	ksec1( 5)=0		! FLAG other data & uncompressed observed -> 64 compress
c	ksec1( 6)= 12		!  ers-1 and the like            
	ksec1( 6)=  2		!  Upper Air Soundings (PILOT etc)
	ksec1( 7)=187           !  91 if Land Pilot used for DWL pseudoPilot
	ksec1( 7)= 91           !  91 if Land Pilot used for DWL pseudoPilot
c
c	ksec1( 8)=isec1(8)	! VERSION NUMBER OF LOCAL TABLE USED
c	ksec1( 9)=isec1(9)
c	ksec1(10)=isec1(10)
c	ksec1(11)=isec1(11)
c	ksec1(12)=isec1(12)
c	ksec1(13)=isec1(13)
c	ksec1(14)=isec1(14)	! BUFR MASTER TABLE( ZERO) FOR METEOROLOGICAL DATA)
c
        ksec1( 8)=1       ! VERSION NUMBER OF LOCAL TABLE USED
        if ( kyear .le. 2000 ) then
           ksec1( 9)=kyear-1900 ! YEAR
        else
           ksec1( 9)=kyear-2000 ! YEAR
        endif
        ksec1(10)=kmonth  ! MONTH
        ksec1(11)=kday    ! DAY
        ksec1(12)=khour   ! HOUR
        ksec1(13)=kminute ! MINUTE
        ksec1(14)=0    ! BUFR MASTER TABLE( ZERO) FOR METEOROLOGICAL DATA)
c
	ksec1(15)=2		! VERSION NUMBER OF MASTER TABLE USED
	ksec1(16)=0		! LOCAL ADP CENTRE INFO.

C       1.5 Section 2 contents  (empty)
	ksec2(1)=0

C       1.6 Section 3 content
	ksec3(1)=0		! TOTAL LENGTH OF SECTION 3
	ksec3(2)=0		! RESERVED
	ksec3(3)=nmess       	! NUMBER OF SUBSETS  (Reports to be bufrized)
	ksec3(4)=64		!  64 FOR COMPRESSION/   0 MANY SUBSETS when other data

C       1.7 Section 4
	ksec4(1)=0
	ksec4(2)=0
	

C       2.2 set kdata
	koutdata(1)=nlev		! number of present levels
	koutdata(2)=nlev*5+khead ! number of non missing values in 1.bitmap
	koutdata(3)=nlev*5+khead ! number of non missing values in 1.bitmap

	kel=iel
	nbuot=nmess*kel

c
c        write(6,*)'Encoding bufrmessage...'
c


c        call bupkey( key,ksec1,ksec2,kerr )
c	if (kerr .ne. 0) then
c	   write(*,*)' key error - kerr = ',kerr
c	endif

c	if ((nprof .eq. 81) .or. (nprof.eq.82))then
c	   write(*,*) 'outvals ',(outvals(kk),kk=1,150)
c	   stop
c	endif
c	write(*,*) outvals(133)
c	write(*,*) 'itot ',itot
c	do ppp=1,kvals
c	  if (outvals(ppp) .lt. -95000) STOP
c	enddo

	call bufren( ksec0,ksec1,ksec2,ksec3,ksec4,
     +   ktdlen,ktdlst,nset,koutdata,kel,
     +   nbuot,outvals,cvals,kbufl,kbufr,kerr)

	if (kerr .ne. 0) then
	   write(*,*)' packing error - kerr = ',kerr
	endif
	
C	write(*,*)'************************************'
C	write(*,*)'******  END OF BUFR ENCODING  ******'
C	write(*,*)'************************************'
	
c	write(*,*) 'kbufl=',kbufl
	
	kbufl = kbufl * 4
c       
c       write out bufrmessage to file
c       

c
c        write(6,*)'Writing bufrmessage to',kunit
c

	call pbwrite(kunit,kbufr,kbufl,kret)
C
c	go to 300	! read next subset done elsewhere

C       exit OK
 1000	continue

	if(kret.eq.-1) then
	   print*,' Error in writing bufred data  ',kret
	else
c	   print*,kbufl,' Bytes written to kunit ',kunit,kret
	endif
	
	if( ngross .gt. 0 ) write(*,*) 'ngross   =',ngross
	if( nMAXSIGMA .gt. 0 ) write(*,*) 'nMAXSIGMA=',nMAXSIGMA

c        write(6,*)'Closing from out_pseudoPilot...',kunit
c	call PBCLOSE(kunit,iret)
	

	end