program compo_ec implicit none integer::fcheck,jlag,ilag,icount1,icount,idmax,id,nmax,j1,j2,fnum,flength,ifront,jfront,i1,i2,mflag,nn,ii,jj,hjet,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,hmn2,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::maxdvdx,maxududx,maxvdvdx,maxvorjet,maxdudx,jetlon,jetlat,lat1,lon1,jlat,jlon,pi,es,ssp,slpup,slplow,lon2,lat2,minslp,hoge integer,dimension(imx,jmx,hmx)::fflag,fflag2,fpoint,fflag4,fflag3,fcold,fcold2,ftype real,dimension(1000)::ijet,jjet real,dimension(1000)::maxuwind,maxvwind,maxmwind,latmax,lonmax real,dimension(1000)::latfr,lonfr real,dimension(nntmx)::trid3,time3 real,dimension(imx,jmx,hmx)::fflag5,u,v,tmp,hgt,dtdx,dtdy,fmag,rh,zeta,vvel,qx,qy,ug,vg,pv,vor real,dimension(jmx)::co,lat real,dimension(imx)::lon real,dimension(imx,jmx,hmx2)::spfh real,dimension(imx,jmx,hhmx)::reh,rehu,rehv,ru,rv 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)::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(imx,jmx)::uz,vz,mwind,mwind1,dsigmdx,dsigmdy,dsigmdn,kthre,dudn,duds,vdvdx,ududx,vorjet,cape,sreh,ehi,mdvdz,mlapse,psurf,slp,tmp2m,dep,u10m,v10m,rh2m,sp2m,pt2m,cu,cv,suz,svz,shear_u,shear_v,mshear real,dimension(hmx)::dp 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::eclon,eclat,maxdudn,maxduds,dlon,delta,ber,mdlon,mdlat,bermax,r,dlam,dphi,shcr,dist,distmax,maxwind2,lonjet,latjet,lonec,latec real:: mdlons,mdlats,trids,times,minslps character::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.) r=6.3712e6 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='DJF' region='JPN' CAT='PO' ob=1 nob=0 nmax=105 undef=99999 do h=hmn,hmx if(h.eq.hmn) then dp(h)=(p(h+1)-p(h))/2. elseif(h.eq.hmx) then dp(h)=(p(h)-p(h-1))/2. else dp(h)=(p(h+1)-p(h-1))/2. endif enddo do j=jmn,jmx lat(j)=80.-1.25*(j-1) co(j)=cos(lat(j)/180.*pi) enddo do i=imn,imx if(region=='JPN') then lon(i)=80.+1.25*(i-1) else lon(i)=220.+1.25*(i-1) endif enddo open(57,file='jetdetect_DRATE05_100_'//CAT//'_'//SEA//'_long.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='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_190E',status='old',form='formatted') open(29,file='./EC_MINSLP/'//year//''//SEA//''//CAT//'_MINSLP_LONLAT'//region//'_190E',status='old',form='formatted') open(38,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_vavelowpass5day_Sjetconf_190E_1000_700',status='unknown',form='formatted') open(48,file='./EC_MINSLP/'//year//''//SEA//''//CAT//'_MINSLP_LONLAT'//region//'_vavelowpass5day_Sjetconf_190E_1000_700',status='unknown',form='formatted') open(39,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_vavelowpass5day_Njetconf_190E_1000_700',status='unknown',form='formatted') open(40,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_vavelowpass5day_Sjetdiff_190E_1000_700',status='unknown',form='formatted') open(41,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_vavelowpass5day_Njetdiff_190E_1000_700',status='unknown',form='formatted') open(42,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_relativeto_jet_250_vavelowpass5day_Sjetconf_190E',status='unknown',form='formatted') open(43,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_relativeto_jet_250_vavelowpass5day_Njetconf_190E',status='unknown',form='formatted') open(44,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_relativeto_jet_250_vavelowpass5day_Sjetdiff_190E',status='unknown',form='formatted') open(45,file='./EC_MAXDEV/'//year//''//SEA//''//CAT//'_MAXDEV_LONLAT'//region//'_relativeto_jet_250_vavelowpass5day_Njetdiff_190E',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 120 continue read(29,*,end=121) mdlons,mdlats,trids,times,minslps !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 ! do itime=time2,time2 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='old',form='formatted') do nnt=1,nntmx read(26,*,end=130) lon1,lat1,time3(nnt),trid3(nnt),ber ! 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') !open(51,file='/home/tochimoto/breeze1/tochimoto/runmean5day_jra55/runmean_'//region//'.'//year//'.'//mmonth//'_5day.data',status='old',form='unformatted') open(51,file='/home/tochimoto/breeze1/tochimoto/lowpass5day_jra55/lowpass_'//region//'.'//year//'.'//mmonth//'_5day.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) !((u(i,j,h),i=imn,imx),j=jmn,jmx) do h=hmn,hmx read(51) ((u(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h=hmn,hmx read(51) ((v(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h=hmn,hmx read(51) ((tmp(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h=hmn,hmx2 read(51) ((spfh(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h=hmn,hmx read(51) ((hgt(i,j,h),i=imn,imx),j=jmn,jmx) enddo fflag=0 fflag2=0 fcold=0 fcold2=0 hjet=12 hmn2=7 !vertical average of horizontal wind uz=0 vz=0 do h=hmn2,hjet do j=jmn,jmx do i=imn,imx uz(i,j)=uz(i,j)+u(i,j,h)*dp(h) vz(i,j)=vz(i,j)+v(i,j,h)*dp(h) enddo enddo enddo uz(:,:)=uz(:,:)/(p(hjet)-p(hmn2)) vz(:,:)=vz(:,:)/(p(hjet)-p(hmn2)) mwind=sqrt(uz*uz+vz*vz) !dudn=100. 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 dlam=r*1.25*2.*pi/180. dphi=r*1.25*2.*pi/180.*cos(jlat/180.*pi) vorjet(i,j)=(vz(i+1,j)-vz(i-1,j))/dlam-(uz(i,j-1)*cos((jlat+1.25)/180.*pi)-uz(i,j+1)*cos((jlat-1.25)/180.*pi))/dphi dphi=r*1.25*2.*pi/180. dlam=r*1.25*2.*pi/180.*cos(jlat/180.*pi) if(i.eq.imn) then ududx(i,j)=(mwind(i+1,j)-mwind(i,j))/dlam*2. elseif(i.eq.imx) then ududx(i,j)=(mwind(i,j)-mwind(i-1,j))/dlam*2. else ududx(i,j)=(mwind(i+1,j)-mwind(i-1,j))/dlam endif if(j.eq.jmn) then vdvdx(i,j)=(mwind(i,j)-mwind(i,j+1))/dphi*2. elseif(j.eq.jmx) then vdvdx(i,j)=(mwind(i,j-1)-mwind(i,j))/dphi*2. else vdvdx(i,j)=(mwind(i,j-1)-mwind(i,j+1))/dphi endif dudn(i,j)=vz(i,j)/mwind(i,j)*ududx(i,j)-uz(i,j)/mwind(i,j)*vdvdx(i,j) duds(i,j)=uz(i,j)/mwind(i,j)*ududx(i,j)+vz(i,j)/mwind(i,j)*vdvdx(i,j) 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 fflag4=0 fflag5=0 it=it+1 !write(6,*) 'it2= ',it if(it.eq.itmax/2+1) then out=out+1 endif if(it.eq.itmax/2+1) then !condition for shear gradient kthre=dudn*dudn+mwind*dsigmdn shcr=-0.3e-09 !find dudn=0. h = hjet id=0 write(6,*) 'id=', id write(57) ((fflag5(i2,j2,hjet),i2=imn,imx),j2=jmn,jmx) write(57) ((mwind(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((dudn(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((duds(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((uz(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((vz(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((kthre(i2,j2),i2=imn,imx),j2=jmn,jmx) write(57) ((dsigmdn(i2,j2),i2=imn,imx),j2=jmn,jmx) write(6,*) it,itime,time3(nnt), iyear,'out' write(6,*) lat2,lon2,ice1,jce1 !write(6,*) it 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(i+ice1,j+jce1) endif enddo enddo enddo do j=-yst,yen sco3(j,it)=sco3(j,it)+1 !write(6,*) it,sco3(j,it),co(j+jce1) enddo mwind1=mwind mflag=1 maxvorjet=vorjet(ice1,jce1) maxdudn=dudn(ice1,jce1) maxduds=duds(ice1,jce1) !maxduds=maxuwind(idmax)/maxmwind(idmax)*ududx(ice1,jce1)+maxvwind(idmax)/maxmwind(idmax)*vdvdx(ice1,jce1) write(6,*) 'maxdudn, maxduds', maxdudn,maxduds !maxdvdx=maxuwind(idmax)*real(ijet(idmax))+maxvwind(idmax)*real(-jjet(idmax)) !maxdudn=maxuwind(idmax)*real(ijet(idmax))+maxvwind(idmax)*real(-jjet(idmax)) !maxduds=maxuwind(idmax)*real(ijet(idmax))+maxvwind(idmax)*real(-jjet(idmax)) if(maxdudn.lt.0.and.maxduds.gt.0) then write(38,*) mdlon,mdlat,trid2,time2,bermax write(48,*) minslps write(42,*) jetlon,jetlat,ijet(1),jjet(1) write(6,*) 'sjetconf' write(6,*) mdlon,mdlat,trid2,time2,bermax endif if(maxdudn.gt.0.and.maxduds.gt.0) then write(39,*) mdlon,mdlat,trid2,time2,bermax write(43,*) jetlon,jetlat,ijet(1),jjet(1) write(6,*) 'njetconf' write(6,*) mdlon,mdlat,trid2,time2,bermax endif if(maxdudn.lt.0.and.maxduds.lt.0) then write(40,*) mdlon,mdlat,trid2,time2,bermax write(44,*) jetlon,jetlat,ijet(1),jjet(1) write(6,*) 'sjetdiff' write(6,*) mdlon,mdlat,trid2,time2,bermax endif if(maxdudn.gt.0.and.maxduds.lt.0) then write(41,*) mdlon,mdlat,trid2,time2,bermax write(45,*) jetlon,jetlat,ijet(1),jjet(1) write(6,*) 'njetdiff' write(6,*) mdlon,mdlat,trid2,time2,bermax endif endif endif ! pick up EC enddo !time loop endif close(50) close(51) enddo !nntloop 130 continue close(26) !write(6,*) n,nt,nnt,iyear enddo ! itimeloop endif !write(6,*) n,nt enddo !ntloop 121 continue !write(6,*) n,nt close(28) enddo !year loop ! write(6,*) 'out=',out,'nout=',nout endprogram compo_ec