      program mapprobwfi8s1
C.......................................................................1 July 2021
C
C .....................................................................27 June 2021 
C Dominant Pos/Neg before Abs(Fipos) > Abs(Fi neg)
C and Vice Versa: Dominant Neg/Pos after Abs(Fipos) > Abs(Fi neg)
C.........................................................................June 2021
C
C Combine W+- levels results ito Proxy GIM-W-index map
C Produce W+- levels results for Wup and Wlo
C Add Driving Phase Wup Wlo for testing or modelling set
C
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,INPHASE,OUTFILE,INSTORM,OUTPROB ! prob_med, prob_max
     +,outmap
	 CHARACTER*2 AYR,AMN,ADY,AUT,AMM,NSTRM(9),ANN
	+,BYR,BMN,BDY,BUT
	DIMENSION prob(0:72),iwf(0:72),iprob(0:72),xmap(0:72,71)
	+,imapwu(0:72,71),imapwl(0:72,71)
     +,irtot(8),imtot(8),iwout(8),imaxpn(8) !imaxpn: +1,+2,+3,+4, -1,-2,-3,-4;
     +,idat(9,3,0:192),rdat(9,2,0:192) !idat(iyr,imn,idy),rdat(rut,wuphi,wlphi)
	+,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)
	+,icp01(0:72,71),icp02(0:72,71),icp03(0:72,71),icp04(0:72,71)
	+,icp11(0:72,71),icp22(0:72,71),icp33(0:72,71),icp44(0:72,71)
	+,iwc(8)
	DATA NSTRM/'12','13','15','16','17','19','20','21','22'/
	DATA IWC/1,2,3,4,-1,-2,-3,-4/
C
      instorm='wulephasNNsm.txt'
      inphase='wuwlphasMM.txt' !Wu/WL phase for 3 days => output proxy GIM-Wx index
       infile='d:\web\gtc\YY\MM\wsMMDDUTMM.sYY'   !JPL-15min Shub median 
	 outfile='d:\web\graf\duc\YR\MN\ut\wxMNDDUT.sYR' !combined proxy W-index 1h step
	outprob='mapfiprobNN.txt'
	outmap='wxYRMNDDUT.txt'
C
   37 format(1X,I4,2I2,1X,I2,2(1X,F6.3))

        write(*,*) ' Enter ANN for output GIM-W: '
	read(*,*) ANN
CTemp
      write(*,*) ' Enter No.NSTRM (0, 1...9) to exclude from results:'
	read(*,*) nexclude   
	inphase(9:10)=ANN
	outprob(10:11)=ANN
	OPEN(12,file=outprob,access='APPEND')
	write(*,*) 'Exclude: ',nexclude
	write(12,*) 'Exclude: ',nexclude
	OPEN(13,file=inphase,action='READ')
C
C input training files:15 min
C
       DO 555 jev=1,9   
	instorm(9:10)=NSTRM(jev)
	OPEN(14,file=instorm,action='READ')
	DO ii=0,192
	READ(14,30,end=554) idat(jev,1,ii),idat(jev,2,ii),idat(jev,3,ii)
	+,rdat(jev,1,ii),wpos,rdat(jev,2,ii),wneg,rdat(jev,3,ii)
	 enddo
  554	CLOSE(unit=14)
  555	CONTINUE
C Cycle on 72 1h GIM-Wu/Wl maps:
      phiup_pre=0.1
	philo_pre=0.
C   
       DO 121 ihr=0,71
       read(13,37,end=49) kyear,kmn,kdy,kut,phiwup,phiwlo  ! for output files
	kyr=kyear-(kyear/100)*100
      call blet2(kyr,BYR)
C-	AYR=BYR
	call blet2(kmn,BMN)
	call blet2(kdy,BDY)
      call blet2(kut,BUT)
        write(*,37) kyear,kmn,kdy,kut,phiwup,phiwlo
	  write(12,37) kyear,kmn,kdy,kut,phiwup,phiwlo
