program est_kdensity
implicit none
integer::k,nnt,nntmx,iyear,syear,eyear,n,nmx,nn,nmax,nnew,nnew2,nnew3,i,j,imx,jmx
parameter(nntmx=10000,syear=1979,eyear=2016,nmx=100000,imx=121,jmx=65)
real,dimension(nmx)::lon,lat,lon2,lat2,long,lati
real,dimension(imx)::rlon,rlat
real,dimension(imx,jmx)::f,kfunc1,kfunc2
real::lon1,lat1,time3,trid3,ber,lonmin,lonmin2,latmin,latmin2,lonq4,lonq34,latq4,latq34,avelon,avelat,savelon,savelat,h1,h2,pi
character::SEA*3, CAT*2, region*3, year*4, jetpos*8
CAT='PO'
region='JPN'
jetpos='Sjetdiff'


do k=3,3
if(k.eq.1) then
SEA='MAM'
elseif(k.eq.2) then
SEA='SON'
else
SEA='DJF'
endif

lonmin2=0
latmin2=0
avelon=0.
avelat=0.
pi=acos(-1.)
nn=1

open(11,file='kernel_maxdev_'//CAT//'_'//SEA//'_'//jetpos//'_190E_1000_700.dat',status='unknown',form='unformatted')

 do iyear=syear,eyear
  write(year,'(i4)') iyear

open(28,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_vavelowpass5day_'//jetpos//'_190E_1000_700',status='old',form='formatted')  
do nnt=1,nntmx
  read(28,*,end=130) lon1,lat1,time3,trid3,ber
  if(lon1.ne.999.and.lat1.ne.999.and.lat1.gt.0) then
  long(nn)=lon1
  lati(nn)=lat1
  write(6,*) long(nn),lati(nn)
  nn=nn+1
  endif
  enddo
130 continue
 enddo ! yearloop

nmax=nn-1

do i=1,imx
  rlon(i)=80.+1.25*(i-1)
  rlat(i)=0.+1.25*(i-1)
enddo

nnew=0
lon=long
lat=lati

do n=1,nmax
 lonmin=1000.
 do nn=n,nmax

  if(lon(nn).eq.lonmin2.and.nn.ne.nnew2) then
   lonmin=lon(nn)
   !lon(nn)=lon(nn)*100
   nnew=nn
   !write(6,*) nn,lonmin
  elseif(lon(nn).lt.lonmin.and.lon(nn).gt.lonmin2) then
   lonmin=lon(nn)
   nnew=nn
  endif

  !if(lon(n).eq.lon(nn).and.n.ne.nn) then
  ! write(6,*) lon(n),lon(nn),n,nn
  !endif

 enddo
lon(nnew)=lon(n)
lon(n)=lonmin
lon2(n)=lonmin
lonmin2=lonmin
nnew2=nnew

enddo
nnew=0
do n=1,nmax
 latmin=1000.
 do nn=n,nmax

  if(lat(nn).eq.latmin2.and.nn.ne.nnew2) then
   latmin=lat(nn)
   !lat(nn)=lat(nn)*100
   nnew=nn
  elseif(lat(nn).lt.latmin.and.lat(nn).gt.latmin2) then
   latmin=lat(nn)
   nnew=nn
  endif
 enddo
lat(nnew)=lat(n)
lat(n)=latmin
lat2(n)=latmin
latmin2=latmin
nnew2=nnew
!write(6,*) lon2(n),lat2(n),lat(n),n
enddo

do n=1,nmax
!write(6,*) lon2(n),lat2(n),lat(n),n
avelon=avelon+lon2(n)/real(nmax)
avelat=avelat+lat2(n)/real(nmax)
enddo

do n=1,nmax
 savelon=savelon+(lon2(n)-avelon)**2/real(nmax-1)
 savelat=savelat+(lat2(n)-avelat)**2/real(nmax-1)
enddo

savelon=sqrt(savelon)
savelat=sqrt(savelat)

lonq34=lon2(nmax-nmax/4)
lonq4=lon2(nmax/4)

latq34=lat2(nmax-nmax/4)
latq4=lat2(nmax/4)

write(6,*) savelon,savelat
write(6,*) lonq34-lonq4,latq34-latq4
write(6,*) min(0.9*savelon,0.75*(lonq34-lonq4)),0.9*savelon,0.75*(lonq34-lonq4)





h1= min(0.9*savelon,0.75*(lonq34-lonq4))/(nmax**0.2)
h2= min(0.9*savelat,0.75*(latq34-latq4))/(nmax**0.2)
write(6,*) h1,h2
!h1=0.5
!h2=0.5
f=0.
do i=1,imx
 do j=1,jmx
do nn=1,nmax
   kfunc1(i,j)=1./sqrt(2*pi)*exp(-((rlon(i)-long(nn))/h1)**2./2.)
   kfunc2(i,j)=1./sqrt(2*pi)*exp(-((rlat(j)-lati(nn))/h2)**2./2.)
   f(i,j)=f(i,j)+1./(real(nmax)*h1*h2)*kfunc1(i,j)*kfunc2(i,j)!*((rlon(i)-long(nn))/h1*(rlat(j)-lati(nn))/h2)
enddo
 enddo
enddo

!f = f * nmax

!write(6,*) f
write(6,*) 'number of EC = ', nmax 
write(11) ((f(i,j),i=1,imx),j=1,jmx)
write(11) ((kfunc1(i,j),i=1,imx),j=1,jmx)
write(11) ((kfunc2(i,j),i=1,imx),j=1,jmx)

enddo !kloop

endprogram est_kdensity