program compo_ec
implicit none
integer::ii,jj,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
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)
real::lat1,lon1,jlat,jlon,pi,es,ssp,slpup,slplow,lon2,lat2,minslp,hoge
real,dimension(nntmx)::trid3,time3
real,dimension(imx,jmx,hmx)::u,v,tmp,hgt,dtdx,dtdy,fmag,rh,zeta,vvel,qx,qy,ug,vg,pv,vor
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(-20:20,-20:20,hmx,itmax)::cou,cov,covor,cotmp,cou2,cov2,covor2,cotmp2,cou3,cov3,covor3,cotmp3,covvel,covvel2,covvel3,cohgt,cohgt2,cohgt3,coqx,coqy,coqx2,coqy2,coqx3,coqy3,coug,covg,coug2,covg2,coug3,covg3,copv,copv2,copv3
real,dimension(-20:20,-20:20,hhmx,itmax)::coreh,coreh2,coreh3,corehu,corehu2,corehu3,corehv,corehv2,corehv3,coru,corv,coru2,corv2,coru3,corv3,couz,covz,couz2,covz2,couz3,covz3,vmag,vmag2,vmag3,angle,angle2,angle3
real,dimension(-20:20,-20:20,hmx2,itmax)::cosp,cosp2,cosp3
real,dimension(-20:20,-20:20,itmax)::cocold,cowarm,cosreh,cosreh2,cosreh3,cocape,cocape2,cocape3,cosuz,cosvz,cosuz2,cosvz2,cosuz3,cosvz3,coshearu,coshearv,&
                              comshear,coshearu2,coshearv2,comshear2,coshearu3,coshearv3,comshear3,cocu,cocv,cocu2,cocv2,cocu3,cocv3
real,dimension(-20:20,itmax)::sco3,sco4,sco5
real,dimension(-20:20,-20:20)::newi
real,dimension(-20:20,-20:20)::newj
real,dimension(-20:20,-20:20)::sreh_rot,cold_rot,warm_rot
real,dimension(imx,jmx)::cold,warm,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::ang,dlon,dlat,ai,aj,ber,mdlon,mdlat,bermax,hoge1,hoge2,hoge3
character::kt1*1,kt*2,mmonth*2,year*4,SEA*3,region*3,CAT*2
logical::fast
fast=.false. ! a little output variables
nmx=20000
ntmx=20000
pi=acos(-1.)
write(6,*) pi
xst=20
yst=20
xen=20
yen=20
xst2=2
yst2=5
xen2=5
yen2=5
xst3=5
yst3=5
xen3=2
yen3=5
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='SON'
region='USA'
CAT='AO'
cosreh=0
cocold=0
cowarm=0

ob=1
nob=0

undef=99999


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


