program falsedetection
implicit none
integer::ii,jj,s,i,j,h,imn,imx,jmn,jmx,syear,eyear,iyear,smonth,emonth,imonth,ice1,jce1,xst,xen,yst,yen,month1,n,nto,nto2,nmx,nt,nnt,out,nout,nto3,it,itmax,tt,kk,tmx1,tmx2,m,nn,num
integer::xst2,xen2,yst2,yen2,xst3,xen3,yst3,yen3,time,time2,time4,itime,t,tmx,tmn,ntmx,nntmx,trid,trid2,xst1,xen1,yst1,yen1,iit,nnew
integer::day1,hour1,hmx,hmn,hhmn,hhmx,hmx2,ob,nob
integer::k,ememax,emem,emem1
parameter(imx=51,jmx=51,hmx=37,syear=1995,eyear=2012,smonth=3,emonth=5,hmn=1,jmn=1,imn=1,tmn=1,hhmn=1,hhmx=40,hmx2=27,nntmx=10000,tmx1=2172,tmx2=1909,itmax=5000,nmx=10000,ememax=2000)
integer,dimension(itmax)::nn2
real,dimension(itmax)::r,rr
real,dimension(nntmx)::trid3,time3
real,dimension(jmx)::co,lat
real,dimension(jmx,itmax)::sco3,sco4,sco5
real,dimension(imx,jmx)::psreh, pbaroc, pdiab
real,dimension(imx,jmx)::pu, pv, ptmp, pspfh, ppv, phgt
real,dimension(imx,jmx,hmx2,itmax)::ssp1,ssp2
real::snum,snum1,snum2,snum3,tnum,tnum2
real,    dimension(1000) :: pfdr, minsreh3, minbaroc3, mindiab3
integer, dimension(1000) :: isreh3, jsreh3
integer, dimension(1000) :: ibaroc3, jbaroc3
integer, dimension(1000) :: idiab3, jdiab3
real::p(hmx)=(/1000., 975., 950., 925., 900., 875., 850., 825., 800., 775.,750., 700., &
650., 600., 550., 500., 450., 400., 350., 300., 250., 225., 200., 175., 150.,125., 100., &
70., 50., 30., 20., 10.,7., 5., 3., 2., 1./)
character::kt1*1,mmonth*2,year*4,SEA*3,region*3,kt*2,CAT*2,CAT2*2,jetpos*8
character:: jetpos1*8, jetpos2*8
logical::fast
integer::imn2,jmn2,imx2,jmx2
integer::nmax,isreh,jsreh,ndata

real :: alpha, minpsreh, minpsreh2, lonsreh, latsreh
real :: minpbaroc, minpbaroc2, lonbaroc, latbaroc
real :: minpdiab, minpdiab2, londiab, latdiab


CAT='PO'
CAT2='PO'
jetpos1='Njetconf'
jetpos2='Njetdiff'
SEA='DJF'
region='JPN'

it=5

do it = 9, 9, 2
if(it.le.9) then
write(kt1,'(i1)') it
kt='0'//kt1
else
write(kt,'(i2)') it
endif