c	 outfilep='d:\web\graf\duc\YR\MN\ut\wxMNDDUT.sYR' !combined proxy W-index 1h step
	outfile(17:18)=BYR
	outfile(36:37)=BYR
	outfile(20:21)=BMN
	outfile(28:29)=BMN
	outfile(30:31)=BDY
	outfile(32:33)=BUT
C
       if (phiwup.ge.0.) iopt=1        ! Option for select priority WU
       if (phiwlo.lt.0.) iopt=2        ! Option for select priority WL
	IF ((phiwup.lt.0.).and.(phiwlo.gt.0.)) THEN
	if (abs(phiwup).lt.abs(phiwlo)) then 
	      iopt=2
	else
            iopt=1
	endif
	ENDIF
C      Cycle on SOLPHI  Main phase
C WUP
        IF (phiwup.le.1.) THEN	  !phi positive.
        iflag_1=0
        phi0up=phiwup-0.025  ! 
	if (phiwup.eq.0.) phi0up=phiwup
	phi1up=phi0up+0.025  !
	if (phi1up.gt.1.) phi1up=1.
	                  ENDIF
	  IF (phiwup.lt.0.) THEN	   !phi negative.
        iflag_1=1
	phi0up=phiwup-0.025
	if (phiwup.lt.-1.0) phi0up=-1.
	phi1up=phi0up-0.025		   
      if (ihr.eq.71) phi1up=phiwup
      endif
C WLO
        IF (phiwlo.le.1.) THEN	  !phi positive.
        iflag_2=0
        phi0l=phiwlo-0.025  ! 
	if (phiwlo.eq.0.) phi0l=phiwlo
	phi1l=phi0l+0.025  !
	if (phi1l.gt.1.) phi1lo=1.
	                  ENDIF
	  IF (phiwlo.lt.0.) THEN	   !phi negative.
        iflag_2=1
	phi0lo=phiwlo-0.025
	if (phiwlo.lt.-1.0) phi0lo=-1.
	phi1lo=phi0lo-0.025		   
      if (ihr.eq.71) phi1lo=phiwlo
      endif
	cmapp=0.0
c-	cmapn=0.0
CREM      DO 551 jphi=1,10
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
	imapwu(k,i)=0
	imapwl(k,i)=0	   !C-?
	enddo
	enddo
C
   17	format(A2)
   30 FORMAT(3X,I2,2(1X,I2),1X,F5.2,2(1X,F6.2,1X,F6.3))
   13 FORMAT(I2)
   32 FORMAT(73(1X,I4))
C
C
C
C Cycle on extreme storms SC23
C
	cmapp=0.
       DO 552 jev=1,9			! All training storms
	if (jev.eq.nexclude) goto 552      ! => exclude one training case
C
C Select range PHI
C
        itup=0
	itlo=0
C
C Wpos     POSSSSSSSSSSSSSSSSSSSS
        DO 1331 itm=1,192	! 
C Temp exclude 1 training set:
      phi_pre=rdat(jev,2,itm-1)
	phi_cur=rdat(jev,2,itm)
	phi_nex=rdat(jev,2,itm+1)
C
	if ((phi_pre.eq.0.).and.(phiwup.eq.0)) then
	itup=itm-1
	GOTO 1332
	endif
C
	if ((phi_cur.eq.1.).and.(phiwup.eq.1)) then
	itup=itm
	GOTO 1332
	endif
C
      IF ((phiwup.ge.phi_pre).and.(phiwup.le.phi_cur)) THEN
	d_pre=phiwup-phi_pre
 	d_nex=phi_nex-phiwup
	if (d_pre.lt.d_nex) then
	itup=itm-1
	                 else
	itup=itm
	                 endif
	                                         ELSE
	GOTO 1331
	                                   ENDIF
 1332	cmapp=cmapp+1.0
      iyr=idat(jev,1,itup)
	call blet2(iyr,AYR)
	imn=idat(jev,2,itup)
	idy=idat(jev,3,itup)
	rut=rdat(jev,1,itup)
	iut=int(rut)
	xmin=RUT-float(iut)/0.25*15.
	imin=int(xmin)