open(57,file='compo_coldfront_DRATE05_100_'//CAT//'_'//SEA//'_rotate.data',status='unknown',form='unformatted')
open(60,file='each_coldfront_DRATE05_100_'//CAT//'_'//SEA//'_rotate.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 iyear=syear,eyear
write(year,'(i4)') iyear


open(28,file='./'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'',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(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='./'//year//''//SEA//''//CAT//'_track_EC_'//region//'',status='unknown',form='formatted')
open(27,file='./'//year//''//SEA//''//CAT//'_movement_EC_'//region//'',status='unknown',form='formatted')
  do nnt=1,nntmx

    read(26,*,end=130) lon1,lat1,time3(nnt),trid3(nnt),ber
    read(27,*,end=150) 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(50,file='/data2/tochimoto/work/jra55/data/'//region//'/NEW_TFPWF1000_0_JRA55_'//region//'.'//year//'.'//mmonth//'.data',status='old',form='unformatted')
open(50,file='/home/tochimoto/breeze1/tochimoto/jra55/FGENE/'//region//'/NEW_TFPCF1000_0_JRA55_'//region//'.'//year//'.'//mmonth//'.data',status='old',form='unformatted')


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(50) !((u(i,j,h),i=imn,imx),j=jmn,jmx)
  enddo

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

   read(50) ((sreh(i,j),i=imn,imx),j=jmn,jmx)
   read(50) !((u(i,j,h),i=imn,imx),j=jmn,jmx)
   read(50) ((cold(i,j),i=imn,imx),j=jmn,jmx)
   read(50) ((warm(i,j),i=imn,imx),j=jmn,jmx)

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

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

ang=atan2(dlat,dlon)
!if(dlat.lt.0) then
! ang=ang+2*pi
!endif

write(6,*) 'cos(ang)',cos(ang),'sin(ang)',sin(ang)
write(6,*) 'ang',ang

 !write(6,*) it

 do i=-xst,xen
 do j=-yst,yen   
   newi(i,j)=real(i)*cos(ang)+real(j)*sin(ang)
   newj(i,j)=real(j)*cos(ang)-real(i)*sin(ang)
 enddo
 enddo

!do i=-xst,xen
! write(6,*) 'newi',newi(i),i
!enddo

!do j=-yst,yen
! write(6,*) 'newj',newj(j),j
!enddo

!interporate
do i=-xst,xen
 do j=-yst,yen
  do ii=-xst,xen
   do jj=-yst,yen
     if(newi(i,j).ge.real(ii).and.newi(i,j).le.real(ii+1).and.newj(i,j).ge.real(jj).and.newj(i,j).le.real(jj+1)) then
      ai = newi(i,j) - real(ii)
      aj = newj(i,j) - real(jj)
      sreh_rot(i,j)=ai*aj*sreh(ii+ice1+1,jj+jce1+1)+(1.-ai)*aj*sreh(ii+ice1,jj+jce1+1)+ai*(1-aj)*sreh(ii+ice1+1,jj+jce1)+(1.-ai)*(1.-aj)*sreh(ii+ice1,jj+jce1) 
      cold_rot(i,j)=ai*aj*cold(ii+ice1+1,jj+jce1+1)+(1.-ai)*aj*cold(ii+ice1,jj+jce1+1)+ai*(1-aj)*cold(ii+ice1+1,jj+jce1)+(1.-ai)*(1.-aj)*cold(ii+ice1,jj+jce1) 
      warm_rot(i,j)=ai*aj*warm(ii+ice1+1,jj+jce1+1)+(1.-ai)*aj*warm(ii+ice1,jj+jce1+1)+ai*(1-aj)*warm(ii+ice1+1,jj+jce1)+(1.-ai)*(1.-aj)*warm(ii+ice1,jj+jce1) 
     endif
   enddo
  enddo
 enddo
enddo


if(it.eq.itmax/2+1) then

 !write(60) ((sreh_rot(i,j),i=-xst,xen),j=-yst,yen)
 !write(60) ((sreh(i,j),i=-xst+ice1,xen+ice1),j=-yst+jce1,yen+jce1)
 !write(60) ((cold_rot(i,j),i=-xst,xen),j=-yst,yen)
 !write(60) ((warm_rot(i,j),i=-xst,xen),j=-yst,yen)

endif



 do h=hmn,hmx
 do i=-xst,xen
 do j=-yst,yen

  if(h.le.hmx2) then
  endif

  if(h.eq.hmn) then
   cosreh(i,j,it)=cosreh(i,j,it)+sreh_rot(i,j)
   cocold(i,j,it)=cocold(i,j,it)+cold_rot(i,j)
   cowarm(i,j,it)=cowarm(i,j,it)+warm_rot(i,j)
  endif

 enddo
 enddo
 enddo


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

open(58,file='variance_front_emem_rotate_DRATE05_100_'//CAT//'_'//kt//'.data',status='unknown',position='append',form='unformatted')


 write(58) ((sreh_rot(i,j),i=-xst,xen),j=-yst,yen)
 write(58) ((cold_rot(i,j),i=-xst,xen),j=-yst,yen)
 write(58) ((warm_rot(i,j),i=-xst,xen),j=-yst,yen)

 do j=-yst,yen
  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(50)

  enddo !nntloop
130 continue
150 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=-xst,xen
 do j=-yst,yen

if(fast) then
else

   if(h.eq.hmn) then
   cosreh(i,j,it)=cosreh(i,j,it)/sco3(j,it)
   cocold(i,j,it)=cocold(i,j,it)/sco3(j,it)
   cowarm(i,j,it)=cowarm(i,j,it)/sco3(j,it)
   endif

endif !fast
   enddo
  enddo
 enddo


enddo

do it=1,itmax

 write(57) ((cosreh(i,j,it),i=-xst,xen),j=-yst,yen)
 write(57) ((cocold(i,j,it),i=-xst,xen),j=-yst,yen)
 write(57) ((cowarm(i,j,it),i=-xst,xen),j=-yst,yen)

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