      program calmedstd
C.........................................................Aug. 2021
C Calculate median, std of WU, WL, WE for SSN=0, kp>=3
      DIMENSION CC(2000),WULE(4,2000),iut(2000),wulems(15,3,3) !WU(med-std,med+std),Wl..WE.. 
	+,xkp(19)
	CHARACTER*64 infile,outfile,outmedstd(15,3)
	CHARACTER*6 YRMNDY(2000)
C     +,YYMMDD
	DATA XKP/3.0,3.3,3.7,4.0,4.3,4.7,5.0,5.3,5.7,6.0,6.3,6.7
	+,7.0,7.3,7.7,8.0,8.3,8.7,9.0/
C
      infile='wulessn0kp3n.txt' !1858 events WU, WL, WE
	outfile='wulessn0kp3r.txt'  ! output events <med+-std> except outliers
	outmedstd='wulemedstd.txt'	!
	OPEN(11,file=outmedstd)
      OPEN(12,file=outfile)
C
      OPEN(10,file=infile,action='READ')
	DO nn=1,2000
	read(10,*,end=49) YRMNDY(nn),iut(nn),(WULE(k,nn),k=1,4)
	nend=nn
	ENDDO
   49	close(unit=10)
C
      DO 52 nkp=1,10
	zkp=xkp(nkp)
C
	DO 53 nw=1,3	 ! WU, WL, WE
C
	icnt=0
      DO 50 nn=1,nend
	wun=WULE(nw,nn)
	wkp=WULE(4,nn)
C
	          IF (nkp.lt.10) THEN
	            if(abs(wkp-zkp).lt.0.1) then
	icnt=icnt+1
	if (nw.ne.2) then
	CC(icnt)=wun
	               else
	CC(icnt)=-wun
	endif
	                endif
	                                       ELSE
	if (wkp.ge.zkp) then
	icnt=icnt+1
	if (nw.ne.2) then
	CC(icnt)=wun
	               else
	CC(icnt)=-wun
	endif
	endif
                               	ENDIF
   50 CONTINUE
C
        if (icnt.gt.1) then
	call smed1std3(icnt,CC,SMED,STD,STD1,STD2)
	  else
	if (icnt.eq.0) wun=0.
	smed=wun
	std=0.
	  endif
C
	bot=smed-STD
	top=smed+std
      if (nw.eq.2) then
	bot1=bot
	bot=-top
	smed=-smed
	top=-bot1
	endif
	wulems(nkp,nw,1)=top
	wulems(nkp,nw,2)=smed
	wulems(nkp,nw,3)=bot
  53  continue
C
C output top, med,bot
      write(*,30) zkp,icnt,(wulems(nkp,1,k),k=1,3)
	+,(wulems(nkp,2,k),k=1,3),(wulems(nkp,3,k),k=1,3)	
	write(11,30) zkp,icnt,(wulems(nkp,1,k),k=1,3)
	+,(wulems(nkp,2,k),k=1,3),(wulems(nkp,3,k),k=1,3)	
  30  format(1X,F4.1,1X,I4,9(1X,F8.2))
  52	CONTINUE
      close(unit=11)
C
C Output events between top and bot
       kout=0
            DO 60 nn=1,nend
		wun=WULE(1,nn)
	wln=WULE(2,nn)
	wen=WULE(3,nn)
 	wkp=WULE(4,nn)
C
      DO 62 nkp=1,10
	zkp=xkp(nkp)
	        IF (nkp.lt.10) THEN ! 1>>>>>>>>>>>>>
	IF (wkp.eq.zkp) THEN		  ! 2
	wutop=wulems(nkp,1,1)
	wubot=wulems(nkp,1,3)
	wltop=wulems(nkp,2,1)
	wlbot=wulems(nkp,2,3)
      IF ((wun.le.wutop).and.(wun.ge.wubot).and.
	+	(wln.le.wltop).and.(wln.ge.wlbot)) THEN  ! 3
	 kout= kout+1
	write(*,31) kout,YRMNDY(nn),iut(nn),(WULE(k,nn),k=1,4)
	write(12,31) kout,YRMNDY(nn),iut(nn),(WULE(k,nn),k=1,4)
   31	format(1X,I4,1X,A6,1X,I2,3(1X,F6.2),1X,F4.1)
                                      ENDIF			!3
	ENDIF             ! 2
	                ELSE ! 1>>>>>>>>>>>>>>>>>>>>
	if (wkp.ge.zkp) then	!               4
      IF ((wun.le.wutop).and.(wun.ge.wubot).and.
	+	(wln.le.wltop).and.(wln.ge.wlbot)) THEN        ! 5
	 kout= kout+1
	write(*,31) kout,YRMNDY(nn),iut(nn),(WULE(k,nn),k=1,4)
	write(12,31) kout,YRMNDY(nn),iut(nn),(WULE(k,nn),k=1,4)
                        ENDIF								! 5
	endif									  ! 4
	                ENDIF !1<<<<<<<<<<<<<<<<<
   62	CONTINUE
   60	CONTINUE
C
	close(unit=12)
C
       pause ' '
	stop
	end