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 :: hlev 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 !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 real,dimension(imx,jmx,hmx) :: ldiab, cdiab, vdiab, tdiab real,dimension(imx,jmx,hmx) :: coslant, codef, coept, cofmag real,dimension(imx,jmx,hmx) :: def,slant,ept,fmag real, dimension(imx, jmx, itmax) :: coldiab, coldiab2, coldiab3 real, dimension(imx, jmx, itmax) :: cotdiab, cotdiab2, cotdiab3 real, dimension(imx, jmx, itmax) :: cocdiab, cocdiab2, cocdiab3 real, dimension(imx, jmx, itmax) :: covdiab, covdiab2, covdiab3 real, dimension(imx, jmx, hmx) :: coldiab4, cocdiab4, covdiab4, cotdiab4 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, ememax) :: cdiab2d, ldiab2d, vdiab2d, tdiab2d 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::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 su2 = 0 sv2 = 0 svor2 = 0 stmp2 = 0 sspfh2 = 0 shgt2 = 0 snum11=0 snum12=0 snum13=0 snum21=0 snum22=0 snum23=0 coldiab = 0 cocdiab2 = 0 cocdiab = 0 cocdiab2 = 0 covdiab = 0 covdiab2 = 0 cotdiab = 0 cotdiab2 = 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 hlev = 12 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_diabatic_rotate_'//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//'_diabatic_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) ((ldiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((cdiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((vdiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo tdiab = ldiab + cdiab + vdiab h = hlev do j = jmn, jmx do i = imn, imx ldiab2d(i,j,emem) = ldiab(i,j,h) cdiab2d(i,j,emem) = cdiab(i,j,h) vdiab2d(i,j,emem) = vdiab(i,j,h) tdiab2d(i,j,emem) = tdiab(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//'_diabatic_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos2//'_160E_190E.data',status='old',form='unformatted') kk=68 do it=1,ememax do h = hmn, hmx read(kk,end=151) ((ldiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((cdiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo do h = hmn, hmx read(kk) ((vdiab(i,j,h),i=imn,imx),j=jmn,jmx) enddo tdiab = ldiab + cdiab + vdiab h = hlev do j = jmn, jmx do i = imn, imx ldiab2d(i,j,emem) = ldiab(i,j,h) cdiab2d(i,j,emem) = cdiab(i,j,h) vdiab2d(i,j,emem) = vdiab(i,j,h) tdiab2d(i,j,emem) = tdiab(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_diabatic_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos1//'_160E_190E.data',status='old',form='unformatted') tnum=snum11+snum12+snum13 write(6,*) snum11,snum12,snum13 do it=1,17 !AO do h=hmn,hmx read(57) ((coldiab4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((cocdiab4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(57) ((covdiab4(i,j,h),i=imn,imx),j=imn,imx) enddo cotdiab4 = coldiab4 + cocdiab4 + covdiab4 h = hlev do j = jmn, jmx do i = imn, imx coldiab3(i,j,it) = coldiab4(i,j,h) cocdiab3(i,j,it) = cocdiab4(i,j,h) covdiab3(i,j,it) = covdiab4(i,j,h) cotdiab3(i,j,it) = cotdiab4(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 coldiab(i,j,it) = coldiab(i,j,it) + coldiab3(i,j,it)*snum/tnum cocdiab(i,j,it) = cocdiab(i,j,it) + cocdiab3(i,j,it)*snum/tnum covdiab(i,j,it) = covdiab(i,j,it) + covdiab3(i,j,it)*snum/tnum cotdiab(i,j,it) = cotdiab(i,j,it) + cotdiab3(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_diabatic_rotate_DRATE05_100_'//CAT//'_'//SEA//'_'//jetpos2//'_160E_190E.data',status='old',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) ((coldiab4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((cocdiab4(i,j,h),i=imn,imx),j=imn,imx) enddo do h=hmn,hmx read(58) ((covdiab4(i,j,h),i=imn,imx),j=imn,imx) enddo cotdiab4 = coldiab4 + cocdiab4 + covdiab4 h = hlev do j = jmn, jmx do i = imn, imx coldiab3(i,j,it) = coldiab4(i,j,h) cocdiab3(i,j,it) = cocdiab4(i,j,h) covdiab3(i,j,it) = covdiab4(i,j,h) cotdiab3(i,j,it) = cotdiab4(i,j,h) enddo enddo do j=jmn,jmx do i=imn,imx coldiab2(i,j,it) = coldiab2(i,j,it) + coldiab3(i,j,it)*snum/tnum2 cocdiab2(i,j,it) = cocdiab2(i,j,it) + cocdiab3(i,j,it)*snum/tnum2 covdiab2(i,j,it) = covdiab2(i,j,it) + covdiab3(i,j,it)*snum/tnum2 cotdiab2(i,j,it) = cotdiab2(i,j,it) + cotdiab3(i,j,it)*snum/tnum2 enddo enddo enddo enddo !season loop it= itime do j=jmn,jmx do i=imn,imx dcou(i,j) = coldiab(i,j,it) - coldiab2(i,j,it) dcov(i,j) = cocdiab(i,j,it) - cocdiab2(i,j,it) dcovor(i,j) = covdiab(i,j,it) - covdiab2(i,j,it) dcotmp(i,j) = cotdiab(i,j,it) - cotdiab2(i,j,it) enddo enddo nu = 0 nv = 0 nvor = 0 ntmp = 0 nspfh = 0 nhgt = 0 do nnt=1,nmx aveu1 = 0. avev1 = 0. avevor1 = 0. avetmp1 = 0. avespfh1 = 0. avehgt1 = 0. aveu2 = 0. avev2 = 0. avevor2 = 0. avetmp2 = 0. avespfh2 = 0. avehgt2 = 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) = ldiab2d(i, j, nn2(it)) sv1(i,j,it) = cdiab2d(i, j, nn2(it)) svor1(i,j,it) = vdiab2d(i, j, nn2(it)) stmp1(i,j,it) = tdiab2d(i, j, nn2(it)) else su2(i,j,it - emem1) = ldiab2d(i, j, nn2(it)) sv2(i,j,it - emem1) = cdiab2d(i, j, nn2(it)) svor2(i,j,it - emem1) = vdiab2d(i, j, nn2(it)) stmp2(i,j,it - emem1) = tdiab2d(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) 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) 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) 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 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)) enddo enddo !psreh=psreh*2. !pcold=pcold*2. !pwarm=pwarm*2. write(6,*) pu(20,20),pv(10,10),pvor(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) endprogram permutation