program compo_ec
implicit none
integer::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
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
integer:: day1,hour1,hmx,hmn,hhmn,hhmx,hmx2,ob,nob,nfs,undef
integer:: xst4, xen4, yst4, yen4
parameter(imx=105,jmx=65,hmx=37,syear=1979,eyear=2016,smonth=3,emonth=5,hmn=1,jmn=1,imn=1,tmn=1,hhmn=1,hhmx=40,hmx2=27,itmax=17,nntmx=20000)
parameter(xst=20, xen=20, yst=20, yen=20, xst2=25, xen2=25, yst2=25, yen2=25, xst3=35, xen3=35, yst3=35, yen3=35)
parameter(xst4=40, xen4=40, yst4=40, yen4=40)
real::lat1,lon1,jlat,jlon,pi,es,ssp,slpup,slplow,lon2,lat2,minslp,hoge
real,dimension(nntmx)::trid3,time3
real,dimension(imx,jmx,hmx)::ept,def,slant,fmag,tmp,hgt,dtdx,dtdy,rh,zeta,vvel,qx,qy,ug,vg,pv,vor,def1,def2,slant1,slant2,fmag1,fmag2,ept1,ept2
real,dimension(jmx)::co,lat
real,dimension(imx,jmx,hmx2)::spfh
real,dimension(imx,jmx,hhmx)::reh,rehu,rehv,ru,rv,uz,vz
real,dimension(-xst4:xen4,-yst4:yen4,hmx,itmax)::codef,coslant,cofmag,coept
real,dimension(-yst4:yen4,itmax)::sco3,sco4,sco5
real,dimension(imx,jmx)::cape,sreh,ehi,mdvdz,mlapse,psurf,slp,tmp2m,dep,u10m,v10m,rh2m,sp2m,pt2m,cu,cv,suz,svz,shear_u,shear_v,mshear
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./)
real::ber,mdlon,mdlat,bermax
character::jetpos*8,nmem1*1,nmem*2,mmonth*2,year*4,SEA*3,region*3,CAT*2
logical::fast

integer::ii, jj, in, jn, irmax, jrmax
real,dimension(-yst4:yen4)::newy,y_ec,ypos
real,dimension(-xst4:xen4, -yst4:yen4)      :: newi, newj
real,dimension(-xst4:xen4,-yst4:yen4,hmx)::ffmag,front,fcold,fwarm
real,dimension(-xst4:xen4,-yst4:yen4,hmx)::ffmag_rot,front_rot,fcold_rot,fwarm_rot
real,dimension(-xst4:xen4)::newx
real,dimension(-yst4:yen4,-yst4:yen4)::x_ec,xpos
real::r,dx,dy,defx,defy,ai,aj,dphi
real::dist,dist2,defx2,defy2
real :: dlat, dlon, ang, wang, magwind
real :: hoge3
integer :: hoge1, hoge2

fast=.false. ! a little output variables


nmx=20000
ntmx=20000
pi=acos(-1.)
write(6,*) pi
xst1=10
xen1=10
yst1=10
yen1=10
!slpup=100500.
!slplow=99500.
slpup=110000.
slplow=50000.
ice1=0
jce1=0
out=0
nout=0
sco3=0.
sco4=0.
!SEA='MAM'
!region='JPN'
!CAT='OJ'

ob=1
nob=0



write(6,*) 'read namelist'
namelist/namdim/CAT,SEA,jetpos,region
write(6,*) 'category'
write(6,*) CAT,' ',SEA,' ',jetpos,' ',region
read(5,namdim)


r=6378.e3
dphi=1.25*pi/180.

dx=1.25e5
dy=1.25e5

undef=99999


do i=-40,40
 newx(i)=dx*real(i)
enddo

do j=-40,40
 newy(j)=dy*real(j)
enddo

do j=jmn,jmx
 lat(j)=80.-1.25*(j-1)
 co(j)=cos(lat(j)/180.*pi)
enddo