C
      if (phiwup.eq.0.) then
	lda=ndoy(iyr,imn,idy)
	lda_pre=lda-1
	call submmdd(lda_pre,iyr,imn,idy)
	iut=23
	imin=45
	endif
	call blet2(imn,AMN)
	call blet2(idy,ADY)
      call blet2(iut,AUT)
	call blet2(imin,AMM)
C
      infile(12:13)=AYR
	infile(30:31)=AYR
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)
      EXIT
 1331 CONTINUE  !POSSSSSSSSSSSSSSS
C Wneg
        DO 1333 itm=1,192	   !NEGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
      phi_pre=rdat(jev,2,itm-1)
	phi_cur=rdat(jev,2,itm)
	phi_nex=rdat(jev,2,itm+1)
C
	if ((phi_pre.eq.0.).and.(phiwlo.eq.0)) then
	itlo=itm-1
	GOTO 1334
	endif
C
	if ((phi_cur.eq.1.).and.(phiwlo.eq.1)) then
	itlo=itm
	GOTO 1334
	endif
C
      IF ((phiwlo.ge.phi_pre).and.(phiwlo.le.phi_cur)) THEN
	d_pre=phiwlo-phi_pre
 	d_nex=phi_nex-phiwlo
	if (d_pre.lt.d_nex) then
	itlo=itm-1
	                 else
	itlo=itm
	                 endif
	                                         ELSE
	GOTO 1333
	                                   ENDIF
 1334	cmapp=cmapp+1.0
      iyr=idat(jev,1,itlo)
	imn=idat(jev,2,itlo)
	idy=idat(jev,3,itlo)
	rut=rdat(jev,1,itlo)
	iut=int(rut)
	xmin=RUT-float(iut)/0.25*15.
	imin=int(xmin)
C
      if (phiwlo.eq.0.) then
	lda=ndoy(iyr,imn,idy)
	lda_pre=lda-1
	call submmdd(lda_pre,iyr,imn,idy)
	iut=23
	imin=45
	endif
	call blet2(iyr,AYR)
	call blet2(imn,AMN)
	call blet2(idy,ADY)
      call blet2(iut,AUT)
	call blet2(imin,AMM)
C
      infile(12:13)=AYR
	infile(30:31)=AYR
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) 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
   36 close(unit=10)
      EXIT
 1333 CONTINUE  !NEGGGGGGGGGGGGGGGGGGGG
C
  552	CONTINUE       ! end 9 storms events
C 
C output results for selected PHI range
C
C =================================================
C POS + NEG results
   38 FORMAT(1X,I2,5(1X,I5))
c      iw=0
C
	cntmon=cmapp       !POS					   W=-1
C       outfilep(7:8)='xx' ! med, max
C       outfilep(7:8)='_1'
	DO lat=1,71			!POS					   W=-1
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icp_1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icp_1(lon,lat)=iprob(lon)  ! Prob
	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,mqua_1,medp_1,max_1)
	imaxpn(5)=max_1
	iw01=-1
c+	write(*,38) iw01,medp_1,max_1
c+	write(12,38) iw,medp_1,max_1
C
	DO lat=1,71    !POS					   W=-2
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icp_2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icp_2(lon,lat)=iprob(lon)  ! Prob
	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,mqua_2,medp_2,max_2)
	imaxpn(6)=max_2
	iw02=-2
c+	write(*,38) iw02,medp_2,max_2
c+	write(12,38) iw,medp_2,max_2
C
	DO lat=1,71				   !POS					   W=-3
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icp_3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icp_3(lon,lat)=iprob(lon)  ! Prob
	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,mqua_3,medp_3,max_3)
	imaxpn(7)=max_3
	iw03=-3
