C.........................................................Aug. 2021
C
c VERIFICATION OF MEDIAN ================
C
      subroutine smed1std3(ndnr,CC,SMED,STD,STD1,STD2)
C
C Median smed, and standard deviation std, std1(neg),std2(pos) for -14+1 days at specified hour
C
C Input file: CC(ndnr) for ndnr data (1 col), 
C ncol=72 for longitudes -180:5:180 or ncol=23 for hours UT=0,1,...,23
C Number days for median : ndnr=15	(14 preceding days+current day)
C Output arrays: SMED(0:72) for longitudes -180:5:180 (or SMED(0:23) for UT hours)
C                 STD1(0:72),STD2(0:72)
C
	DIMENSION CC(2000)
C	+,SMED,STD,STD1,STD2
	+,X1(2000),X2(2000)
C "MEDIAN:"
C       DO 1030 M=0,ncol
      MM=0
C Add calculation of average:
	aver=0.
      DO N=1,ndnr
      X1(N)=CC(N)
	aver=aver+x1(n)
      ENDDO 
	aver=aver/ndnr
  880 XX=0.
CREM MAX X1(N)
  890 DO N=1,ndnr
  900 IF (XX.lt.X1(N)) THEN 
      XX=X1(N)
	ENDIF
      ENDDO      !N
  920 IF (XX.eq.0) GOTO 980
C
  930 DO 960 N=1,ndnr
      IF (X1(N).lt.XX) GOTO 960
      MM=MM+1 
      X2(MM)=XX 
	X1(N)=0.
  960 CONTINUE      ! N
  970 IF (MM.lt.NDNR) GOTO 880
  980 XMM=FLOAT(MM)
      XD1=XMM/2.
      ND1=INT(XD1)
	CD1=FLOAT(ND1) 
      ND2=ND1+1 
      IF (CD1.LT.XD1) THEN 
	XMD=X2(ND2)
	else
	XMD=(X2(ND1)+X2(ND2))/2.
	ENDIF
C
      SMED=XMD
C Standard deviation:
	sumneg=0.
	sumpos=0.
	icntpos=0
	icntneg=0
	avneg=0.
	avpos=0.
	iavneg=0
	iavpos=0
	do ii=1,ndnr
	if (x2(ii).gt.xmd) then
	sumpos=sumpos+(x2(ii)-xmd)**2
	icntpos=icntpos+1
	endif
	if (x2(ii).lt.xmd) then
	sumneg=sumneg+(x2(ii)-xmd)**2
	icntneg=icntneg+1
	endif
	if (x2(ii).gt.aver) then
	avpos=avpos+(x2(ii)-aver)**2
	iavpos=iavpos+1
	endif
	if (x2(ii).le.aver) then
	avneg=avneg+(x2(ii)-aver)**2
	iavneg=iavneg+1
	endif
	enddo
	STD=0.0
	STD1=0.0
	STD2=0.0
	if ((icntneg+icntpos).gt.0) then
	STD=sqrt((sumneg+sumpos)/(icntneg+icntpos))
	else
	if((iavneg+iavpos).gt.0) STD=sqrt((avneg+avpos)/
	&(iavneg+iavpos))
	endif
	if (icntneg.gt.0) then 
	if (iavneg.gt.0) STD1=sqrt(sumneg/icntneg)
	else
	if (iavneg.gt.0) STD1=sqrt(avneg/iavneg)
	endif
	if (icntpos.gt.0) then 
	STD2=sqrt(sumpos/icntpos)
	else
	if (iavpos.gt.0) STD2=sqrt(avpos/iavpos)
	endif
	if (std1.eq.0.0) STD1=STD
	if (std2.eq.0.0) STD2=STD
 1030 CONTINUE
C 
c
 1050 RETURN
      END