	program probhwpjc
C.....................................................Mar. 2013
C Probability of occurrence of W+ , W- storm
C produce counts of CWp=3 & 4, /5112 grids on map * 100% (Occur.%)
C produce counts of CWn=-3 & -4, /5112 grids on map * 100% (Occur.%)
C CWt=cwp+cwn
C
C.....................................................Nov. 2012
C JPL maps
C.....................................................Apr. 2012
C
C IONEX maps
C for glat=-87.5:2.5:87.5, glong=-180:5:165  (5112 => 840 grids)
C Correction of (W+max-W-max)*(1+cntstrm/5112)
C.....................................................Feb. 2012
C
C.....................................................Aug.2011
C
C Input: maps of W-index: tcMNDYUT.wYR files 
C             for glat=-85:5:85, glong=-180:15:165  (=840 grids)
C Correction of (W+max-W-max)*(1+cntstrm/840)
C
C........................................................Apr.2008
C...... Count RANGE only when W+=>3, and/or W-<=-3 ..............
C
c........................................................Oct.2007
C Correction of (W+max-W-max)*(1+cntstrm/600)
C cntstrm=total number of W+=3,W+=4, W-=-3,W-=-4 for the IONEX-map
C
c .......................................................Sep. 2007
C Range of TEC indices from 35 latitude results TEC.fYR.........
C
c from hourly IONEX maps
c
	DIMENSION IRESP(0:71),IRESN(0:71),DHR(0:71),IM(12)
C 
     +,IRES(0:71)
	 INTEGER*2 IYR,IMN,IDY,idy1,idy2,iut
	CHARACTER(10) DD
      CHARACTER(10) TT
      CHARACTER(5) ZZ
	CHARACTER*64,infile,outfile,outfilep,outfilen,title
	CHARACTER*133 wpast(0:365),wwpos(0:365),wwneg(0:365),txt
	CHARACTER*2 AYR,AMN,AUT,ADY
	+,aa(24),bb
	CHARACTER*6 YRMNDY
	CHARACTER*4 AYEAR,AMEAN
      DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	DATA aa/'00','01','02','03','04','05','06','07','08','09','10'
	+,'11','12','13','14','15','16','17','18','19','20','21','22','23'/
	bb='UT'
	amean='mean'
	title='YYMMDD UT CNTNEG CNTPOS CNTTOT   TOT%    Wp'
	CALL DATE_AND_TIME(DATE=DD,TIME=TT,ZONE=ZZ)
      TT(5:)='      '
      WRITE(*,*) 'PC Date: Year,Month,Day = ',DD,'Time = ',TT

C start
      WRITE(*,'(\A\)') ' ENTER YEAR OF INPUT FILE='
      READ(*,'(A)') AYEAR
CREM	AYEAR='2012'
	ayr=ayear(3:4)
	read(ayr,*) year
	iyr=int(year)
C
      WRITE(*,'(\A\)') ' ENTER MONTH OF INPUT FILE='
      READ(*,'(A)') AMN
	read(amn,*) rmn
	mn=int(rmn)
	if(int(year/4.)*4.eq.iyr) THEN 
	ndnr=366
	im(2)=29
	                          ELSE
	ndnr=365
	im(2)=28
	ENDIF
C Cycle on month
c=	DO 100 imn=1,mn
	imn=mn
	call blet2(imn,AMN)

c     WRITE(*,'(\A\)') ' ENTER DY1 & DY2 ='
c     READ(*,13) idy1,idy2
	idy1=1
	idy2=im(imn)
 13   FORMAT(I2)

CC CALCULATION OF DAYS OF YEAR	
C?      infile='c:\web\dws\YR\MN\wcMNDYUT.jYR' ! JPL W(glat & glong map)
      infile='c:\web\dws\YR\MN\wiMNDYUT.jYR' ! JPL W(mlat & mlong map) -15d
  	infile(12:13)=ayr
	infile(15:16)=amn
  	infile(28:29)=ayr
	infile(20:21)=amn
      outfile='wjticYR.txt'	!W-indices from W-maps geogr.coord.<oc>=occurrence
	outfile(6:7)=AYR
 	outfilep=outfile
	outfilen=outfile
	outfilep(1:1)='p'
	outfilen(1:1)='n'
   14	format(1X,A43)
C
  	OPEN(112,FILE=OUTFILE)
	OPEN(113,FILE=OUTFILEP)
	OPEN(114,FILE=OUTFILEN)
 198	format(4X,A2,24(3X,A2),3X,A4)
C
	lda1=ndoy(iyr,imn,idy1)  ! day-of-year
	mndy=im(imn)
  130	format(A133)

	if((imn.eq.1).and.(idy1.eq.1)) then 

	write(112,198) bb,aa,amean
	write(113,198) bb,aa,amean
	 write(114,198) bb,aa,amean
      write(*,198) bb,aa,amean
	goto 197 ! avoid read-re-write
	endif