c+	write(*,38) iw03,medp_3,max_3
c+	write(12,38) iw,medp_3,max_3
C
	DO lat=1,71					 !POS					   W=-4
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icp_4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icp_4(lon,lat)=iprob(lon)  ! Prob
	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,mqua_4,medp_4,max_4)
	imaxpn(8)=max_4
	iw04=-4
c+	write(*,38) iw04,medp_4,max_4
c+	write(12,38) iw,medp_4,max_4
C
	DO lat=1,71			   !POS+NEG					   W=+1
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icpx1(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icpx1(lon,lat)=iprob(lon)  ! Prob
	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,mqua11,medp11,maxp11)
	imaxpn(1)=maxp11
	iw11=1
c+	write(*,38) iw,medp11,maxp11
c+	write(12,38) iw,medp11,maxp11
C
	DO lat=1,71				  !POS+neg					   W=+2
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icpx2(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icpx2(lon,lat)=iprob(lon)  ! Prob
	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,mqua22,medp22,maxp22)
	imaxpn(2)=maxp22
	iw22=2
C+	write(*,38) iw,medp22,maxp22
c+	write(12,38) iw,medp22,maxp22
C
	DO lat=1,71				   !POS+NEG					   W=+3
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icpx3(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icpx3(lon,lat)=iprob(lon)  ! Prob
	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,mqua33,medp33,maxp33)
	imaxpn(3)=maxp33
       iw33=3
C+	write(*,38) iw,medp33,maxp33
c+	write(12,38) iw,medp33,maxp33
C
	DO lat=1,71				  !POS+NEG					   W=+4
	  do lon=0,72
	prob(lon)=0.
      if (cntmon.gt.0.) prob(lon)=icpx4(lon,lat)/cntmon*10000.
	iprob(lon)=nint(prob(lon))
	icpx4(lon,lat)=iprob(lon)  ! Prob
	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,mqua44,medp44,maxp44)
	imaxpn(4)=maxp44
	iw44=4
c+	write(*,38) iw,medp44,maxp44
c+	write(12,38) iw,medp44,maxp44
C
C+++++++++++++++++++++++++++++++++++++++++
C
       do lat=1,71
	   do lon=0,72
      kprob=icp_1(lon,lat)
	coef=(mqua_1-medp_1)*(1.0-abs(phiwlo))/2.  !+
	nqua_1=nint(coef+mqua_1)				!+
C+	if ((kprob.ge.medp_1).and.(kprob.gt.icpx1(lon,lat))) then
C=	if (kprob.ge.mqua_1) then
	if (kprob.ge.nqua_1) then			    !+