open(57,file='compo_ptfgene_diaba_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos//'_160E_190E.data',status='unknown',form='unformatted')

!open(57,file='compo_ec_'//CAT//'_'//SEA//'.data',status='unknown',form='unformatted')
!open(58,file='compo_ec_nout_'//SEA//'_50year_oriv1abz0f0_min992.data',status='unknown',form='unformatted')

do it = 1,itmax

 if(it.le.9) then
  write(nmem1,'(i1)') it
  nmem = '0'//nmem1
 else
  write(nmem,'(i2)') it
 endif

  open(80+it,file='./'//CAT//'/emem_'//nmem//'_ptfgene_diaba_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos//'_160E_190E.data',status='unknown',form='unformatted')


enddo

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='unknown',form='formatted')

do nt=1,ntmx
  !write(6,*) n
  !read(22,*,end=120) nto2,trid2,time2
  read(28,*,end=120) mdlon,mdlat,trid2,time2,bermax
 write(6,*) 'KT',iyear,trid2,time2,bermax
!if(trid.eq.trid2.and.minslp.ge.slplow.and.minslp.le.slpup) then 
minslp=100000.
if(mdlon .ge. 160. .and. minslp.ge.slplow.and.minslp.le.slpup.and.bermax.ge.0.5.and.bermax.le.10.0.and.trid2.ne.undef) then 

it=0
iit=-4
 do itime=time2-itmax/2,time2+itmax/2
write(year,'(i4)') iyear
write(6,*) 'itime is from',time2-itmax/2,'to',time2+itmax/2
write(6,*) 'keytime=',time2
!write(6,*) 'iyear=', iyear,'itime=',itime
if(SEA.eq.'MAM') then

if(itime.le.31*4) then
month1=3
mmonth='03'
tmx=31*4
if(mod(itime,4)==0) then
day1=itime/4
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1
hour1=12
endif

elseif(itime.le.31*4+30*4) then
month1=4
mmonth='04'
tmx=30*4
if(mod(itime,4)==0) then
day1=itime/4-31
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-31
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-31
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-31
hour1=12
endif



write(6,*) 'pass',iyear,trid,minslp,time2


else
month1=5
mmonth='05'
tmx=31*4
if(mod(itime,4)==0) then
day1=itime/4-31-30
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-31-30
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-31-30
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-31-30
hour1=12
endif

endif

endif
!write(6,*) month1


if(SEA.eq.'SON') then

if(itime.le.30*4) then
month1=9
mmonth='09'
tmx=30*4
if(mod(itime,4)==0) then
day1=itime/4
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1
hour1=12
endif

elseif(itime.le.31*4+30*4) then
month1=10
mmonth='10'
tmx=31*4
if(mod(itime,4)==0) then
day1=itime/4-30
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-30
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-30
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-30
hour1=12
endif



write(6,*) 'pass',iyear,trid,minslp,time2


else
month1=11
mmonth='11'
tmx=30*4
if(mod(itime,4)==0) then
day1=itime/4-31-30
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-31-30
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-31-30
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-31-30
hour1=12
endif

endif

endif
!write(6,*) month1

if(SEA.eq.'DJF') then

if(itime.le.31*4) then
month1=0
mmonth='12'
tmx=31*4
if(mod(itime,4)==0) then
day1=itime/4
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1
hour1=12
endif

elseif(itime.le.31*4+31*4) then
month1=1
mmonth='01'
tmx=31*4
if(mod(itime,4)==0) then
day1=itime/4-31
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-31
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-31
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-31
hour1=12
endif

else
month1=2
mmonth='02'

if(mod(iyear,4).eq.0) then
tmx=29*4
else
tmx=28*4
endif

if(mod(itime,4)==0) then
day1=itime/4-31-31
hour1=18
elseif(mod(itime,4)==1) then
day1=itime/4+1-31-31
hour1=0
elseif(mod(itime,4)==2) then
day1=itime/4+1-31-31
hour1=6
elseif(mod(itime,4)==3) then
day1=itime/4+1-31-31
hour1=12
endif
endif

endif

!write(6,*) 'month=', mmonth,'  ',year
iit=iit+1
!write(6,*) month1


open(26,file='./EC_TRACK/'//year//''//SEA//''//CAT//'_track_EC_'//region//'_190E',status='unknown',form='formatted')
open(27,file='./EC_TRACK/'//year//''//SEA//''//CAT//'_movement_EC_'//region//'_190E',status='unknown',form='formatted')

  do nnt=1,nntmx

    read(26,*,end=130) lon1,lat1,time3(nnt),trid3(nnt),ber
    read(27,*) dlon,dlat,hoge1,hoge2,hoge3
    
    !    write(6,*) 'read track',lon1,lat1,time3(nnt),trid3(nnt),ber
!if(trid2.ne.trid3(nnt).and.time3(nnt).eq.itime.and.month1.ne.3.) then
!it=it+1
!endif
!write(6,*) 'trid2=',trid2,'trid3=',trid3(nnt)
if(trid2.eq.trid3(nnt)) then

!write(6,*) 'pass2',trid2

if(itmax.ge.2) then
do tt=1,itmax/2
!do tt=1,8
if(time3(nnt).eq.time2-itmax/2+tt.and.it.eq.0) then
 it=it+tt
endif
enddo
endif

!write(6,*) 'it1= ', it
!if(time3(nnt).eq.time2-itmax+3.and.it.eq.0) then
! it=it+2
!endif

!if(time3(nnt).eq.time2-itmax+4.and.it.eq.0) then
! it=it+3
!endif

!if(time3(nnt).eq.time2-itmax+5.and.it.eq.0) then
! it=it+4
!endif
  
!if(trid2.eq.trid3.and.nto.ge.1) then
!write(6,*) 'pass2'


if(month1.eq.0) then
write(year,'(i4)') iyear-1
else
write(year,'(i4)') iyear
endif


open(53,file='/home/e_tochi/data2/breeze1/tochimoto/jra55/FGENE/'//region//'/NEW_FGENE_DIABA_PT1000_0_JRA55_'//region//'.'//year//'.'//mmonth//'.data', status='old', form='unformatted', access = 'sequential', action = 'READ')


do t=tmn,tmx

if(month1.eq.0) then
time4=t
endif

if(month1.eq.1) then
time4=t+31*4
endif

if(month1.eq.2) then
time4=t+31*4+31*4
endif

if(month1.eq.3) then
time4=t
endif

if(month1.eq.4) then
time4=t+31*4
endif

if(month1.eq.5) then
time4=t+30*4+31*4
endif

if(month1.eq.9) then
time4=t
endif

if(month1.eq.10) then
time4=t+30*4
endif

if(month1.eq.11) then
time4=t+30*4+31*4
endif


do h=hmn,hmx
 read(53) ((def1(i,j,h),i=imn,imx),j=jmn,jmx)
enddo

do h=hmn,hmx
 read(53) ((slant1(i,j,h),i=imn,imx),j=jmn,jmx)
enddo

do h=hmn,hmx
 read(53) ((ept1(i,j,h),i=imn,imx),j=jmn,jmx)
enddo

do h=hmn,hmx
 read(53) ((fmag1(i,j,h),i=imn,imx),j=jmn,jmx)
enddo

do h=hmn,hmx
 read(53) !((fmag1(i,j,h),i=imn,imx),j=jmn,jmx)
enddo


!def=(def1-def2)*1.e10
!slant=(slant1-slant2)*1.e10
!def=def1*1.e5*3600*6.
!slant=slant1*1.e5*3600.*6.
def=def1
slant=slant1
ept=ept1
fmag=fmag1


do i=imn,imx
do j=jmn,jmx

 jlat=80.-1.25*real(j-1)
 if(region=='USA') then
 jlon=220.+1.25*real(i-1)
 else
 jlon=80.+1.25*real(i-1)
 endif

 if(lat1.ge.jlat-1.25*0.5.and.lat1.le.jlat+1.25*0.5.and.lon1.ge.jlon-1.25*0.5.and.lon1.le.jlon+1.25*0.5) then
  ice1=i
  jce1=j
  lon2=jlon
  lat2=jlat
 endif

enddo
enddo
!write(6,*) 'itime==time4?', itime,time4,time3(nnt)
    if(itime.eq.time4.and.itime.eq.time3(nnt).and.minslp.ge.slplow.and.minslp.le.slpup) then
 
 it=it+1
!write(6,*) 'it2= ',it
if(it.eq.itmax/2+1) then
out=out+1
endif


!低気圧中心に相対的なxyを計算
do j=-yst2,yen2
 do i=-xst2,xen2
  x_ec(i,j) = r*co(j+jce1)*dphi*i
  y_ec(j) = r*dphi*j
 enddo
enddo


!xy座標に内挿
do h=1,hmx
 do j = -yst3,yen3
  do i = -xst3,xen3
   dist2=1.e30
   dist2=1.e30
   do jj = -yst2,yen2
    do ii = -xst2,xen2


    ! defx=(x_ec(ii,jj)-newx(i))*(x_ec(ii+1,jj)-newx(i))
    ! defy=(y_ec(jj)-newy(j))*(y_ec(jj+1)-newy(j))

    defx=min(abs(x_ec(ii,jj)-newx(i)),abs(x_ec(ii+1,jj)-newx(i)))
    defy=min(abs(y_ec(jj)-newy(j)),abs(y_ec(jj+1)-newy(j)))
    dist=sqrt(defx*defx+defy*defy)

     !if(x_ec(ii,jj).le.newx(i).and.x_ec(ii+1,jj).ge.newx(i).and.y_ec(jj).le.newy(j).and.y_ec(jj+1).ge.newy(j).and.dist.lt.dist2) then
     if(x_ec(ii,jj).le.newx(i).and.x_ec(ii+1,jj).ge.newx(i).and.y_ec(jj).le.newy(j).and.y_ec(jj+1).ge.newy(j)) then
      ai = abs((x_ec(ii,jj)-newx(i))/(x_ec(ii+1,jj)-x_ec(ii,jj)))
      aj = abs((y_ec(jj)-newy(j))/(y_ec(jj+1)-y_ec(jj)))
      write(6,*) ai,aj
      if(ii+ice1.ge.imn.and.jj+jce1.ge.jmn.and.ii+1+ice1.le.imx.and.jj+1+jce1.le.jmx) then
      
      front(i,j,h) = ai*aj*def(ii+1+ice1,jj+1+jce1,h)+(1.-ai)*aj*def(ii+ice1,jj+1+jce1,h)+ai*(1-aj)*def(ii+1+ice1,jj+jce1,h)+(1.-ai)*(1.-aj)*def(ii+ice1,jj+jce1,h)
      fcold(i,j,h) = ai*aj*slant(ii+1+ice1,jj+1+jce1,h)+(1.-ai)*aj*slant(ii+ice1,jj+1+jce1,h)+ai*(1-aj)*slant(ii+1+ice1,jj+jce1,h)+(1.-ai)*(1.-aj)*slant(ii+ice1,jj+jce1,h)
      fwarm(i,j,h) = ai*aj*ept(ii+1+ice1,jj+1+jce1,h)+(1.-ai)*aj*ept(ii+ice1,jj+1+jce1,h)+ai*(1-aj)*ept(ii+1+ice1,jj+jce1,h)+(1.-ai)*(1.-aj)*ept(ii+ice1,jj+jce1,h)
      ffmag(i,j,h) = ai*aj*fmag(ii+1+ice1,jj+1+jce1,h)+(1.-ai)*aj*fmag(ii+ice1,jj+1+jce1,h)+ai*(1-aj)*fmag(ii+1+ice1,jj+jce1,h)+(1.-ai)*(1.-aj)*fmag(ii+ice1,jj+jce1,h)

       defx2=min(abs(x_ec(ii,jj)-newx(i)),abs(x_ec(ii+1,jj)-newx(i)))
       defy2=min(abs(y_ec(jj)-newy(j)),abs(y_ec(jj+1)-newy(j)))
       dist2=sqrt(defx2*defx2+defy2*defy2)

      endif 
     endif

    enddo
   enddo
  enddo
 enddo
enddo


! rotation


ang=atan2(dlat,dlon)



 do i = -xst4,xen4
  do j = -yst4,yen4
   newi(i,j) = real(i) * cos(ang) + real(j) * sin(ang)
   newj(i,j) = real(j) * cos(ang) - real(i) * sin(ang)
  enddo
 enddo

irmax = 0
jrmax = 0


do h = hmn, hmx
 call interp2d(front(:, :, h), front_rot(:, :, h), newi, newj, xen2, yen2, irmax, jrmax, xen4, yen4)
enddo

do h = hmn, hmx
 call interp2d(fcold(:, :, h), fcold_rot(:, :, h), newi, newj, xen2, yen2, irmax, jrmax, xen4, yen4)
enddo

do h = hmn, hmx
 call interp2d(fwarm(:, :, h), fwarm_rot(:, :, h), newi, newj, xen2, yen2, irmax, jrmax, xen4, yen4)
enddo

do h = hmn, hmx
 call interp2d(ffmag(:, :, h), ffmag_rot(:, :, h), newi, newj, xen2, yen2, irmax, jrmax, xen4, yen4)
enddo




do h=hmn,hmx
 write(80+it) ((front_rot(i,j,h),i=-xst2,xen2),j=-yst2,yen2)
enddo

do h=hmn,hmx
 write(80+it) ((fcold_rot(i,j,h),i=-xst2,xen2),j=-yst2,yen2)
enddo

do h=hmn,hmx
 write(80+it) ((fwarm_rot(i,j,h),i=-xst2,xen2),j=-yst2,yen2)
enddo

do h=hmn,hmx
 write(80+it) ((ffmag_rot(i,j,h),i=-xst2,xen2),j=-yst2,yen2)
enddo


 write(6,*) it,itime,time3(nnt), iyear,'out'
 write(6,*) lat2,lon2,ice1,jce1

 !write(6,*) it
 do h=hmn,hmx
   do i=-xst2,xen2
    do j=-yst2,yen2
     codef(i,j,h,it)   = codef(i,j,h,it)   + front_rot(i,j,h)
     coslant(i,j,h,it) = coslant(i,j,h,it) + fcold_rot(i,j,h)
     coept(i,j,h,it)   = coept(i,j,h,it)   + fwarm_rot(i,j,h)
     cofmag(i,j,h,it)  = cofmag(i,j,h,it)  + ffmag_rot(i,j,h)
    enddo
   enddo
  enddo


 do j=-yst2,yen2
  sco3(j,it)=sco3(j,it)+1
  !write(6,*) it,sco3(j,it),co(j+jce1)
 enddo

    endif ! pick up EC

  

enddo !time loop

endif



close(53)

  enddo !nntloop
130 continue
close(26)
close(27)
!write(6,*) n,nt,nnt,iyear
 enddo ! itimeloop
endif

!write(6,*) n,nt
 enddo  !ntloop
  
120 continue
!write(6,*) n,nt
close(28)

enddo !year loop

 do it=1,itmax
  do h=hmn,hmx
   do i=-xst2,xen2
    do j=-yst2,yen2

    codef(i,j,h,it)=codef(i,j,h,it)/sco3(j,it)
    coslant(i,j,h,it)=coslant(i,j,h,it)/sco3(j,it)
    coept(i,j,h,it)=coept(i,j,h,it)/sco3(j,it)
    cofmag(i,j,h,it)=cofmag(i,j,h,it)/sco3(j,it)
 
   enddo
  enddo
 enddo


enddo

do it=1,itmax



do h=hmn,hmx
 write(57) ((codef(i,j,h,it),i=-xst2,xen2),j=-yst2,yen2) !ldiab
enddo

do h=hmn,hmx
 write(57) ((coslant(i,j,h,it),i=-xst2,xen2),j=-yst2,yen2) !cdiab
enddo


do h=hmn,hmx
 write(57) ((coept(i,j,h,it),i=-xst2,xen2),j=-yst2,yen2) !vdiab
enddo

do h=hmn,hmx
 write(57) ((cofmag(i,j,h,it),i=-xst2,xen2),j=-yst2,yen2) !sum of diaba
enddo

enddo
!
write(6,*) 'out=',out,'nout=',nout
endprogram compo_ec

subroutine interp2d(var, varrot, newi, newj, ilag, jlag, irmax, jrmax, ilag2, jlag2)

integer             :: ii, jj
integer, intent(in) :: ilag, jlag, ilag2, jlag2, irmax, jrmax
real, intent(in),  dimension(-ilag2 : ilag2, -jlag2 : jlag2)    :: newi, newj
real, intent(in),  dimension(-ilag2 : ilag2, -jlag2 : jlag2)    :: var
real, intent(out), dimension(-ilag2 : ilag2, -jlag2 : jlag2)    :: varrot
real              :: ai, aj

            do j = - jlag, jlag
              do i = - ilag, ilag
                ii = nint(newi(i, j))
                jj = nint(newj(i, j))
                if(ii .ge. -ilag + 1 .and. jj .ge. -jlag + 1 .and. ii .le. ilag - 1 .and. jj .le. jlag - 1) then
                  ai = newi(i, j) - real(ii)
                  aj = newj(i, j) - real(jj)
                   varrot(i, j) = ai * aj * var(ii + irmax + 1, jj + jrmax + 1) &
                           + (1. - ai)*aj*var(ii + irmax, jj + jrmax + 1) + &
                           ai * (1 - aj) * var(ii + irmax + 1, jj + jrmax) + &
                           (1. - ai) * (1. - aj)*var(ii + irmax, jj + jrmax)
                 endif
              enddo
             enddo


end subroutine interp2d