!open(69,file='pvalue_for_permutation_frontxy_KT'//kt//'_'//CAT//'_'//CAT2//'_'//jetpos//'.data',status='unknown',form='unformatted')

open(69,file='pvalue_for_permutation_ecjet_vvel_pv_rotate_700hPa_'//CAT//'_'//kt//'_'//jetpos1//'_'//jetpos2//'.data',status='unknown',form='unformatted')

open(70,file='false_detection_for_permutation_ano_VVEL_rotate_700hPa_KT'//kt//'_'//CAT//'_'//jetpos1//'_'//jetpos2//'.txt',status='unknown',form='formatted')
open(71,file='false_detection_for_permutation_ano_UG_rotate_700hPa_KT'//kt//'_'//CAT//'_'//jetpos1//'_'//jetpos2//'.txt',status='unknown',form='formatted')
open(72,file='false_detection_for_permutation_ano_VG_rotate_700hPa_KT'//kt//'_'//CAT//'_'//jetpos1//'_'//jetpos2//'.txt',status='unknown',form='formatted')
open(73,file='false_detection_for_permutation_ano_qy_rotate_700hPa_KT'//kt//'_'//CAT//'_'//jetpos1//'_'//jetpos2//'.txt',status='unknown',form='formatted')
open(74,file='false_detection_for_permutation_ano_EPV_rotate_700hPa_KT'//kt//'_'//CAT//'_'//jetpos1//'_'//jetpos2//'.txt',status='unknown',form='formatted')

! read data

 read(69) ((ppv(i, j),  i = imn, imx), j = jmx, jmn, -1)
 read(69) ((pv(i, j),  i = imn, imx), j = jmx, jmn, -1)
 read(69) ((pu(i, j), i = imn, imx), j = jmx, jmn, -1)
 read(69) !((ptmp(i,j), i = imn, imx),j = jmx, jmn, -1)
 read(69) ((pspfh(i,j), i = imn, imx),j = jmx, jmn, -1) 
 read(69) ((phgt(i, j), i = imn, imx), j = jmx, jmn, -1)


!write(6,*) psreh


! sort pvalues

!imn2=10
!imx2=42
!jmn2=10
!jmx2=42
!imn2=9
!imx2=43
!jmn2=9
!jmx2=43
imn2=18
imx2=34
jmn2=18
jmx2=34
!imn2=14
!imx2=38
!jmn2=14
!jmx2=38

call fdr(imx, jmx, imn2, imx2, jmn2, jmx2, ppv, ibaroc3, jbaroc3, nmax, pfdr, minbaroc3)

write(6, *) nmax

do n = 1, nmax
  lonbaroc = 230.+ (ibaroc3(n)-1) * 1.5
  latbaroc = 10.+ (jbaroc3(n)-1) * 1.5
  !write(70,*)  minsreh3(n),isreh3(n),jsreh3(n),pfdr(n)
  write(70,*)  lonbaroc, latbaroc, pfdr(n)
  write(6,*)  minbaroc3(n), lonbaroc, latbaroc, pfdr(n)
enddo

nmax = 0
pfdr = 0.
call fdr(imx, jmx, imn2, imx2, jmn2, jmx2, pu, idiab3, jdiab3, nmax, pfdr, mindiab3)

write(6, *) nmax

do n = 1, nmax
  londiab = 230. + (idiab3(n)-1) * 1.5
  latdiab = 10. + (jdiab3(n)-1) * 1.5
  !write(70,*)  minsreh3(n),isreh3(n),jsreh3(n),pfdr(n)
  write(71,*)  londiab, latdiab, pfdr(n)
  !write(6,*)  mindiab3(n), londiab, latdiab, pfdr(n)
enddo

nmax = 0
pfdr = 0.
call fdr(imx, jmx, imn2, imx2, jmn2, jmx2, pv, idiab3, jdiab3, nmax, pfdr, mindiab3)

write(6, *) nmax

do n = 1, nmax
  londiab = 230. + (idiab3(n)-1) * 1.5
  latdiab = 10. + (jdiab3(n)-1) * 1.5
  !write(70,*)  minsreh3(n),isreh3(n),jsreh3(n),pfdr(n)
  write(72,*)  londiab, latdiab, pfdr(n)
  !write(6,*)  mindiab3(n), londiab, latdiab, pfdr(n)
enddo



nmax = 0
pfdr = 0.
call fdr(imx, jmx, imn2, imx2, jmn2, jmx2, pspfh, idiab3, jdiab3, nmax, pfdr, mindiab3)

write(6, *) nmax

do n = 1, nmax
  londiab = 230. + (idiab3(n)-1) * 1.5
  latdiab = 10. + (jdiab3(n)-1) * 1.5
  !write(70,*)  minsreh3(n),isreh3(n),jsreh3(n),pfdr(n)
  write(73,*)  londiab, latdiab, pfdr(n)
  !write(6,*)  mindiab3(n), londiab, latdiab, pfdr(n)
enddo


nmax = 0
pfdr = 0.
call fdr(imx, jmx, imn2, imx2, jmn2, jmx2, phgt, idiab3, jdiab3, nmax, pfdr, mindiab3)

write(6, *) nmax

do n = 1, nmax
  londiab = 230. + (idiab3(n)-1) * 1.5
  latdiab = 10. + (jdiab3(n)-1) * 1.5
  !write(70,*)  minsreh3(n),isreh3(n),jsreh3(n),pfdr(n)
  write(74,*)  londiab, latdiab, pfdr(n)
  !write(6,*)  mindiab3(n), londiab, latdiab, pfdr(n)
enddo


write(6,*) it
!write(6,*) (jmx2-jmn2+1)*(imx2-imn2+1)


enddo

endprogram falsedetection


subroutine fdr(imx, jmx, imn2, imx2, jmn2, jmx2, var, ivar3, jvar3, nmax, pfdr, minvar3)
integer                                  :: i, j, ii, jj
integer                                  :: nmax
integer, intent(in)                      :: imx, jmx
integer, intent(in)                      :: imn2, imx2, jmn2, jmx2
real, intent(in), dimension(imx, jmx)    :: var
integer                                  :: n
integer                                  :: num, ndata
real                                     :: alpha
real                                     :: minvar, minvar2
integer                                  :: ivar, jvar
integer, intent(out), dimension(1000)       :: ivar3, jvar3 
real, intent(out), dimension(1000)                    :: minvar3, pfdr 

ndata=(jmx2 - jmn2 + 1)*(imx2 - imn2 + 1)

!write(6, *) var

! sort pvalues
alpha = 0.1
num=0
minvar2 = -1000
do jj = jmn2, jmx2
 do ii = imn2, imx2
  minvar = 1.
  do j = jmn2, jmx2
   do i = imn2, imx2
    if(var(i, j) .gt. minvar2 .and. var(i, j) .lt. minvar) then
     minvar = var(i,j)
     ivar = i
     jvar = j
    endif
   enddo
  enddo

   num = num+1
   minvar2 = minvar
   minvar3(num) = minvar
   ivar3(num) = ivar
   jvar3(num) = jvar
   pfdr(num) = num*alpha/real(ndata)
!   write(6,*) minvar3(num),ivar3(num),jvar3(num),pfdr(num)
 enddo
enddo

nmax = 0
do n = 1, num
 if( minvar3(n) .lt. pfdr(n) .and. n.gt.nmax ) then
  nmax = n
 endif
enddo

write(6, *) num

end subroutine fdr