c-      if (kprob.gt.0) then
	iprob(lon)=-1
	                    else
	iprob(lon)=0
	endif
	icp01(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icp_2(lon,lat)
	coef=(mqua_2-medp_2)*(1.0-abs(phiwlo))/2.  !+
	nqua_2=nint(coef+mqua_2)				!+
C+	if ((kprob.ge.medp_2).and.(kprob.gt.icpx2(lon,lat))) then
C+	if (kprob.ge.mqua_2) then
	if (kprob.ge.nqua_2) then               !+
c-      if (kprob.gt.0) then
	iprob(lon)=-2
	                    else
	iprob(lon)=0
	endif
	icp02(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icp_3(lon,lat)
	coef=(mqua_3-medp_3)*(1.0-abs(phiwlo))/2.  !+
	nqua_3=nint(coef+mqua_3)				!+
C+	if ((kprob.ge.medp_3).and.(kprob.gt.icpx3(lon,lat))) then
C+	if (kprob.ge.mqua_3) then
	if (kprob.ge.nqua_3) then               !+
c-      if (kprob.gt.0) then
	iprob(lon)=-3
	                    else
	iprob(lon)=0
	endif
	icp03(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icp_4(lon,lat)
	coef=(mqua_4-medp_4)*(1.0-abs(phiwlo))/2.  !+
	nqua_4=nint(coef+mqua_4)				!+
C	if ((kprob.ge.medp_4).and.(kprob.gt.icpx4(lon,lat))) then
C+	if (kprob.ge.mqua_4) then
	if (kprob.ge.nqua_4) then                 !+
c-      if (kprob.gt.0) then
	iprob(lon)=-4
	                    else
      iprob(lon)=0
	endif
	icp04(lon,lat)=iprob(lon)
	    enddo
	    enddo
C
         do lat=1,71
	   do lon=0,72
      kprob=icpx1(lon,lat)
	coef=(mqua11-medp11)*(1.0-abs(phiwup))/2.  !+
	nqua11=nint(coef+mqua11)				!+
C+	if ((kprob.ge.medp11).and.(kprob.gt.icp_1(lon,lat))) then
C+	if (kprob.ge.mqua11) then
	if (kprob.ge.nqua11) then               !+
c-      if (kprob.gt.0) then
	iprob(lon)=1
	                    else
	iprob(lon)=0
	endif
	icp11(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icpx2(lon,lat)
	coef=(mqua22-medp22)*(1.0-abs(phiwup))/2.  !+
	nqua22=nint(coef+mqua22)				!+
C+	if ((kprob.ge.medp22).and.(kprob.gt.icp_2(lon,lat))) then
C+	if (kprob.ge.mqua22) then
	if (kprob.ge.nqua22) then               !+
c      if (kprob.gt.0) then
	iprob(lon)=2
	                    else
	iprob(lon)=0
	endif
	icp22(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icpx3(lon,lat)
	coef=(mqua33-medp33)*(1.0-abs(phiwup))/2.  !+
	nqua33=nint(coef+mqua33)				!+
C+	if ((kprob.ge.medp33).and.(kprob.gt.icp_3(lon,lat))) then
c+	if (kprob.ge.mqua33) then
	if (kprob.ge.nqua33) then               !+
c-      if (kprob.gt.0) then
	iprob(lon)=3
	                    else
	iprob(lon)=0
	endif
	icp33(lon,lat)=iprob(lon)
	    enddo
	    enddo
       do lat=1,71
	   do lon=0,72
      kprob=icpx4(lon,lat)
	coef=(mqua44-medp44)*(1.0-abs(phiwup))/2.  !+
	nqua44=nint(coef+mqua44)				!+
C+	if ((kprob.ge.medp44).and.(kprob.gt.icp_4(lon,lat))) then
C+	if (kprob.ge.mqua44) then
	if (kprob.ge.nqua44) then                !+
c      if (kprob.gt.0) then
	iprob(lon)=4
	                    else
	iprob(lon)=0
	endif
	icp44(lon,lat)=iprob(lon)
	    enddo
	    enddo
C Combine GIM-W-index Proxy +++++++++++++++++++++++
C
C Rank output W indices
C
        ip1=0
	ind1=0
        ip5=0
Crem	ind5=0
	ind5=9999				 ! Select min WLx
C 1st max rank Pos/Neg AND 5th max rank Neg/Pos:
C=      IF((phiup_pre.gt.philo_pre).and.(abs(phiwup).gt.abs(phiwlo))
C=	+.or.(phiwup.eq.0.)) THEN
      IF(iopt.eq.1) THEN
C        max => ind1
	do ip=1,4
	if (imaxpn(ip).gt.ind1) then
	ind1=imaxpn(ip)
	ip1=ip  
	endif
	enddo
C
        irtot(1)=ip1
	imtot(1)=ind1
	imaxpn(ip1)=0
	iwout(1)=iwc(ip1)
C			 Max => ind5  <= min 
	do ip=5,8
Crem	if (imaxpn(ip).gt.ind5) then
	if (imaxpn(ip).lt.ind5) then
	ind5=imaxpn(ip)	! min WLx
	ip5=ip  
	endif
	enddo
C
        irtot(5)=ip5
	imtot(5)=ind5
	imaxpn(ip5)=0
	iwout(5)=iwc(ip5)
C
                                  ELSE
C 1st max rank Neg/Pos AND 5th max rank Pos/Neg:
C iopt=2

	do ip=5,8
Crem	if (imaxpn(ip).gt.ind1) then
	if (imaxpn(ip).lt.ind1) then
	ind1=imaxpn(ip)
	ip1=ip  
	endif
	enddo
C
        irtot(1)=ip1
	imtot(1)=ind1
	imaxpn(ip1)=0
	iwout(1)=iwc(ip1)
C
	do ip=1,4
	if (imaxpn(ip).gt.ind5) then
	ind5=imaxpn(ip)
	ip5=ip  
	endif
	enddo
C
        irtot(5)=ip5
	imtot(5)=ind5
	imaxpn(ip5)=0
	iwout(5)=iwc(ip5)
                                   ENDIF								 
C 2nd rank:max rank Pos/Neg AND 6th max rank Neg/Pos:
        ip2=0
	ind2=0
        ip6=0
	ind6=0
c=      IF((phiup_pre.gt.philo_pre).and.(abs(phiwup).gt.abs(phiwlo))
c=	+.or.(phiwup.eq.0.)) THEN
       IF (iopt.eq.1) THEN
C        max => ind2
	do ip=1,4
	if (imaxpn(ip).gt.ind2) then
	ind2=imaxpn(ip)
	ip2=ip  
	endif
   	enddo
C
        irtot(2)=ip2
	imtot(2)=ind2
	imaxpn(ip2)=0
	iwout(2)=iwc(ip2)
C 		max => ind6
	do ip=5,8
	if (imaxpn(ip).gt.ind6) then
	ind6=imaxpn(ip)
	ip6=ip  
	endif
   	enddo
C
        irtot(6)=ip6
	imtot(6)=ind6
	imaxpn(ip6)=0
	iwout(6)=iwc(ip6)
	                                       ENDIF
C 3d rank: :max rank Pos/Neg AND 7th max rank Neg/Pos:
C iopt=2
        ip3=0
	ind3=0
        ip7=0
	ind7=0
c=      IF((phiup_pre.gt.philo_pre).and.(abs(phiwup).gt.abs(phiwlo))
c=	+.or.(phiwup.eq.0.)) THEN
         IF (iopt.eq.1) THEN
C  max => ind3
	do ip=1,4
	if (imaxpn(ip).gt.ind3) then
	ind3=imaxpn(ip)
	ip3=ip  
	endif
   	enddo
C
        irtot(3)=ip3
	imtot(3)=ind3
	imaxpn(ip3)=0
	iwout(3)=iwc(ip3)
C max => 	ind7
	do ip=5,8
	if (imaxpn(ip).gt.ind7) then
	ind7=imaxpn(ip)
	ip7=ip  
	endif
   	enddo
C
        irtot(7)=ip7
	imtot(7)=ind7
	imaxpn(ip7)=0
	iwout(7)=iwc(ip7)
	                                                ELSE
C  max => ind3
C iopt=2
	do ip=5,8
	if (imaxpn(ip).gt.ind3) then
	ind3=imaxpn(ip)
	ip3=ip  
	endif
   	enddo
C
        irtot(3)=ip3
	imtot(3)=ind3
	imaxpn(ip3)=0
	iwout(3)=iwc(ip3)
C  max => ind7
	do ip=1,4
	if (imaxpn(ip).gt.ind7) then
	ind7=imaxpn(ip)
	ip7=ip  
	endif
   	enddo
C
        irtot(7)=ip7
	imtot(7)=ind7
	imaxpn(ip7)=0
	iwout(7)=iwc(ip7)
	                                         ENDIF
C 4th rank: Pos/Neg AND 5th max rank Neg/Pos:
        ip4=0
	ind4=0
        ip8=0
	ind8=0
c=      IF((phiup_pre.gt.philo_pre).and.(abs(phiwup).gt.abs(phiwlo))
c=	+.or.(phiwup.eq.0.)) THEN
       IF (iopt.eq.1) THEN
C  max => ind4
	do ip=1,4
	if (imaxpn(ip).gt.ind4) then
	ind4=imaxpn(ip)
	ip4=ip  
	endif
   	enddo
C
        irtot(4)=ip4
	imtot(4)=ind4
	imaxpn(ip4)=0
	iwout(4)=iwc(ip4)
C    max => ind8
	do ip=5,8
	if (imaxpn(ip).gt.ind8) then
	ind8=imaxpn(ip)
	ip8=ip  
	endif
   	enddo
C
        irtot(8)=ip8
	imtot(8)=ind8
	imaxpn(ip8)=0
	iwout(8)=iwc(ip8)
	                                             ELSE
C  max => ind4
C iopt=2
	do ip=5,8
	if (imaxpn(ip).gt.ind4) then
	ind4=imaxpn(ip)
	ip4=ip  
	endif
   	enddo
C
        irtot(4)=ip4
	imtot(4)=ind4
	imaxpn(ip4)=0
	iwout(4)=iwc(ip4)
C    max => ind8
	do ip=1,4
	if (imaxpn(ip).gt.ind8) then
	ind8=imaxpn(ip)
	ip8=ip  
	endif
   	enddo
C
        irtot(8)=ip8
	imtot(8)=ind8
	imaxpn(ip8)=0
	iwout(8)=iwc(ip8)
                                             ENDIF

C
  334	OPEN(11,file=outfile)
C Lati cycle
C
C
      DO 122 lat=1,71
C Longi cycle
	   DO 123 lon=0,72
C	 
      DO 125 iopt=1,8
      ncycle=irtot(iopt)
	iwxes=iwout(iopt)
      select case (ncycle)
	case (1)
   1	if((imapwu(lon,lat).eq.0).and.(icp11(lon,lat).eq.1)) then 
	 imapwu(lon,lat)=1                 !W 1st rank at W=+1
	 goto 125
	endif
	case (2)
   2	 if((imapwu(lon,lat).eq.0).and.(icp22(lon,lat).eq.2)) then 
	 imapwu(lon,lat)=2                 !W 1st rank at W=+1
	 goto 125
	endif
	case (3)
   3	 if((imapwu(lon,lat).eq.0).and.(icp33(lon,lat).eq.3)) then 
	 imapwu(lon,lat)=3                 !W 1st rank at W=+1
	 goto 125
	endif
	case (4)
   4	 if((imapwu(lon,lat).eq.0).and.(icp44(lon,lat).eq.4)) then 
	 imapwu(lon,lat)=4                 !W 1st rank at W=+1
	 goto 125
	endif
	case (5)
   5	 if((imapwu(lon,lat).eq.0).and.(icp01(lon,lat).eq.-1)) then 
	 imapwu(lon,lat)=-1                 !W 1st rank at W=+1
	 goto 125
	endif
	case (6)
   6	 if((imapwu(lon,lat).eq.0).and.(icp02(lon,lat).eq.-2)) then 
	 imapwu(lon,lat)=-2                 !W 1st rank at W=+1
	 goto 125
	endif
	case (7)
   7	if((imapwu(lon,lat).eq.0).and.(icp03(lon,lat).eq.-3)) then 
	 imapwu(lon,lat)=-3                 !W = +1 
	 goto 125
	endif
	case (8)
   8	if((imapwu(lon,lat).eq.0).and.(icp04(lon,lat).eq.-4)) then 
	 imapwu(lon,lat)=-4                 !W = +1 
	 goto 125
	endif
	END SELECT
  125	CONTINUE     ! end ranks cycle

  123 CONTINUE	 ! end longs
C
	write(*,32) (imapwu(k,lat),k=0,72)
	write(11,32) (imapwu(k,lat),k=0,72)
  122 CONTINUE     ! end lats
C
      ind01=0		 
  	ind02=0
 	ind03=0
	ind04=0
	ind11=0
	ind22=0
	ind33=0
	ind44=0
      iwx01=0		 
  	iwx02=0
 	iwx03=0
	iwx04=0
	iwx11=0
	iwx22=0
	iwx33=0
	iwx44=0
	do ilat=1,71
	do ilon=0,72
	if (icp01(ilon,ilat).eq.-1) ind01=ind01+1
	if (icp02(ilon,ilat).eq.-2) ind02=ind02+1
	if (icp03(ilon,ilat).eq.-3) ind03=ind03+1
	if (icp04(ilon,ilat).eq.-4) ind04=ind04+1
	if (icp11(ilon,ilat).eq.1) ind11=ind11+1
	if (icp22(ilon,ilat).eq.2) ind22=ind22+1
	if (icp33(ilon,ilat).eq.3) ind33=ind33+1
	if (icp44(ilon,ilat).eq.4) ind44=ind44+1
	if (imapwu(ilon,ilat).eq.-1) iwx01=iwx01+1
	if (imapwu(ilon,ilat).eq.-2) iwx02=iwx02+1
	if (imapwu(ilon,ilat).eq.-3) iwx03=iwx03+1
	if (imapwu(ilon,ilat).eq.-4) iwx04=iwx04+1
	if (imapwu(ilon,ilat).eq.1) iwx11=iwx11+1
	if (imapwu(ilon,ilat).eq.2) iwx22=iwx22+1
	if (imapwu(ilon,ilat).eq.3) iwx33=iwx33+1
	if (imapwu(ilon,ilat).eq.4) iwx44=iwx44+1
	enddo
	enddo
	write(*,38) iw01,nqua_1,medp_1,max_1,ind01,iwx01
	write(*,38) iw02,nqua_2,medp_2,max_2,ind02,iwx02
	write(*,38) iw03,nqua_3,medp_3,max_3,ind03,iwx03
	write(*,38) iw04,nqua_4,medp_4,max_4,ind04,iwx04
	write(*,38) iw11,nqua11,medp11,maxp11,ind11,iwx11
	write(*,38) iw22,nqua22,medp22,maxp22,ind22,iwx22
	write(*,38) iw33,nqua33,medp33,maxp33,ind33,iwx33
	write(*,38) iw44,nqua44,medp44,maxp44,ind44,iwx44
	write(12,38) iw01,nqua_1,medp_1,max_1,ind01,iwx01
	write(12,38) iw02,nqua_2,medp_2,max_2,ind02,iwx02
	write(12,38) iw03,nqua_3,medp_3,max_3,ind03,iwx03
	write(12,38) iw04,nqua_4,medp_4,max_4,ind04,iwx04
	write(12,38) iw11,nqua11,medp11,maxp11,ind11,iwx11
	write(12,38) iw22,nqua22,medp22,maxp22,ind22,iwx22
	write(12,38) iw33,nqua33,medp33,maxp33,ind33,iwx33
	write(12,38) iw44,nqua44,medp44,maxp44,ind44,iwx44
C
C
       close(unit=11)
C++	if (phiwup.ne.0.) phiup_pre=abs(phiwup)
C++	philo_pre=abs(phiwlo)
	if (phiwup.ne.0.) phiup_pre=phiwup
	philo_pre=phiwlo

C
  121 CONTINUE  ! End GIM-W maps for selected event
   49	close(unit=13)
	close(unit=12)
C
C 551	CONTINUE
C
	pause ' '
	STOP
	END
