      program mapprobwfi1s
C.........................................................................June 2021
C maps WGTEC for selected SOLPHASE with GTEC-Shub median for 9 extreme storms SAC23
C    .....................................................................Mar 2021
C Probability GIM-PWf of Wf = -2,-3,-4,+2,+3,+4 for a month
C at each grid on a map
C
      CHARACTER*128 INFILE,OUTFILEP,OUTFILEN,INSTORM
	      CHARACTER*2 AYR,AMN,ADY,AUT,AMM,NSTRM(9)
	CHARACTER*3 PHI
	DIMENSION prob(0:72),iwf(0:72),iprob(0:72)
	+,icp_2(0:72,71),icp_3(0:72,71),icp_4(0:72,71)
	+,icp2(0:72,71),icp3(0:72,71),icp4(0:72,71)
	+,icn_2(0:72,71),icn_3(0:72,71),icn_4(0:72,71)
	+,icn2(0:72,71),icn3(0:72,71),icn4(0:72,71)
C	DATA IM/31,28,31,30,31,30,31,31,30,31,30,31/
	DATA NSTRM/'12','13','15','16','17','19','20','21','22'/
C
      instorm='wulephasNNsm.txt'
C      infile='d:\web\dws\YY\MM\wcMMDDUT.jYY'   ! month for analysys
C       infile='d:\web\graf\dwu\YY\MM\ut\wcMMDDUTMM.uYY'   !UPC month for analysys
       infile='d:\web\gtc\YY\MM\wsMMDDUTMM.sYY'   !JPL-15min Shub median 
	      outfilep='wsmphixxp.txt'
	outfilen='wsmphixxn.txt'
C      Cycle on SOLPHI  Main phase
CREM       phi0=0.0
        phi0=0.9  ! peak
	cmapp=0.0
	cmapn=0.0
CREM      DO 551 jphi=1,10
	phi1=phi0+0.1  !
	phi10=phi1*10
	iphi10=int(phi10)
	call blet3(iphi10,PHI)
	outfilep(4:6)=PHI
	outfilen(4:6)=PHI
C
	do i=1,71
	do k=0,72
	icp_2(k,i)=0
	icp_3(k,i)=0
	icp_4(k,i)=0
	icp2(k,i)=0
	icp3(k,i)=0
	icp4(k,i)=0
	icn_2(k,i)=0
	icn_3(k,i)=0
	icn_4(k,i)=0
	icn2(k,i)=0
	icn3(k,i)=0
	icn4(k,i)=0
	enddo
	enddo
C
   17	format(A2)
   30 FORMAT(1X,I4,2(1X,I2),1X,F5.2,2(1X,F6.2,1X,F6.3))
   13 FORMAT(I2)
   32 FORMAT(73(1X,I4))
C
C Cycle on extreme storms SC23
C
       DO 552 jev=1,9
C	call blet2(jev,NN)
	instorm(9:10)=NSTRM(jev)
	OPEN(12,file=instorm,action='READ')
C
C Select range PHI
C
	DO 553 itime=0,96
	READ(12,30,end=554) iyyyy,imn,idy,RUT,wpos,fipos,wneg,fineg
C Wpos
	IF ((fipos.ge.phi0).and.(fipos.le.phi1)) THEN  !POSSSSSSSSSSSSSSSSSSSS
	cmapp=cmapp+1.0
	iyr=iyyyy-(iyyyy/100)*100
      call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
	iut=int(rut)
      call blet2(iut,AUT)
	xmin=RUT-float(iut)/0.25*15.
	imin=int(xmin)
	call blet2(imin,AMM)

C
      infile(12:13)=AYR
	infile(30:31)=AYR
C
c1333	CONTINUE
C
	infile(15:16)=AMN
	infile(20:21)=AMN
C
       infile(22:23)=ADY
C
	    infile(24:25)=AUT
C
      infile(26:27)=AMM
      OPEN(10,file=infile,action='READ')
C
       do lat=1,71
	  read(10,32,end=35) (iwf(k),k=0,72) 
         do lon=0,72
	if (iwf(lon).eq.-2) icp_2(lon,lat)=icp_2(lon,lat)+1
      if (iwf(lon).eq.-3) icp_3(lon,lat)=icp_3(lon,lat)+1
	if (iwf(lon).eq.-4) icp_4(lon,lat)=icp_4(lon,lat)+1
	if (iwf(lon).eq.2) icp2(lon,lat)=icp2(lon,lat)+1
	if (iwf(lon).eq.3) icp3(lon,lat)=icp3(lon,lat)+1
	if (iwf(lon).eq.4) icp4(lon,lat)=icp4(lon,lat)+1
	   enddo
       enddo
   35 close(unit=10)
                            ENDIF  !POSSSSSSSSSSSSSSS
