      program mapprobwfi2s
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 = -1, -2,-3,-4,+1, +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),xmap(0:72,71)
	+,icp_1(0:72,71),icp_2(0:72,71),icp_3(0:72,71),icp_4(0:72,71)
	+,icpx1(0:72,71),icpx2(0:72,71),icpx3(0:72,71),icpx4(0:72,71)
	+,icn_1(0:72,71),icn_2(0:72,71),icn_3(0:72,71),icn_4(0:72,71)
	+,icnx1(0:72,71),icnx2(0:72,71),icnx3(0:72,71),icnx4(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.95  ! peak
	cmapp=0.0
	cmapn=0.0
CREM      DO 551 jphi=1,10
	phi1=phi0+0.05  !
	phi3=-1.*phi1		   !TEMP only for storm peak
	phi4=phi3+0.05		   !TEMP only for storm peak
	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_1(k,i)=0
	icp_2(k,i)=0
	icp_3(k,i)=0
	icp_4(k,i)=0
	icpx1(k,i)=0
	icpx2(k,i)=0
	icpx3(k,i)=0
	icpx4(k,i)=0
	icn_1(k,i)=0
	icn_2(k,i)=0
	icn_3(k,i)=0
	icn_4(k,i)=0
	icnx1(k,i)=0
	icnx2(k,i)=0
	icnx3(k,i)=0
	icnx4(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)).or.
	+((fipos.ge.phi3).and.(fipos.le.phi4))) 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.-1) icp_1(lon,lat)=icp_1(lon,lat)+1
	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.1) icpx1(lon,lat)=icpx1(lon,lat)+1
	if (iwf(lon).eq.2) icpx2(lon,lat)=icpx2(lon,lat)+1
	if (iwf(lon).eq.3) icpx3(lon,lat)=icpx3(lon,lat)+1
	if (iwf(lon).eq.4) icpx4(lon,lat)=icpx4(lon,lat)+1
	   enddo
       enddo
   35 close(unit=10)
                            ENDIF  !POSSSSSSSSSSSSSSS
C Wneg
	IF (((fineg.ge.phi0).and.(fineg.le.phi1)).or.
	+((fineg.ge.phi3).and.(fineg.le.phi4))) 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.-1) icn_1(lon,lat)=icn_1(lon,lat)+1
	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.1) icnx1(lon,lat)=icnx1(lon,lat)+1
	if (iwf(lon).eq.2) icnx2(lon,lat)=icnx2(lon,lat)+1
	if (iwf(lon).eq.3) icnx3(lon,lat)=icnx3(lon,lat)+1
	if (iwf(lon).eq.4) icnx4(lon,lat)=icnx4(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 POS results
	cntmon=cmapp       !POS
       outfilep(7:8)='xx'
	OPEN(11,file=outfilep,access='APPEND')
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icp_1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp_1,max_1)
	write(*,*) 'up_1: ',medp_1,max_1
	write(11,*) 'up_1: ',medp_1,max_1
C
       outfilep(7:8)='_2'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icp_2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp_2,max_2)
	write(*,*) 'up_2: ',medp_2,max_2
	write(11,*) 'up_2: ',medp_2,max_2
C
       outfilep(7:8)='_3'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icp_3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp_3,max_3)
	write(*,*) 'up_3: ',medp_3,max_3
	write(11,*) 'up_3: ',medp_3,max_3
C
       outfilep(7:8)='_4'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icp_4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp_4,max_4)
	write(*,*) 'up_4: ',medp_4,max_4
	write(11,*) 'up_4: ',medp_4,max_4
C
       outfilep(7:8)='p1'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icpx1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp11,maxp11)
	write(*,*) 'up11: ',medp11,maxp11
	write(11,*) 'up11: ',medp11,maxp11
C
       outfilep(7:8)='p2'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icpx2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp22,maxp22)
	write(*,*) 'up22: ',medp22,maxp22
	write(11,*) 'up22: ',medp22,maxp22
C
       outfilep(7:8)='p3'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icpx3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp33,maxp33)
	write(*,*) 'up33: ',medp33,maxp33
	write(11,*) 'up33: ',medp33,maxp33
C
       outfilep(7:8)='p4'
c	OPEN(11,file=outfilep)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icpx4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medp44,maxp44)
	write(*,*) 'up44: ',medp44,maxp44
	write(11,*) 'up44: ',medp44,maxp44
C
C NEG results
	cntmon=cmapn       !NEG
       outfilen(7:8)='_1'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icn_1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn_1,maxn_1)
	write(*,*) 'lo_1: ',medn_1,maxn_1
	write(11,*) 'lo_1: ',medn_1,maxn_1
C
       outfilen(7:8)='_2'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icn_2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn_2,maxn_2)
	write(*,*) 'lo_2: ',medn_2,maxn_2
	write(11,*) 'lo_2: ',medn_2,maxn_2
C
       outfilen(7:8)='_3'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icn_3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn_3,maxn_3)
	write(*,*) 'lo_3: ',medn_3,maxn_3
	write(11,*) 'lo_3: ',medn_3,maxn_3
C
       outfilen(7:8)='_4'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icn_4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn_4,maxn_4)
	write(*,*) 'lo_4: ',medn_4,maxn_4
	write(11,*) 'lo_4: ',medn_4,maxn_4
C
       outfilen(7:8)='p1'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icnx1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn11,maxn11)
	write(*,*) 'lo11: ',medn11,maxn11
	write(11,*) 'lo11: ',medn11,maxn11
C
       outfilen(7:8)='p2'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icnx2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn22,maxn22)
	write(*,*) 'lo22: ',medn22,maxn22
	write(11,*) 'lo22: ',medn22,maxn22
C
       outfilen(7:8)='p3'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icnx3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
c	close(unit=11)
	call median5183(xmap,medn33,maxn33)
	write(*,*) 'lo33: ',medn33,maxn33
	write(11,*) 'lo33: ',medn33,maxn33
C
       outfilen(7:8)='p4'
c	OPEN(11,file=outfilen)
	DO lat=1,71
	  do lon=0,72
      prob(lon)=icnx4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	xmap(lon,lat)=prob(lon)
	  enddo
c	write(*,32) (iprob(k),k=0,72)
c	write(11,32) (iprob(k),k=0,72)
	ENDDO
  	call median5183(xmap,medn44,maxn44)
	write(*,*) 'lo44: ',medn44,maxn44
	write(11,*) 'lo44: ',medn44,maxn44
	close(unit=11)
551	CONTINUE
C
	pause ' '
	STOP
	END
