! $Id: struct_sph.f90,v 1.27 2004/04/21 10:07:17 nilshau Exp $ ! !*********************************************************************** program struct_sph ! ! This program calculates the correlation function along the local magnetic ! field, normal to the local field and binormal to it. ! ! lSF_long=T : first to sixth order magnetic structure functions ! : considering only the part of the vector in the direction ! lSF_long_v=T : first to sixth order kinetic structure functions ! : considering only the part of the vector in the direction ! lSF_dot=T : first to sixth order magnetic structure functions where we ! : the whole vector ! lSF_dot_v=T : first to sixth order magnetic structure functions where we ! : the whole vector ! ! The program reads in the data from ALL the different processors, ! and stores it in a global array. If the runs are to large the use of ! shared memory computers are necessary. In IRIX machines OpenMP is ! used as default, on othere operating systems one can invoke ! OpenMP by adding the line: ! FLAGS_OMP=-mp ! in Makefile.local ! You run the program with the script run_cal.sh. ! !----------------------------------------------------------------------- ! 05-sept-03/nils: adapted from frenet.f90 ! use Cdata use Cparam use Sub use Deriv use General use Mpicomm use Common ! implicit none ! ! ! define parameters ! character (len=130) :: file,file_ext,input_file character (len=4) :: ch integer :: cpu,mom,exp1,exp2,lb_ll,counter real :: dummy real :: x1p,y1p,step_length real, dimension ((maxstep+1-minstep)**3) :: step_length_arr real, dimension ((maxstep+1-minstep)**3,3) :: step_arr real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (:,:,:,:), allocatable :: bb_loc,fluxp_loc,fluxm_loc real, dimension (:,:,:,:), allocatable :: bb_glob,vv_glob real, dimension (:,:,:,:), allocatable :: fluxp_glob,fluxm_glob real, dimension (nx,3) :: bb,BnB real, dimension (nx) :: BnB2,l_d,Cx,Cy,Cz,bx2 double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_long,SF_long_local,SF_long_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_long_v,SF_long_v_local,SF_long_v_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_trans,SF_trans_local,SF_trans_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_trans_v,SF_trans_v_local,SF_trans_v_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_dot,SF_dot_local,SF_dot_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_dot_v,SF_dot_v_local,SF_dot_v_local_mm double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: & SF_flux,SF_flux_local,SF_flux_local_mm real, dimension (:,:), allocatable :: B0,B_vec,tmpB,V_vec,V0 real, dimension (:,:), allocatable :: fluxp0,fluxm0,fluxp_vec,fluxm_vec real, dimension (:), allocatable :: sf1 integer, dimension (:), allocatable :: z1p,zzz,zzp1 integer :: l1_glob,l2_glob,m1_glob,m2_glob,n1_glob,n2_glob integer :: xxp1,yyp1 integer :: ll,mmm,nnn,i,dir_nr,dir_count,ii integer :: subbox_meshpoints,sub_l,sub_m,sub_n,slp,smp,snp,sln,smn,snn integer :: xstepp,xstep,ystep,zstep,step2,nr_dir logical :: lSF_long,lSF_long_v,lext,lSF_dot,lSF_dot_v,olddir,newdir logical :: lSF_trans,lSF_trans_v,lSF_flux logical :: lmag=.false. logical :: lkin=.false. logical :: lflux=.false. ! ! Stuff for OpenMP ! integer :: my_number=0,total_number=1,ipy_loc,ipz_loc integer :: ny_OMP,nz_OMP,ipy_OMP,ipz_OMP,yprocs,zprocs,m_OMP,n_OMP ! ! ! Setting some constants. (Should have been set by invoking the different ! initilize routines.) ! iux=1 iuy=2 iuz=3 ilnrho=4 iax=5 iay=6 iaz=7 ! ip=10 ! ! Reading frenet.in ! call safe_character_assign(input_file,'struct_sph.in') if (ip<8) print*,'Reading '//trim(input_file) open(1,FILE=input_file,FORM='formatted') read(1,*) snap_start,snap_stop read(1,*) lSF_long,lSF_long_v,lSF_trans,lSF_trans_v,lSF_dot,lSF_dot_v read(1,*) lSF_flux close(1) ! ! Allocating allocatable arrays. The large global arrays must be ! defined with allocate and NOT in the normal way, in order to ! make sure they go into the heap and not the stack. This is ! important on shared memory computers where there normally is ! an upper limit on the stack but not on the heap. ! !Allocating array for velocities and B-field only when necessary if (lSF_dot .or. lSF_long .or. lSF_trans) lmag=.true. if (lSF_dot_v .or. lSF_long_v .or. lSF_trans_v) lkin=.true. if (lSF_flux) then lflux=.true. lmag=.true. endif if (lmag) allocate(bb_glob(nxgrid,nygrid,nzgrid,3)) if (lkin) allocate(vv_glob(nxgrid,nygrid,nzgrid,3)) if (lflux) allocate(fluxp_glob(nxgrid,nygrid,nzgrid,3)) if (lflux) allocate(fluxm_glob(nxgrid,nygrid,nzgrid,3)) ! ! Loop over snapshots ! do snap=snap_start,snap_stop ! ! Loop over number of CPU's to get global array ! ! Beginning parallel session ! !$OMP PARALLEL DEFAULT(none) & !$OMP firstprivate(datadir,snap_stop,ip,lshear,lmag,lkin,lflux) & !$OMP firstprivate(lSF_long,lSF_long_v,lSF_dot,lSF_dot_v) & !$OMP firstprivate(lSF_trans,lSF_trans_v,lSF_flux) & !$OMP firstprivate(total_number,yprocs,zprocs,ny_OMP,nz_OMP) & !$OMP firstprivate(snap_start,snap) & !$OMP private(my_number,cpu,ipy,ipz,ch) & !$OMP private(file,f,x,y,z,deltay,ipy_loc,ipz_loc) & !$OMP private(f,bb,bb_loc,fluxp_loc,fluxm_loc) & !$OMP private(l1_glob,l2_glob,m1_glob,m2_glob) & !$OMP private(n1_glob,n2_glob) & !$OMP private(counter) & !$OMP private(SF_long_local_mm,SF_long_v_local_mm,V_vec,V0) & !$OMP private(SF_trans_local_mm,SF_trans_v_local_mm) & !$OMP private(SF_flux_local,SF_flux_local_mm,fluxp_vec,fluxm_vec) & !$OMP private(mmm,nnn,ll,B0,i,ii,zzz,fluxp0,fluxm0) & !$OMP private(B_vec,m_OMP,n_OMP,nr_dir,xstepp) & !$OMP private(x1p,y1p,z1p,xxp1,yyp1,zzp1) & !$OMP private(ipy_OMP,ipz_OMP) & !$OMP private(SF_long_local,SF_long_v_local,mom,sf1,exp1,exp2,lb_ll) & !$OMP private(SF_dot_local,SF_dot_local_mm,tmpB) & !$OMP private(SF_trans_local,SF_trans_v_local) & !$OMP private(SF_dot_v_local,SF_dot_v_local_mm) & !$OMP private(sub_l,sub_m,sub_n,bx2,xstep,ystep,zstep,step_length) & !$OMP private(slp,smp,snp,sln,smn,snn,dummy,step2,olddir,newdir) & !$OMP shared(iaz,iay,iax,iux,iuy,iuz,csnap) & !$OMP shared(bb_glob,t,dx,dy,dz,SF_long,SF_long_v,vv_glob) & !$OMP shared(SF_trans,SF_trans_v,SF_flux,fluxp_glob,fluxm_glob) & !$OMP shared(SF_dot,SF_dot_v,nr_dir_glob,step_length_arr,step_arr) ! !$omp barrier ! ! ! Stop if ncpus is not divideable by total_number ! !$ total_number=OMP_GET_NUM_THREADS() if (mod(ncpus,total_number) .ne. 0) then call stop_it('mod(ncpus,total_number) .ne. 0') endif ! ! Determine distribution of processors ! yprocs=1 if (total_number .ge. 4) yprocs=2 if (total_number .ge. 16) yprocs=4 if (total_number .ge. 64) yprocs=8 zprocs=total_number/yprocs if (yprocs*zprocs .ne. total_number) then call stop_it('(yprocs*zprocs) .ne. total_number)') endif ny_OMP=nygrid/yprocs nz_OMP=nzgrid/zprocs ! ! Allocate allocatable arrays ! allocate(fluxp_vec(nz_OMP,3)) allocate(fluxm_vec(nz_OMP,3)) allocate(B0(nz_OMP,3)) allocate(fluxp0(nz_OMP,3)) allocate(fluxm0(nz_OMP,3)) allocate(B_vec(nz_OMP,3)) allocate(tmpB(nz_OMP,3)) allocate(V_vec(nz_OMP,3)) allocate(V0(nz_OMP,3)) allocate(z1p(nz_OMP)) allocate(zzp1(nz_OMP)) allocate(zzz(nz_OMP)) allocate(sf1(nz_OMP)) if (lmag) allocate(bb_loc(mx,my,mz,3)) if (lflux) allocate(fluxp_loc(mx,my,mz,3)) if (lflux) allocate(fluxm_loc(mx,my,mz,3)) ! ! Who am I? ! !$ my_number=OMP_GET_THREAD_NUM() ! ipy_OMP = modulo(my_number, yprocs) ipz_OMP = my_number/(yprocs) ! ! Read data from all processors ! do ipz_loc=0,nprocz/zprocs-1 do ipy_loc=0,nprocy/yprocs-1 ipy=ipy_OMP*nprocy/yprocs+ipy_loc ipz=ipz_OMP*nprocz/zprocs+ipz_loc cpu=ipy+ipz*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) print*,'Reading ',trim(file) ! Use 'my_number' instead of e.g. 1 as file identifier. open(my_number,file=file,form='unformatted') read(my_number) f(1:mx,1:my,1:mz,1:mvar) if (lshear) then read(my_number) t,x,y,z,dx,dy,dz,deltay else read(my_number) t,x,y,z,dx,dy,dz endif close(my_number) ! ! Find magnetic field (if necessary) ! if (lmag) then do m_OMP=m1,m2 do n_OMP=n1,n2 call curl_OMP(f,iax,bb,m_OMP,n_OMP) bb_loc(l1:l2,m_OMP,n_OMP,:)=bb enddo enddo endif ! ! Find flux (if necessary) ! if (lflux) then fluxp_loc(l1:l2,m1:m2,n1:n2,:)=f(l1:l2,m1:m2,n1:n2,iux:iuz) & +bb_loc(l1:l2,m1:m2,n1:n2,:) fluxm_loc(l1:l2,m1:m2,n1:n2,:)=f(l1:l2,m1:m2,n1:n2,iux:iuz) & -bb_loc(l1:l2,m1:m2,n1:n2,:) endif ! ! Put local variables in global array ! l1_glob=1 l2_glob=nxgrid m1_glob=1+ipy*ny m2_glob=(1+ipy)*ny n1_glob=1+ipz*nz n2_glob=(1+ipz)*nz ! if (lmag) then bb_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & bb_loc(l1:l2,m1:m2,n1:n2,:) endif if (lkin) then vv_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & f(l1:l2,m1:m2,n1:n2,iux:iuz) endif if (lflux) then fluxp_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & fluxp_loc(l1:l2,m1:m2,n1:n2,iux:iuz) fluxm_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & fluxm_loc(l1:l2,m1:m2,n1:n2,iux:iuz) endif ! ! Finnishing reading data from all processors ! enddo enddo ! ! --------------------------------------------------- ! Start calculation of the functions of interest ! --------------------------------------------------- ! SF_long=0 SF_long_local=0 SF_long_v=0 SF_long_v_local=0 SF_dot=0 SF_dot_local=0 SF_dot_v=0 SF_dot_v_local=0 SF_flux=0 SF_flux_local=0 ! ! Define which part of the global array to work on (In case of OpenMP) ! l1_glob=1 l2_glob=nxgrid m1_glob=1+ipy_OMP*ny_OMP m2_glob=(1+ipy_OMP)*ny_OMP n1_glob=1+ipz_OMP*nz_OMP n2_glob=(1+ipz_OMP)*nz_OMP ! ! Need a barrier such that the whole global array is produced ! before we start the calculations ! !$omp barrier ! if (my_number .eq. 0) then print*,'------------------------------------' print*,'Start calculation of structure functions' print*,'------------------------------------' endif ! ! Define zzz (should probably also do the same with xxx and yyy) ! do ii=1,nz_OMP; zzz(ii)=n1_glob-1+ii; enddo ! ! Loop over all individual directions ! nr_dir=0 do xstepp=minstep,maxstep do ystep=minstep,maxstep do zstep=minstep,maxstep ! ! If maxstep is zero; calculate stucture functions in the ! x direction ! if (maxstep == 0) then xstep=1 else xstep=xstepp end if ! ! Make sure not to calculate along the same direction several times ! (e.g. 1,1,1 and 2,2,2) ! newdir=.true. if ((xstep .eq. 0).and.(ystep .eq. 0).and.(zstep .eq. 0)) then newdir=.false. endif do step2=2,maxstep if (mod(xstep,step2)==0 .and. mod(ystep,step2)==0 .and. & mod(zstep,step2)==0) newdir=.false. enddo ! ! If this is a new direction; start calculation ! if (newdir) then nr_dir=nr_dir+1 ! ! Find length of step ! step_length=sqrt(xstep**2.+ystep**2.+zstep**2.) !$OMP master step_length_arr(nr_dir)=step_length step_arr(nr_dir,1)=xstep step_arr(nr_dir,2)=ystep step_arr(nr_dir,3)=zstep print*,'xstep,ystep,zstep=',xstep,ystep,zstep !$OMP end master ! ! Loop over all grid points in *my* part of the global array ! do ll=l1_glob,l2_glob if ((my_number .eq. 0) .and. (mod(ll,10).eq.0)) print*,'ll,l2_glob=',ll,l2_glob ! if (my_number .eq. 0) print*,'ll,l2_glob=',ll,l2_glob SF_long_local_mm(:,:,nr_dir) = 0 SF_long_v_local_mm(:,:,nr_dir) = 0 SF_trans_local_mm(:,:,nr_dir) = 0 SF_trans_v_local_mm(:,:,nr_dir)= 0 SF_dot_local_mm(:,:,nr_dir) = 0 SF_dot_v_local_mm(:,:,nr_dir) = 0 SF_flux_local_mm(:,:,nr_dir) = 0 do mmm=m1_glob,m2_glob if (lmag) B0 = bb_glob(ll,mmm,n1_glob:n2_glob,:) if (lkin) V0 = vv_glob(ll,mmm,n1_glob:n2_glob,:) if (lflux) fluxp0 = fluxp_glob(ll,mmm,n1_glob:n2_glob,:) if (lflux) fluxm0 = fluxm_glob(ll,mmm,n1_glob:n2_glob,:) ! ! Loop over all separations ! counter=0 do lb_ll=0,lb_nxgrid*2-2 counter=counter+1 exp2=mod((lb_ll),2) if (lb_ll .eq. 1) exp2=0 exp1=int((lb_ll)/2)-exp2 i=(2**exp1)*(3**exp2) if (lb_ll .eq. 0) i=0 ! ! Find coordinates ! x1p= ll+i*xstep y1p=mmm+i*ystep z1p=zzz+i*zstep ! ! wrap coordinates back into array if outside ! xxp1=mod(int(x1p+3*nx-1),nx)+1 yyp1=mod(int(y1p+3*nx-1),nx)+1 zzp1=mod(int(z1p+3*nx-1),nx)+1 ! ! Field at grid point. ! if (lmag) B_vec = bb_glob(xxp1,yyp1,zzp1,:) if (lkin) V_vec = vv_glob(xxp1,yyp1,zzp1,:) if (lflux) fluxp_vec = fluxp_glob(xxp1,yyp1,zzp1,:) if (lflux) fluxm_vec = fluxm_glob(xxp1,yyp1,zzp1,:) ! ! Add up the results ! ! Magnetic structure functions. ! Only considering component in the direction ! of interest. (We do get very much noise for the higher ! order SF for large resolutions. Don't know why.) if (lSF_long) call struct(SF_long_local_mm,xstep,ystep,& zstep,B_vec,B0,counter,nr_dir,nz_OMP,.false.,.false.) ! Kinetic structure function. Only considering component ! in the direction of interest. if (lSF_long_v) call struct(SF_long_v_local_mm,xstep,ystep,& zstep,V_vec,V0,counter,nr_dir,nz_OMP,.false.,.false.) ! Magnetic structure functions. ! Only considering component transversal to the direction ! of interest. (We do get very much noise for the higher ! order SF for large resolutions. Don't know why.) if (lSF_trans) call struct(SF_trans_local_mm,xstep,ystep,& zstep,B_vec,B0,counter,nr_dir,nz_OMP,.false.,.true.) ! Kinetic structure function. Only considering component ! transversal to the direction of interest. if (lSF_trans_v) call struct(SF_trans_v_local_mm,xstep,ystep,& zstep,V_vec,V0,counter,nr_dir,nz_OMP,.false.,.true.) ! Magnetic structure function. Here we ! consider the whole vector and not just the part along ! the direction of interest. if (lSF_dot) call struct(SF_dot_local_mm,xstep,ystep,& zstep,B0-B_vec,V0,counter,nr_dir,nz_OMP,.true.,.false.) ! Kinetic second order structure function. Consider the ! whole vector as for lSF_dot. if (lSF_dot_v) call struct(SF_dot_v_local_mm,xstep,ystep,& zstep,V0-V_vec,V0,counter,nr_dir,nz_OMP,.true.,.false.) ! Flux like structure function. (Politano & Pouquet) if (lSF_flux) call struct_flux(SF_flux_local_mm,xstep,ystep,& zstep,fluxp_vec,fluxp0,fluxm_vec,fluxm0,counter,nr_dir,nz_OMP) ! ! end do end do SF_long_local(:,:,nr_dir) = SF_long_local(:,:,nr_dir) & + SF_long_local_mm(:,:,nr_dir) SF_long_v_local(:,:,nr_dir) = SF_long_v_local(:,:,nr_dir) & + SF_long_v_local_mm(:,:,nr_dir) SF_trans_local(:,:,nr_dir) = SF_trans_local(:,:,nr_dir) & + SF_trans_local_mm(:,:,nr_dir) SF_trans_v_local(:,:,nr_dir) = SF_trans_v_local(:,:,nr_dir) & + SF_trans_v_local_mm(:,:,nr_dir) SF_dot_local(:,:,nr_dir) = SF_dot_local(:,:,nr_dir) & + SF_dot_local_mm(:,:,nr_dir) SF_dot_v_local(:,:,nr_dir) = SF_dot_v_local(:,:,nr_dir) & + SF_dot_v_local_mm(:,:,nr_dir) SF_flux_local(:,:,nr_dir) = SF_flux_local(:,:,nr_dir) & + SF_flux_local_mm(:,:,nr_dir) end do ! ! Summing contribution from all threads ! !$OMP critical if (lSF_long) SF_long(:,:,nr_dir) = SF_long(:,:,nr_dir) & + SF_long_local(:,:,nr_dir) if (lSF_long_v) SF_long_v(:,:,nr_dir) = SF_long_v(:,:,nr_dir) & + SF_long_v_local(:,:,nr_dir) if (lSF_trans) SF_trans(:,:,nr_dir) = SF_trans(:,:,nr_dir) & + SF_trans_local(:,:,nr_dir) if (lSF_trans_v) SF_trans_v(:,:,nr_dir) = SF_trans_v(:,:,nr_dir) & + SF_trans_v_local(:,:,nr_dir) if (lSF_dot) SF_dot(:,:,nr_dir) = SF_dot(:,:,nr_dir) & + SF_dot_local(:,:,nr_dir) if (lSF_dot_v) SF_dot_v(:,:,nr_dir) = SF_dot_v(:,:,nr_dir) & + SF_dot_v_local(:,:,nr_dir) if (lSF_flux) SF_flux(:,:,nr_dir) = SF_flux(:,:,nr_dir) & + SF_flux_local(:,:,nr_dir) !$OMP end critical ! ! Finish loop over all directions ! endif enddo enddo enddo nr_dir_glob=nr_dir ! ! Deallocating ! deallocate(B0) deallocate(B_vec) deallocate(tmpB) deallocate(V_vec) deallocate(V0) deallocate(z1p) deallocate(zzp1) deallocate(zzz) deallocate(sf1) ! deallocate(fluxp_vec) deallocate(fluxm_vec) deallocate(fluxp0) deallocate(fluxm0) if (lmag) deallocate(bb_loc) if (lflux) deallocate(fluxp_loc) if (lflux) deallocate(fluxm_loc) ! ! ! Finalizinig parallel session ! !$OMP END PARALLEL ! ! Normalizing ! if (lSF_long) SF_long = SF_long/(nw*ncpus) if (lSF_long_v) SF_long_v = SF_long_v/(nw*ncpus) if (lSF_trans) SF_trans = SF_trans/(nw*ncpus) if (lSF_trans_v) SF_trans_v = SF_trans_v/(nw*ncpus) if (lSF_dot) SF_dot = SF_dot/(nw*ncpus) if (lSF_dot_v) SF_dot_v = SF_dot_v/(nw*ncpus) if (lSF_flux) SF_flux = SF_flux/(nw*ncpus) ! ! Saving the results ! if (lSF_long) & call save_struct('SF_long',step_length_arr,step_arr,SF_long) if (lSF_long_v) & call save_struct('SF_long_v',step_length_arr,step_arr,SF_long_v) if (lSF_trans) & call save_struct('SF_trans',step_length_arr,step_arr,SF_trans) if (lSF_trans_v)& call save_struct('SF_trans_v',step_length_arr,step_arr,SF_trans_v) if (lSF_dot) & call save_struct('SF_dot',step_length_arr,step_arr,SF_dot) if (lSF_dot_v) & call save_struct('SF_dot_v',step_length_arr,step_arr,SF_dot_v) if (lSF_flux) & call save_struct('SF_flux',step_length_arr,step_arr,SF_flux) ! ! Finnishing loop over snapshots (if snap_stop .ne. 0) ! enddo ! ! Deallocating the allocateable arrays ! if (lmag) deallocate(bb_glob) if (lkin) deallocate(vv_glob) if (lflux) deallocate(fluxp_glob) if (lflux) deallocate(fluxm_glob) ! end program struct_sph !**************************************************************************** subroutine save_struct(varr,step_length,step_arr,SF) ! use Cdata use General use Common ! implicit none ! character (len=*), intent(in):: varr character (len=50) :: result_file character (len=30) :: var real, dimension ((maxstep+1-minstep)**3) :: step_length real, dimension ((maxstep+1-minstep)**3,3) :: step_arr double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: SF ! if (snap_stop .eq. 0) call safe_character_assign(result_file, & trim(datadir)//'/'//trim(varr)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(result_file, & trim(datadir)//'/'//trim(varr)//'_VAR'//trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(result_file) open(1,file=result_file) write(1,'(1e12.4)') t write(1,'(1I8)') nr_dir_glob write(1,'(3e12.4)') step_length write(1,'(3e12.4)') step_arr write(1,'(3e12.4)') SF close(1) ! end subroutine save_struct !**************************************************************************** subroutine struct(SF,xstep,ystep,zstep,B_vec,B0,counter,nr_dir,nz_OMP,dot,trans) ! use Common ! implicit none ! integer :: mom,counter,xstep,ystep,zstep,nr_dir,ii,xs,ys,zs,nz_OMP real :: step_length,dotprod double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: SF real, dimension (nz_OMP,3) :: B0,B_vec real, dimension (nz_OMP) :: sf1 logical :: dot,trans ! ! If we want dot product of vector ! if (dot) then sf1=sqrt(B_vec(:,1)**2+B_vec(:,2)**2+B_vec(:,3)**2) do mom=1,maxmom do ii=1,nz_OMP SF(counter,mom,nr_dir)= & SF(counter,mom,nr_dir)+sf1(ii)**mom enddo enddo else ! ! If we want longitudinal or transversal parts ! if (trans) then xs=ystep*(maxstep+2)-zstep ys=zstep-xstep*(maxstep+2) zs=xstep-ystep dotprod=xstep*xs+ystep*ys+zstep*zs if (dotprod .ne. 0) print*,'Error in caclulation of trans. dir.' else xs=xstep ys=ystep zs=zstep endif step_length=sqrt(xs**2.+ys**2.+zs**2.) sf1=abs(xs*B_vec(:,1)+ys*B_vec(:,2)+zs*B_vec(:,3) & -(xs*B0(:,1)+ys*B0(:,2)+zs*B0(:,3)))/step_length do mom=1,maxmom do ii=1,nz_OMP SF(counter,mom,nr_dir)= & SF(counter,mom,nr_dir)+sf1(ii)**mom enddo enddo end if ! end subroutine struct !**************************************************************************** subroutine struct_flux(SF,xstep,ystep,zstep,fluxp_vec,fluxp0,fluxm_vec,fluxm0,counter,nr_dir,nz_OMP) ! use Common ! implicit none ! integer :: mom,counter,xstep,ystep,zstep,nr_dir,ii,xs,ys,zs,nz_OMP real :: step_length,dotprod double precision, dimension (imax,maxmom,(maxstep+1-minstep)**3) :: SF real, dimension (nz_OMP,3) :: fluxp0,fluxp_vec,fluxm0,fluxm_vec,dzm,dzp real, dimension (nz_OMP) :: sf1,dotpart,longpart ! ! Defining some auxillary variables ! dzp=fluxp_vec-fluxp0 dzm=fluxm_vec-fluxm0 xs=xstep ys=ystep zs=zstep step_length=sqrt(xs**2.+ys**2.+zs**2.) ! ! Find longitudinal part and dot product part ! longpart=abs(xs*dzp(:,1)+ys*dzp(:,2)+zs*dzp(:,3))/step_length dotpart=dzm(:,1)**2+dzm(:,2)**2+dzm(:,3)**2 sf1=longpart*dotpart ! ! Small moments ! do mom=1,2 do ii=1,nz_OMP SF(counter,mom,nr_dir)= & SF(counter,mom,nr_dir)+sf1(ii)**(mom/3.) enddo enddo ! ! Third moment ! do ii=1,nz_OMP SF(counter,3,nr_dir)=SF(counter,3,nr_dir)+sf1(ii) enddo ! ! Large moments ! if (maxmom .gt. 3) then do mom=4,maxmom do ii=1,nz_OMP SF(counter,mom,nr_dir)= & SF(counter,mom,nr_dir)+sf1(ii)**(mom/3.) enddo enddo endif ! end subroutine struct_flux !****************************************************************************