! $Id: wavenumbers.f90,v 1.5 2003/12/10 21:23:51 nilshau Exp $ ! !*********************************************************************** program wavenumbers ! ! Program to find lots of diffenrent interesting wavenumbers. ! Be aware that there are some numberes that MUST be specified in ! wavenumbers.in ! !----------------------------------------------------------------------- ! 19-june-03/nils: code ! use Cdata use Cparam use Sub use Deriv use General ! implicit none ! ! ! define parameters ! character (len=130) :: file,wave_out_file,dimfile,wave_in_file character (len=130) :: wave_result_file character (len=4) :: ch,csnap integer :: i,j,k,itx=1,ity=1,itz=1 integer :: cpu real :: urms,brms,eps,epsM,epsK,kf,eta real :: dummy real :: b2_mean,b4_mean,u2_mean,nB2_mean,nu2_mean,BnB2_mean,BxJ2_mean real :: k_rms, k_parallel, k_BxJ,k_lambda,k_nu,k_eta real :: Re,Rm,Re_taylor,Rm_taylor real, dimension (mx,my,mz,mvar+maux) :: f double precision, dimension (my,mz) :: mean_B2,mean_B4,mean_nB2,mean_BnB2 double precision, dimension (my,mz) :: mean_BxJ2,mean_nu2,mean_u2 double precision, dimension (ncpus) :: cpumean_B2,cpumean_B4,cpumean_nB2 double precision, dimension (ncpus) :: cpumean_BnB2 double precision, dimension (ncpus) :: cpumean_BxJ2,cpumean_nu2,cpumean_u2 real, dimension (nx,3) :: bb,BnB,jj,BxJ real, dimension (nx) :: b2,b4,nB2,BnB2,u2,BxJ2 real, dimension (nx) :: bxx,bxy,bxz,byx,byy,byz,bzx,bzy,bzz,a1,a2 integer :: snap,snap_start=0,snap_stop=0 ! Just hardwired since it will rarely be used ! ! Setting some constants ! iux=1 iuy=2 iuz=3 ilnrho=4 iax=5 iay=6 iaz=7 ! ip=10 ! ! ! Loop over snapshots ! do snap=snap_start,snap_stop ! ! Loop over number of CPU's in original run ! do cpu=0,ncpus-1 print*,'Doing processor',cpu+1,'of',ncpus,'processors' ipy = modulo(cpu, nprocy) ipz = cpu/(nprocy) ! ! Read var.dat ! call chn(cpu,ch) call chn(snap,csnap) if (snap_stop .eq. 0) call safe_character_assign(file,trim(datadir) & //'/proc'//trim(ch)//'/var.dat') if (snap_stop .ne. 0) call safe_character_assign(file,trim(datadir) & //'/proc'//trim(ch)//'/VAR'//trim(csnap)) if (ip<8) print*,'Reading '//trim(file) open(1,file=file,form='unformatted') read(1) f(1:mx,1:my,1:mz,1:mvar) if (lshear) then read(1) t,x,y,z,dx,dy,dz,deltay else read(1) t,x,y,z,dx,dy,dz endif close(1) ! ! Read urms, brms, etc...... ! call safe_character_assign(wave_in_file,'wavenumbers.in') if (ip<8) print*,'Reading '//trim(wave_in_file) open(1,FILE=wave_in_file,FORM='formatted') read(1,*) urms,brms read(1,*) eps,epsK,epsM read(1,*) nu,eta read(1,*) kf close(1) ! ! Loop over all pencils ! do m=m1,m2 do n=n1,n2 ! ! Some usefull variables ! call curl(f,iax,bb) call dot2_mn(bb,b2) mean_B2(m,n)=sum(b2)/nx b4=b2*b2 mean_b4(m,n)=sum(b4)/nx call derij(f,iaz,a1,2,1) call derij(f,iay,a2,3,1) bxx=a1-a2 ! bxx=xder(yder(f(:,:,:,iaz)))-xder(zder(f(:,:,:,iay))) call der2(f,iaz,a1,2) call derij(f,iay,a2,3,2) ! bxy=yder2(f(:,:,:,iaz))-yder(zder(f(:,:,:,iay))) bxy=a1-a2 call derij(f,iaz,a1,2,3) call der2(f,iay,a2,3) ! bxz=zder(yder(f(:,:,:,iaz)))-zder2(f(:,:,:,iay)) bxz=a1-a2 call der2(f,iaz,a1,1) call derij(f,iax,a2,3,1) byx=-a1+a2 ! byx=-xder2(f(:,:,:,iaz))+xder(zder(f(:,:,:,iax))) call derij(f,iaz,a1,1,2) call derij(f,iax,a2,3,2) byy=-a1+a2 ! byy=-yder(xder(f(:,:,:,iaz)))+yder(zder(f(:,:,:,iax))) call derij(f,iaz,a1,1,3) call der2(f,iax,a2,3) byz=-a1+a2 ! byz=-zder(xder(f(:,:,:,iaz)))+zder2(f(:,:,:,iax)) call der2(f,iay,a1,1) call derij(f,iax,a2,2,1) bzx=a1-a2 ! bzx=xder2(f(:,:,:,iay))-xder(yder(f(:,:,:,iax))) call derij(f,iay,a1,1,2) call der2(f,iax,a2,2) bzy=a1-a2 ! bzy=yder(xder(f(:,:,:,iay)))-yder2(f(:,:,:,iax)) call derij(f,iay,a1,1,3) call derij(f,iax,a2,2,3) bzz=a1-a2 !bzz=zder(xder(f(:,:,:,iay)))-zder(yder(f(:,:,:,iax))) ! ! Finding k_rms2 ! nB2=bxx**2.+bxy**2.+bxz**2.+byx**2.+byy**2.+byz**2.+bzx**2.+bzy**2.+bzz**2. mean_nB2(m,n)=sum(nB2)/nx ! ! Finding k_parallel ! BnB(:,1)=bb(:,1)*bxx+bb(:,2)*bxy+bb(:,3)*bxz BnB(:,2)=bb(:,1)*byx+bb(:,2)*byy+bb(:,3)*byz BnB(:,3)=bb(:,1)*bzx+bb(:,2)*bzy+bb(:,3)*bzz call dot2_mn(BnB,BnB2) mean_BnB2(m,n) = sum(BnB2)/nx ! ! Finding k_BxJ ! call del2v_etc(f,iax,curlcurl=jj) call cross_mn(bb,jj,BxJ) call dot2_mn(BxJ,BxJ2) mean_BxJ2(m,n)=sum(BxJ2)/nx ! ! Finnishing magnetic stuff ! Starting kinetic stuff ! ! ! Some usefull variables ! call dot2(f(:,:,:,iux:iuz),u2) call der_main(f,iux,bxx,1) call der_main(f,iux,bxy,2) call der_main(f,iux,bxz,3) !bxx=xder(f(:,:,:,iux)) !bxy=yder(f(:,:,:,iux)) !bxz=zder(f(:,:,:,iux)) call der_main(f,iuy,byx,1) call der_main(f,iuy,byy,2) call der_main(f,iuy,byz,3) !byx=xder(f(:,:,:,iuy)) !byy=yder(f(:,:,:,iuy)) !byz=zder(f(:,:,:,iuy)) call der_main(f,iuz,bzx,1) call der_main(f,iuz,bzy,2) call der_main(f,iuz,bzz,3) !bzx=xder(f(:,:,:,iuz)) !bzy=yder(f(:,:,:,iuz)) !bzz=zder(f(:,:,:,iuz)) ! ! Finding k_lambda ! nB2=bxx**2.+bxy**2.+bxz**2.+byx**2.+byy**2.+byz**2.+bzx**2.+bzy**2.+bzz**2. mean_nu2(m,n)=sum(nB2)/nx mean_u2(m,n) =sum(u2)/nx ! ! Finnishing loop over pencils ! enddo enddo ! ! Finding mean values ! b2_mean=sum(mean_b2(m1:m2,n1:n2))/(ny*nz) b4_mean=sum(mean_b4(m1:m2,n1:n2))/(ny*nz) u2_mean=sum(mean_u2(m1:m2,n1:n2))/(ny*nz) nB2_mean=sum(mean_nB2(m1:m2,n1:n2))/(ny*nz) nu2_mean=sum(mean_nu2(m1:m2,n1:n2))/(ny*nz) BnB2_mean=sum(mean_BnB2(m1:m2,n1:n2))/(ny*nz) BxJ2_mean=sum(mean_BxJ2(m1:m2,n1:n2))/(ny*nz) ! ! Saving data for each processor ! call safe_character_assign(wave_out_file,trim(datadir)//'/proc'//trim(ch)//'/wave.dat') if (ip<8) print*,'Writing '//trim(wave_out_file) open(91,file=wave_out_file,form='unformatted') write(91) B2_mean,nB2_mean,BnB2_mean,B4_mean,BxJ2_mean,u2_mean,nu2_mean close(91) ! ! Finnishing loop over processors ! enddo ! ! Calculating the interesting wavenumbers ! do cpu=0,ncpus-1 !print*,'Doing processor',cpu+1,'of',ncpus,'processors' ! ! Read wave.dat ! call chn(cpu,ch) call safe_character_assign(file,trim(datadir)//'/proc'//trim(ch)//'/wave.dat') if (ip<8) print*,'Reading '//trim(file) open(1,file=file,form='unformatted') read(1) B2_mean,nB2_mean,BnB2_mean,B4_mean,BxJ2_mean,u2_mean,nu2_mean close(1) cpumean_b2(cpu+1)=B2_mean cpumean_b4(cpu+1)=B4_mean cpumean_nB2(cpu+1)=nB2_mean cpumean_BnB2(cpu+1)=BnB2_mean cpumean_u2(cpu+1)=u2_mean cpumean_nu2(cpu+1)=nu2_mean cpumean_BxJ2(cpu+1)=BxJ2_mean enddo ! ! Finding global means ! B2_mean=sum(cpumean_b2)/ncpus B4_mean=sum(cpumean_b4)/ncpus u2_mean=sum(cpumean_u2)/ncpus nu2_mean=sum(cpumean_nu2)/ncpus BnB2_mean=sum(cpumean_BnB2)/ncpus nB2_mean=sum(cpumean_nB2)/ncpus BxJ2_mean=sum(cpumean_BxJ2)/ncpus ! ! Calculating, saving and printing wavenumbers and Reynolds numbers ! k_rms=sqrt(nB2_mean/B2_mean) k_parallel=sqrt(BnB2_mean/B4_mean) k_BxJ=sqrt(BxJ2_mean/B4_mean) k_lambda=sqrt(nu2_mean/u2_mean) k_nu=(epsK/nu**3.)**(1./4) k_eta=(epsM/eta**3.)**(1./4) Re=urms/(kf*nu) Rm=urms/(kf*eta) Re_taylor=urms*sqrt(5.)/(nu*k_lambda) Rm_taylor=urms*sqrt(5.)/(eta*k_lambda) ! ! Printing the results ! print*,'k_rms=',k_rms print*,'k_parallel=',k_parallel print*,'k_BxJ=',k_BxJ print*,'k_nu=',k_nu print*,'k_eta=',k_eta print*,'k_lambda=',k_lambda print*,'Re=',Re print*,'Rm=',Rm print*,'Re_taylor=',Re_taylor print*,'Rm_taylor=',Rm_taylor ! ! Saving the results ! if (snap_stop .eq. 0) call safe_character_assign(wave_result_file, & trim(datadir)//'/wave'//'.dat') if (snap_stop .ne. 0) call safe_character_assign(wave_result_file, & trim(datadir)//'/wave_VAR'//trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(wave_result_file) open(1,file=wave_result_file) write(1,'(3e12.4)') k_rms,k_parallel,k_BxJ write(1,'(3e12.4)') k_lambda,k_eta,k_nu write(1,'(4e12.4)') Re,Rm,Re_taylor,Rm_taylor close(1) end do ! end program wavenumbers