program permutation implicit none integer::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,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::kst, ken integer::hmx,hmn,hhmn,hhmx,hmx2 integer::k,ememax,emem,emem1 parameter(imx=51,jmx=51,hmx=37,hmn=1,jmn=1,imn=1,tmn=1,hhmn=1,hhmx=40,hmx2=27,nntmx=10000,itmax=5000,nmx=100000,ememax=5000) integer,dimension(itmax)::nn2 integer, dimension(imx, jmx) :: nu, nv, nvor, ntmp, nspfh, nhgt, nwind !integer,dimension(imx,jmx,hmx)::nv,nvvel,ntmp,npv,nhgt,nvor !integer,dimension(imx,jmx,hmx2)::nsp real::lat1,lon1,jlat,jlon,pi,es,ssp,slpup,slplow,lon2,lat2,minslp,hoge,varvn,varsrehn,s12,varcapen,varpvn,varspn,varhgtn,varvveln,vartmpn,rmin,rmin2 real,dimension(itmax)::r,rr real,dimension(nntmx)::trid3,time3 real,dimension(jmx)::co,lat real,dimension(imx,jmx,hmx)::spfh, u, v, vor, tmp, hgt, mwind real,dimension(imx,jmx,hmx)::coslant, codef, coept, cofmag real,dimension(imx,jmx,hmx)::def,slant,ept,fmag real, dimension(imx, jmx, itmax) :: cosreh, cosreh2, cosreh3 real, dimension(imx, jmx, itmax) :: cou, cou2, cou3 real, dimension(imx, jmx, itmax) :: cov, cov2, cov3 real, dimension(imx, jmx, itmax) :: covor, covor2, covor3 real, dimension(imx, jmx, itmax) :: cotmp, cotmp2, cotmp3 real, dimension(imx, jmx, itmax) :: cospfh, cospfh2, cospfh3 real, dimension(imx, jmx, itmax) :: cohgt, cohgt2, cohgt3 real, dimension(imx, jmx, itmax) :: comwind, comwind2, comwind3 real, dimension(imx, jmx, hmx) :: cou4, cov4, covor4, cotmp4, cospfh4, cohgt4, comwind4 real,dimension(jmx,itmax)::sco3,sco4,sco5 real,dimension(imx,jmx)::avesreh1,avesreh2,dsreh,dcosreh,psreh,avecold1,avecold2,dcold,dcocold,pcold,avewarm1,avewarm2,dwarm,dcowarm,pwarm real, dimension(imx, jmx) :: aveu1, aveu2, du, dcou, pu real, dimension(imx, jmx) :: avev1, avev2, dv, dcov, pv real, dimension(imx, jmx) :: avevor1, avevor2, dvor, dcovor, pvor real, dimension(imx, jmx) :: avetmp1, avetmp2, dtmp, dcotmp, ptmp real, dimension(imx, jmx) :: avespfh1, avespfh2, dspfh, dcospfh, pspfh real, dimension(imx, jmx) :: avehgt1, avehgt2, dhgt, dcohgt, phgt real, dimension(imx, jmx) :: avemwind1, avemwind2, dmwind, dcomwind, pwind real, dimension(imx, jmx, ememax) :: u2d, v2d, vor2d, tmp2d, spfh2d, hgt2d, mwind2d real, dimension(imx, jmx, ememax) :: su1,su2, sv1, sv2, svor1, svor2, stmp1, stmp2 real, dimension(imx, jmx, ememax) :: svvel1, svvel2, sspfh1, sspfh2, spv1, spv2 ,shgt1,shgt2 real, dimension(imx, jmx, ememax) :: smwind1,smwind2 real::snum,snum1,snum2,snum3,tnum,tnum2,ss,snum11,snum12,snum13,snum21,snum22,snum23 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 character :: jetpos1*8, jetpos2*8 logical::fast !itime = 13 !CAT='PO' !CAT2='OJ' SEA='SON' region='JPN' su1 = 0 sv1 = 0 svor1 = 0 stmp1 = 0 sspfh1 = 0 shgt1 = 0 smwind1 = 0 su2 = 0 sv2 = 0 svor2 = 0 stmp2 = 0 sspfh2 = 0 smwind2 = 0 snum11 = 0 snum12 = 0 snum13 = 0 snum21 = 0 snum22 = 0 snum23 = 0 cou = 0 cou2 = 0 cov = 0 cov2 = 0 covor = 0 covor2 = 0 cotmp = 0 cotmp2 = 0 cospfh = 0 cospfh2 = 0 cohgt = 0 cohgt2 = 0 comwind = 0 comwind2 = 0 jetpos1 = 'Njetconf' jetpos2 = 'Njetdiff' !write(6,*) 'read namelist' !namelist/namdim/CAT,CAT2,itime !write(6,*) 'category' !write(6,*) CAT,' ',CAT2,' ',itime !read(5,namdim) CAT='PO' itime = 9 kst = 3 ken = 3 it = itime if(it.le.9) then write(kt1,'(i1)') it kt='0'//kt1 else write(kt,'(i2)') it endif open(69,file='pvalue_for_permutation_ecjet_env_rotate_850hPa_'//CAT//'_'//kt//'_'//jetpos1//'_'//jetpos2//'.data',status='unknown',form='unformatted') emem=1 do k=kst,ken ss = 1 if(k.eq.1) then SEA = 'MAM' elseif(k.eq.2) then SEA = 'SON' elseif(k.eq.3) then SEA = 'DJF' endif write(6, *) SEA open(68,file='./'//CAT//'/emem_'//kt//'_env_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos1//'_160E_190E.data',status='unknown',form='unformatted') !open(68,file=''//CAT//'/emem_exspfront600_ktvar_DRATE05_100_'//CAT//'_'//SEA//'_'//kt//'.data',status='unknown',form='unformatted') kk=68 do it=1,ememax do h = hmn, hmx read(kk,end=150) ((u(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((v(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((tmp(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx2 read(kk) ((spfh(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((hgt(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((vor(i,j,h),i=imn,imx),j=jmn,jmx) enddo h = 7 do j = jmn, jmx do i = imn, imx u2d(i,j,emem) = u(i,j,h) v2d(i,j,emem) = v(i,j,h) tmp2d(i,j,emem) = tmp(i,j,h) spfh2d(i,j,emem) = spfh(i,j,h) hgt2d(i,j,emem) = hgt(i,j,h) vor2d(i,j,emem) = vor(i,j,h) mwind2d(i,j,emem) = sqrt(u(i,j,h)*u(i,j,h) + v(i,j,h)*v(i,j,h)) enddo enddo !read(kk,end=150) ((sreh(i,j,emem),i=imn,imx),j=jmn,jmx) !read(kk) ((cold(i,j,emem),i=imn,imx),j=jmn,jmx) !read(kk) ((warm(i,j,emem),i=imn,imx),j=jmn,jmx) emem = emem + 1 ss = ss + 1 enddo 150 continue if(k.eq.1) then snum11 = ss endif if(k.eq.2) then snum12 = ss endif if(k.eq.3) then snum13 = ss endif enddo !kloop emem1 = emem write(6,*) 'emem1', emem1 do k=kst,ken ss = 1 if(k.eq.1) then SEA = 'MAM' elseif(k.eq.2) then SEA = 'SON' elseif(k.eq.3) then SEA = 'DJF' endif !open(68,file=''//CAT//'/emem_exnspfront600_ktvar_DRATE05_100_'//CAT//'_'//SEA//'_'//kt//'.data',status='unknown',form='unformatted') open(68,file='./'//CAT//'/emem_'//kt//'_env_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos2//'_160E_190E.data',status='unknown',form='unformatted') kk=68 do it=1,ememax do h = hmn, hmx read(kk,end=151) ((u(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((v(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((tmp(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx2 read(kk) ((spfh(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((hgt(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((vor(i,j,h),i=imn,imx),j=jmn,jmx) enddo h = 7 do j = jmn, jmx do i = imn, imx u2d(i,j,emem) = u(i,j,h) v2d(i,j,emem) = v(i,j,h) tmp2d(i,j,emem) = tmp(i,j,h) spfh2d(i,j,emem) = spfh(i,j,h) hgt2d(i,j,emem) = hgt(i,j,h) vor2d(i,j,emem) = vor(i,j,h) mwind2d(i,j,emem) = sqrt(u(i,j,h)*u(i,j,h) + v(i,j,h)*v(i,j,h)) enddo enddo !read(kk,end=151) ((sreh(i,j,emem),i=imn,imx),j=jmn,jmx) !read(kk) ((cold(i,j,emem),i=imn,imx),j=jmn,jmx) !read(kk) ((warm(i,j,emem),i=imn,imx),j=jmn,jmx) emem = emem + 1 ss = ss + 1 enddo 151 continue if(k.eq.1) then snum21 = ss endif if(k.eq.2) then snum22 = ss endif if(k.eq.3) then snum23 = ss endif enddo !kloop write(6,*) 'emem', emem do s=kst,ken if(s.eq.1) then SEA = 'MAM' snum = snum11 elseif(s.eq.2) then SEA = 'SON' snum = snum12 else SEA = 'DJF' snum = snum13 endif open(57,file='compo_env_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos1//'_160E_190E.data',status='unknown',form='unformatted') tnum=snum11+snum12+snum13 write(6,*) snum11,snum12,snum13 do it=1,17 !AO do h=hmn,hmx read(57) ((cou4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((cov4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((covor4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((cotmp4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx2 read(57) ((cospfh4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((cohgt4(i,j,h),i=imn,imx),j=imn,imx) enddo h = 7 do j = jmn, jmx do i = imn, imx cou3(i,j,it) = cou4(i,j,h) cov3(i,j,it) = cov4(i,j,h) covor3(i,j,it) = covor4(i,j,h) cotmp3(i,j,it) = cotmp4(i,j,h) cospfh3(i,j,it) = cospfh4(i,j,h) cohgt3(i,j,it) = cohgt4(i,j,h) comwind3(i,j,it) = sqrt(cou4(i,j,h)*cou4(i,j,h) + cov4(i,j,h)*cov4(i,j,h)) enddo enddo !read(57) ((cosreh3(i,j,it),i=imn,imx),j=jmn,jmx) !read(57) ((cocold3(i,j,it),i=imn,imx),j=jmn,jmx) !read(57) ((cowarm3(i,j,it),i=imn,imx),j=jmn,jmx) do j=jmn,jmx do i=imn,imx cou(i,j,it)=cou(i,j,it)+cou3(i,j,it)*snum/tnum cov(i,j,it)=cov(i,j,it)+cov3(i,j,it)*snum/tnum covor(i,j,it)=covor(i,j,it)+covor3(i,j,it)*snum/tnum cotmp(i,j,it)=cotmp(i,j,it)+cotmp3(i,j,it)*snum/tnum cospfh(i,j,it)=cospfh(i,j,it)+cospfh3(i,j,it)*snum/tnum cohgt(i,j,it)=cohgt(i,j,it)+cohgt3(i,j,it)*snum/tnum comwind(i,j,it)=comwind(i,j,it)+comwind3(i,j,it)*snum/tnum enddo enddo enddo if(s.eq.1) then SEA='MAM' snum=snum21 elseif(s.eq.2) then SEA='SON' snum=snum22 else SEA='DJF' snum=snum23 endif open(58,file='compo_env_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos2//'_160E_190E.data',status='unknown',form='unformatted') !write(6, *) CAT, SEA tnum2=snum21+snum22+snum23 write(6,*) snum21,snum22,snum23 do it=1, 17 !PO do h=hmn,hmx read(58) ((cou4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((cov4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((covor4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((cotmp4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx2 read(58) ((cospfh4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((cohgt4(i,j,h),i=imn,imx),j=imn,imx) enddo h = 7 do j = jmn, jmx do i = imn, imx cou3(i,j,it) = cou4(i,j,h) cov3(i,j,it) = cov4(i,j,h) covor3(i,j,it) = covor4(i,j,h) cotmp3(i,j,it) = cotmp4(i,j,h) cospfh3(i,j,it) = cospfh4(i,j,h) cohgt3(i,j,it) = cohgt4(i,j,h) comwind3(i,j,it) = sqrt(cou4(i,j,h)*cou4(i,j,h) + cov4(i,j,h)*cov4(i,j,h)) enddo enddo do j=jmn,jmx do i=imn,imx cou2(i,j,it)=cou2(i,j,it)+cou3(i,j,it)*snum/tnum2 cov2(i,j,it)=cov2(i,j,it)+cov3(i,j,it)*snum/tnum2 covor2(i,j,it)=covor2(i,j,it)+covor3(i,j,it)*snum/tnum2 cotmp2(i,j,it)=cotmp2(i,j,it)+cotmp3(i,j,it)*snum/tnum2 cospfh2(i,j,it)=cospfh2(i,j,it)+cospfh3(i,j,it)*snum/tnum2 cohgt2(i,j,it)=cohgt2(i,j,it)+cohgt3(i,j,it)*snum/tnum2 comwind2(i,j,it)=comwind2(i,j,it)+comwind3(i,j,it)*snum/tnum2 enddo enddo enddo enddo !season loop it= itime do j=jmn,jmx do i=imn,imx dcou(i,j) = cou(i,j,it) - cou2(i,j,it) dcov(i,j) = cov(i,j,it) - cov2(i,j,it) dcovor(i,j) = covor(i,j,it) - covor2(i,j,it) dcotmp(i,j) = cotmp(i,j,it) - cotmp2(i,j,it) dcospfh(i,j) = cospfh(i,j,it) - cospfh2(i,j,it) dcohgt(i,j) = cohgt(i,j,it) - cohgt2(i,j,it) dcomwind(i,j) = comwind(i,j,it) - comwind2(i,j,it) enddo enddo nu = 0 nv = 0 nvor = 0 ntmp = 0 nspfh = 0 nhgt = 0 nwind = 0 do nnt=1,nmx write(6, *) 'repeat ', nnt aveu1 = 0. avev1 = 0. avevor1 = 0. avetmp1 = 0. avespfh1 = 0. avehgt1 = 0. avemwind1 = 0. aveu2 = 0. avev2 = 0. avevor2 = 0. avetmp2 = 0. avespfh2 = 0. avehgt2 = 0. avemwind2 = 0. rmin2=0. do n=1,emem call random_seed call random_number(r(1:emem)) rmin=1. do nn=1,emem if(r(nn).lt.rmin.and.r(nn).gt.rmin2) then rmin=r(nn) nnew=nn endif enddo rr(n)=rmin nn2(n)=nnew rmin2=rmin enddo !do m=1,itmax !write(6,*) rr(m),nn2(m) !enddo do it=1,emem do j=jmn,jmx do i=imn,imx if(it.le.emem1) then su1(i,j,it) = u2d(i,j,nn2(it)) sv1(i,j,it) = v2d(i,j,nn2(it)) svor1(i,j,it) = vor2d(i,j,nn2(it)) stmp1(i,j,it) = tmp2d(i,j,nn2(it)) sspfh1(i,j,it) = spfh2d(i,j,nn2(it)) shgt1(i,j,it) = hgt2d(i,j,nn2(it)) smwind1(i,j,it) = mwind2d(i,j,nn2(it)) else su2(i,j,it - emem1) = u2d(i,j,nn2(it)) sv2(i,j,it - emem1) = v2d(i,j,nn2(it)) svor2(i,j,it - emem1) = vor2d(i,j,nn2(it)) stmp2(i,j,it - emem1) = tmp2d(i,j,nn2(it)) sspfh2(i,j,it - emem1) = spfh2d(i,j,nn2(it)) shgt2(i,j,it - emem1) = hgt2d(i,j,nn2(it)) smwind2(i,j,it - emem1) = mwind2d(i,j,nn2(it)) endif enddo enddo enddo !do it=1,itmax !write(6,*) sreh(20,20,nn2(it)),ssreh2(20,20,it),it !enddo do it = 1, emem1 do j = jmn, jmx do i = imn, imx aveu1(i,j) = aveu1(i,j) + su1(i,j,it)/real(emem1) avev1(i,j) = avev1(i,j) + sv1(i,j,it)/real(emem1) avevor1(i,j) = avevor1(i,j) + svor1(i,j,it)/real(emem1) avetmp1(i,j) = avetmp1(i,j) + stmp1(i,j,it)/real(emem1) avespfh1(i,j) = avespfh1(i,j) + sspfh1(i,j,it)/real(emem1) avehgt1(i,j) = avehgt1(i,j) + shgt1(i,j,it)/real(emem1) avemwind1(i,j) = avemwind1(i,j) + smwind1(i,j,it)/real(emem1) enddo enddo enddo do it=1,emem-emem1+1 do j=jmn,jmx do i=imn,imx aveu2(i,j) = aveu2(i,j)+su2(i,j,it)/real(emem-emem1+1) avev2(i,j) = avev2(i,j)+sv2(i,j,it)/real(emem-emem1+1) avevor2(i,j) = avevor2(i,j)+svor2(i,j,it)/real(emem-emem1+1) avetmp2(i,j) = avetmp2(i,j)+stmp2(i,j,it)/real(emem-emem1+1) avespfh2(i,j) = avespfh2(i,j)+sspfh2(i,j,it)/real(emem-emem1+1) avehgt2(i,j) = avehgt2(i,j)+shgt2(i,j,it)/real(emem-emem1+1) avemwind2(i,j) = avemwind2(i,j)+smwind2(i,j,it)/real(emem-emem1+1) enddo enddo enddo do j=jmn,jmx do i=imn,imx du(i,j) = aveu1(i,j) - aveu2(i,j) dv(i,j) = avev1(i,j) - avev2(i,j) dvor(i,j) = avevor1(i,j) - avevor2(i,j) dtmp(i,j) = avetmp1(i,j) - avetmp2(i,j) dspfh(i,j) = avespfh1(i,j) - avespfh2(i,j) dhgt(i,j) = avehgt1(i,j) - avehgt2(i,j) dmwind(i,j) = avemwind1(i,j) - avemwind2(i,j) if(abs(du(i,j)).ge.abs(dcou(i,j))) then nu(i,j) = nu(i,j) + 1 endif if(abs(dv(i,j)).ge.abs(dcov(i,j))) then nv(i,j) = nv(i,j) + 1 endif if(abs(dvor(i,j)).ge.abs(dcovor(i,j))) then nvor(i,j) = nvor(i,j) + 1 endif if(abs(dtmp(i,j)).ge.abs(dcotmp(i,j))) then ntmp(i,j) = ntmp(i,j) + 1 endif if(abs(dspfh(i,j)).ge.abs(dcospfh(i,j))) then nspfh(i,j) = nspfh(i,j) + 1 endif if(abs(dhgt(i,j)).ge.abs(dcohgt(i,j))) then nhgt(i,j) = nhgt(i,j) + 1 endif if(abs(dmwind(i,j)) .ge. abs(dcomwind(i,j))) then nwind(i,j) = nwind(i,j) + 1 endif enddo enddo !write(6,*) avesreh1(20,20),avesreh2(20,20),dsreh(20,20),dcosreh(20,20) !write(6,*) nnt enddo !nntloop do j=jmn,jmx do i=imn,imx pu(i,j) = (1.+real(nu(i,j)))/(real(nmx+1)) pv(i,j) = (1.+real(nv(i,j)))/(real(nmx+1)) pvor(i,j) = (1.+real(nvor(i,j)))/(real(nmx+1)) ptmp(i,j) = (1.+real(ntmp(i,j)))/(real(nmx+1)) pspfh(i,j) = (1.+real(nspfh(i,j)))/(real(nmx+1)) phgt(i,j) = (1.+real(nhgt(i,j)))/(real(nmx+1)) pwind(i,j) = (1.+real(nwind(i,j)))/(real(nmx+1)) enddo enddo !psreh=psreh*2. !pcold=pcold*2. !pwarm=pwarm*2. write(6,*) pvor(20,20),phgt(10,10),ptmp(20,20) write(6,*) du(20,20),dhgt(20,20),avetmp1(10,20), avetmp2(10, 20), dspfh(20, 20) write(6,*) avehgt1(10,20), avehgt2(10, 20), avespfh1(20, 20) write(69) ((pu(i,j),i=imn,imx),j=jmn,jmx) write(69) ((pv(i,j),i=imn,imx),j=jmn,jmx) write(69) ((pvor(i,j),i=imn,imx),j=jmn,jmx) write(69) ((ptmp(i,j),i=imn,imx),j=jmn,jmx) write(69) ((pspfh(i,j),i=imn,imx),j=jmn,jmx) write(69) ((phgt(i,j),i=imn,imx),j=jmn,jmx) write(69) ((pwind(i,j),i=imn,imx),j=jmn,jmx) endprogram permutation