C
C read-re-write outfile for days from 1 to lda1-1:
	ncnt=0
	do 149 ipast=0,lda1-1
	read(112,130,END=151,ERR=2) txt
	wpast(ipast)=txt
	YRMNDY=txt(1:6)
  151	read(113,130,END=152,ERR=2) txt
 	wwpos(ipast)=txt
  152	read(114,130,END=153,ERR=2) txt
 	wwneg(ipast)=txt
	if (yrmndy.eq.'      ') goto 149
	ncnt=ncnt+1
  149	continue
  153 close(unit=112)
      close(unit=113)
	close(unit=114)
C
  	OPEN(112,FILE=OUTFILE)
	OPEN(113,FILE=OUTFILEP)
	OPEN(114,FILE=OUTFILEN)
C re-write past data:
	do ipast=0,ncnt-1
	write(112,130) wpast(ipast)
cc	write(*,130) wpast(ipast)  ! test
	write(113,130) wwpos(ipast)
	write(114,130) wwneg(ipast)
	enddo
      close(unit=112)
      close(unit=113)
	close(unit=114)
CC  
  197 continue
  	OPEN(112,FILE=OUTFILE,ACCESS='APPEND')
	OPEN(113,FILE=OUTFILEP,ACCESS='APPEND')
	OPEN(114,FILE=OUTFILEN,ACCESS='APPEND')

C
C  day-to-day cycle
C ========================================

       DO 500 idy=idy1,idy2
	lda=ndoy(iyr,imn,idy)  ! day-of-year
	call blet2(idy,ady)
	infile(22:23)=ady
C
      do k=0,23
	ires(k)=0
	iresp(k)=0
	iresn(k)=0
	enddo
C
 	sum=0.
	sumn=0.
	sump=0.

C cycle on UT:
C -----------------------------------------------------
      DO 801 iut=0,23		! cycle on UT
	call blet2(iut,aut)
	infile(24:25)=aut
	OPEN(111,FILE=INFILE)
c
	cnttot=0.       !!!new for W-map
	cntneg=0.       !!!new for W-map
	cntpos=0.       !!!new for W-map
C ===== Input datafile 
 1311	DO 1799 jj=1,71	  ! glats= -87.5:2.5:87.5
c	pause ' '
    1	READ (111,12,END=3,ERR=2) (DHR(k),k=0,71) ! W(longs=-180:5:175)
  12	FORMAT(72(1X,F4.0))
C
	DO k=0,71				! Cycle on Longi	 <<<<<<<<
C++	Count stormy indices W-:
	if (dhr(k).lt.-2) then 
      cntneg=cntneg+1.  !!!new for W map
	cnttot=cnttot+1.  !!!new for W map
	 endif
C++	Count stormy indices W+:
		if (dhr(k).gt.2) then 
      cnttot=cnttot+1.  !!!new for W map
	cntpos=cntpos+1.  !!!new for W map
	 endif

	enddo		    ! Cycle on Longi	 >>>>>>>>>>>>>


 1799	CONTINUE				   ! Cycle jj on glats (total W map)
  3   CLOSE(unit=111)

C----------------------------------------------------------
	cres=cnttot/0.5112
	cresp=cntpos/0.5112
	cresn=cntneg/0.5112
	ires(iut)=nint(cres) !prob*10000. = occur*100,%
	iresn(iut)=nint(cresn)
	iresp(iut)=nint(cresp)
c   
  801 continue            ! end cycle on UT 

C  Daily mean:
C
	sum=0.
	sump=0.
	sumn=0.
		do m=0,23      ! cycle on UT		======
 175	sum=sum+float(ires(m))
	sump=sump+float(iresp(m))
	sumn=sumn+float(iresn(m))
	enddo							 !=========
	isum=nint(sum/24.)
	isump=nint(sump/24.)
	isumn=nint(sumn/24.)
	yrmndy='000000'
	yrmndy(1:2)=AYR
	yrmndy(3:4)=AMN
	yrmndy(5:6)=ADY	
	write(112,72) yrmndy,(ires(k),k=0,23),isum
		write(113,72) yrmndy,(iresp(k),k=0,23),isump
		write(114,72) yrmndy,(iresn(k),k=0,23),isumn
	write(*,72) yrmndy,(ires(k),k=0,23),isum
 	write(*,72) yrmndy,(iresp(k),k=0,23),isump
 	write(*,72) yrmndy,(iresn(k),k=0,23),isumn
C
  500	CONTINUE ! next day
C
  100	CONTINUE ! next	month
   72 format(A6,24(1X,I4),1X,I6)
	CLOSE(unit=112)
	CLOSE(unit=113)
	CLOSE(unit=114)
	goto 900
    2 CONTINUE
           write(*,*) 'INPUT FILE IS NOT IN YOUR DIRECTORY '
  900 pause ' '
      STOP
      END