C Wneg
	IF ((fineg.ge.phi0).and.(fineg.le.phi1)) THEN  !NEGGGGGGGGGGGGGGGGGG
	cmapn=cmapn+1.0
	iyr=iyyyy-(iyyyy/100)*100
      call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
	iut=int(rut)
      call blet2(iut,AUT)
	xmin=RUT-float(iut)/0.25*15.
	imin=int(xmin)
	call blet2(imin,AMM)

C
      infile(12:13)=AYR
	infile(30:31)=AYR
C
c1333	CONTINUE
C
	infile(15:16)=AMN
	infile(20:21)=AMN
C
       infile(22:23)=ADY
C
	    infile(24:25)=AUT
C
      infile(26:27)=AMM
      OPEN(10,file=infile,action='READ')
C
       do lat=1,71
	  read(10,32,end=36) (iwf(k),k=0,72) 
         do lon=0,72
	if (iwf(lon).eq.-2) icn_2(lon,lat)=icn_2(lon,lat)+1
      if (iwf(lon).eq.-3) icn_3(lon,lat)=icn_3(lon,lat)+1
	if (iwf(lon).eq.-4) icn_4(lon,lat)=icn_4(lon,lat)+1
	if (iwf(lon).eq.2) icn2(lon,lat)=icn2(lon,lat)+1
	if (iwf(lon).eq.3) icn3(lon,lat)=icn3(lon,lat)+1
	if (iwf(lon).eq.4) icn4(lon,lat)=icn4(lon,lat)+1
	   enddo
       enddo
   36 close(unit=10)
                            ENDIF  !NEGGGGGGGGGGGGGGGGGGGGGG
  553	CONTINUE       ! end 96 UT+MM times for NN event
  554	CLOSE(unit=12)
C
  552	CONTINUE       ! end 9 storms events
C 
C output results for selected PHI range

C??        cntmon=96.*(lda2-lda1+1.)
C      outfile='mapYYMMD1wsxx.txt'
c==	outfile(4:5)=AYR
c==	outfile(6:7)=AMN
c==      outfile(8:9)=ADY
C POS results
	cntmon=cmapp       !POS
       outfilep(7:8)='_2'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C???      prob(lon)=icn_2(lon,lat)/cntmon*10000000.
      prob(lon)=icp_2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilep(7:8)='_3'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn_3(lon,lat)/cntmon*10000000.
      prob(lon)=icp_3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilep(7:8)='_4'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn_4(lon,lat)/cntmon*10000000.
      prob(lon)=icp_4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilep(7:8)='p2'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn2(lon,lat)/cntmon*10000000.
      prob(lon)=icp2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilep(7:8)='p3'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn3(lon,lat)/cntmon*10000000.
      prob(lon)=icp3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilep(7:8)='p4'
	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn4(lon,lat)/cntmon*10000000.
      prob(lon)=icp4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
C NEG results
	cntmon=cmapn       !NEG
       outfilen(7:8)='_2'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C???      prob(lon)=icn_2(lon,lat)/cntmon*10000000.
      prob(lon)=icn_2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilen(7:8)='_3'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn_3(lon,lat)/cntmon*10000000.
      prob(lon)=icn_3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilen(7:8)='_4'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn_4(lon,lat)/cntmon*10000000.
      prob(lon)=icn_4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilen(7:8)='p2'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn2(lon,lat)/cntmon*10000000.
      prob(lon)=icn2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilen(7:8)='p3'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn3(lon,lat)/cntmon*10000000.
      prob(lon)=icn3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
C
       outfilen(7:8)='p4'
	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
C??      prob(lon)=icn4(lon,lat)/cntmon*10000000.
      prob(lon)=icn4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	  enddo
	write(*,32) (iprob(k),k=0,72)
	write(11,32) (iprob(k),k=0,72)
	ENDDO
	close(unit=11)
  551	CONTINUE
C
	pause ' '
	STOP
	END
