! $Id: frenet.f90,v 1.12 2003/12/10 21:23:51 nilshau Exp $ ! !*********************************************************************** program frenet ! ! This program calculates the correlation function along the local magnetic ! field, normal to the local field and binormal to it. ! ! lcorr=T : magnetic correlation function ! lcorr_v=T : kinetic correlation function ! lSF=T : first to sixth order magnetic structure functions ! : considering only the part of the vector in the direction ! : of interest (i.e. in the normal, binormal or prallel direction) ! lSF2=T : second order magnetic structure function where we consider ! : the whole vector ! lSF2_v=T : second order kinetic structure function where we consider ! : 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. ! !----------------------------------------------------------------------- ! 22-july-03/nils: coded ! use Cdata use Cparam use Sub use Deriv use General use Mpicomm ! implicit none ! ! ! define parameters ! character (len=130) :: file,corr_result_file,file_ext,input_file character (len=130) :: corr_v_result_file character (len=130) :: SF_result_file,SF2_result_file,SF2_v_result_file character (len=30) :: direction,direction_ character (len=4) :: ch,csnap,file_subbox integer :: cpu,mom,exp1,exp2,lb_ll,counter integer, parameter :: max_moment=6 integer, parameter :: imax=lb_nxgrid*2-1 real :: dummy,sf1 real :: x1p,y1p,z1p real :: x1frac,x2frac,y1frac,y2frac,z1frac,z2frac real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz,3) :: bb_loc ,dir_loc real, dimension (:,:,:,:), allocatable :: bb_glob,dir_glob,vv_glob real, dimension (nx,3) :: bb,BnB,b_parallel,normal,binormal real, dimension (nx) :: b2,BnB2,l_d,Cx,Cy,Cz,b4 real, dimension (nx) :: xxx,xxy,xxz,yyx,yyy,yyz,zzx,zzy,zzz real, dimension (nx) :: bxx,bxy,bxz,byx,byy,byz,bzx,bzy,bzz,a1,a2 double precision, dimension (imax) :: corr,corr_local,corr_local_nn double precision, dimension (imax) :: corr_local_mm double precision, dimension (imax) :: corr_v,corr_v_local,corr_v_local_nn double precision, dimension (imax) :: corr_v_local_mm double precision, dimension (imax,max_moment) :: SF,SF_local,SF_local_nn double precision, dimension (imax,max_moment) :: SF_local_mm double precision, dimension (imax) :: SF2,SF2_local,SF2_local_nn double precision, dimension (imax) :: SF2_local_mm double precision, dimension (imax) :: SF2_v,SF2_v_local,SF2_v_local_nn double precision, dimension (imax) :: SF2_v_local_mm double precision :: meanB2,meanB2_local,meanB2_local_mm,meanB2_local_nn real, dimension (3) :: B0,B_vec,dir0,B_ext,tmpB,V_vec,V0,dir_average integer :: snap,snap_start,snap_stop integer :: l1_glob,l2_glob,m1_glob,m2_glob,n1_glob,n2_glob integer :: xxp1,xxp2,yyp1,yyp2,zzp1,zzp2 integer :: ll,mmm,nnn,i,dir_nr,dir_count integer :: subbox_meshpoints,sub_l,sub_m,sub_n,slp,smp,snp,sln,smn,snn logical :: lSF,lcorr,lcorr_v,lpri,lext,lSF2,lSF2_v ! ! 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 ! ! 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. ! allocate(dir_glob(nxgrid,nygrid,nzgrid,3)) allocate(bb_glob(nxgrid,nygrid,nzgrid,3)) ! ! 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,'frenet.in') if (ip<8) print*,'Reading '//trim(input_file) open(1,FILE=input_file,FORM='formatted') read(1,*) direction_ read(1,*) snap_start,snap_stop read(1,*) lcorr,lcorr_v,lSF,lSF2,lSF2_v read(1,*) B_ext read(1,*) subbox_meshpoints close(1) ! !Allocating array for velocities only when necessary if (lSF2_v .or. lcorr_v) allocate(vv_glob(nxgrid,nygrid,nzgrid,3)) ! ! Doing all directions (parallel, normal and binormal) if direction_='all', ! otherwise doing only direction=direction_ ! if (direction_ .eq. 'all') then dir_count=3 else dir_count=1 endif do dir_nr=1,dir_count if (direction_ .eq. 'all') then if (dir_nr .eq. 1) direction='parallel' if (dir_nr .eq. 2) direction='normal' if (dir_nr .eq. 3) direction='binormal' else direction=direction_ end if ! print*,'Calculating the ',trim(direction),' direction' ! ! 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,direction,ip,lshear) & !$OMP firstprivate(lcorr,lSF,lSF2,B_ext,lSF2_v,lcorr_v) & !$OMP firstprivate(subbox_meshpoints) & !$OMP shared(iaz,iay,iax,iux,iuy,iuz) & !$OMP private(my_number,total_number,cpu,ipy,ipz,ch,snap) & !$OMP private(csnap,file,f,x,y,z,deltay,ipy_loc,ipz_loc) & !$OMP private(f,bb,b2,b4,b_parallel,bb_loc,dir_loc) & !$OMP private(a1,a2,bxx,bxy,bxz,byx,byy,byz,bzx,bzy,bzz ) & !$OMP private(BnB,Cx,Cy,Cz,normal,BnB2,binormal) & !$OMP private(xxx,yyx,zzx,xxy,yyy,zzy,xxz,yyz,zzz) & !$OMP private(l1_glob,l2_glob,m1_glob,m2_glob,corr_local_nn) & !$OMP private(n1_glob,n2_glob,lpri,lext,SF_local_nn) & !$OMP private(SF_local_mm,corr_local_mm,V_vec,V0) & !$OMP private(l_d,corr_local,mmm,nnn,ll,B0,i,x2frac,x1frac) & !$OMP private(y2frac,y1frac,z2frac,z1frac,B_vec,m_OMP,n_OMP) & !$OMP private(x1p,y1p,z1p,xxp1,xxp2,yyp1,yyp2,zzp1,zzp2) & !$OMP private(ny_OMP,nz_OMP,ipy_OMP,ipz_OMP,yprocs,zprocs) & !$OMP private(SF_local,mom,sf1,dir0,exp1,exp2,lb_ll,counter) & !$OMP private(SF2_local,SF2_local_mm,SF2_local_nn,tmpB) & !$OMP private(SF2_v_local,SF2_v_local_mm,SF2_v_local_nn) & !$OMP private(corr_v_local,corr_v_local_mm,corr_v_local_nn) & !$OMP private(meanB2_local,meanB2_local_mm,meanB2_local_nn) & !$OMP private(dir_average,sub_l,sub_m,sub_n) & !$OMP private(slp,smp,snp,sln,smn,snn,dummy) & !$OMP shared(bb_glob,dir_glob,t,dx,dy,dz,corr,SF,SF2,meanB2,SF2_v) & !$OMP shared(vv_glob,corr_v) ! ! Stopping 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 ! ! Who am I? ! !$ my_number=OMP_GET_THREAD_NUM() ! 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 ipy_OMP = modulo(my_number, yprocs) ipz_OMP = my_number/(yprocs) ! ! Read data from all processors ! lext=((B_ext(1).ne.0).or.(B_ext(2).ne.0).or.(B_ext(3).ne.0)) 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) ! ! Finding direction for calculation of the function of interest ! do m_OMP=m1,m2 do n_OMP=n1,n2 lpri=((my_number.eq.0).and.(m_OMP.eq.m1).and.(n_OMP.eq.n1)) call curl_OMP(f,iax,bb,m_OMP,n_OMP) ! ! Adding imposed field (for calculation of direction only) ! if (lext) then if (lpri) print*,'Adding imposed field' bb(:,1)=B_ext(1)+bb(:,1) bb(:,2)=B_ext(2)+bb(:,2) bb(:,3)=B_ext(3)+bb(:,3) endif call dot2_mn(bb,b2) b4=b2*b2 b_parallel(:,1)=bb(:,1)/sqrt(b2) b_parallel(:,2)=bb(:,2)/sqrt(b2) b_parallel(:,3)=bb(:,3)/sqrt(b2) bb_loc(l1:l2,m_OMP,n_OMP,:)=bb if (direction .eq. 'parallel') then dir_loc(l1:l2,m_OMP,n_OMP,:)=b_parallel else if ((direction .eq. 'normal') .or. & (direction .eq. 'binormal') ) then ! ! The normal is in the direction of the curvature, ! i.e. pointing towards the center of the "circle". ! (normal = b_parallel . nabla b_parallel) ! The binormal is the cross product of the local ! magnetic field and the normal. ! (binormal = normal x b_parallel) ! call derij_OMP(f,iaz,a1,2,1,m_OMP,n_omp) call derij_OMP(f,iay,a2,3,1,m_OMP,n_omp) bxx=a1-a2 call der2_OMP(f,iaz,a1,2,m_OMP,n_omp) call derij_OMP(f,iay,a2,3,2,m_OMP,n_omp) bxy=a1-a2 call derij_OMP(f,iaz,a1,2,3,m_OMP,n_omp) call der2_OMP(f,iay,a2,3,m_OMP,n_omp) bxz=a1-a2 call der2_OMP(f,iaz,a1,1,m_OMP,n_omp) call derij_OMP(f,iax,a2,3,1,m_OMP,n_omp) byx=-a1+a2 call derij_OMP(f,iaz,a1,1,2,m_OMP,n_omp) call derij_OMP(f,iax,a2,3,2,m_OMP,n_omp) byy=-a1+a2 call derij_OMP(f,iaz,a1,1,3,m_OMP,n_omp) call der2_OMP(f,iax,a2,3,m_OMP,n_omp) byz=-a1+a2 call der2_OMP(f,iay,a1,1,m_OMP,n_omp) call derij_OMP(f,iax,a2,2,1,m_OMP,n_omp) bzx=a1-a2 call derij_OMP(f,iay,a1,1,2,m_OMP,n_omp) call der2_OMP(f,iax,a2,2,m_OMP,n_omp) bzy=a1-a2 call derij_OMP(f,iay,a1,1,3,m_OMP,n_omp) call derij_OMP(f,iax,a2,2,3,m_OMP,n_omp) bzz=a1-a2 ! BnB(:,1)=(bb(:,1)*bxx+bb(:,2)*bxy+bb(:,3)*bxz)/b2 BnB(:,2)=(bb(:,1)*byx+bb(:,2)*byy+bb(:,3)*byz)/b2 BnB(:,3)=(bb(:,1)*bzx+bb(:,2)*bzy+bb(:,3)*bzz)/b2 ! ! This part is due to derivative of 1/|B| ! Cx=bb(:,1)*bxx+bb(:,2)*byx+bb(:,3)*bzx Cy=bb(:,1)*bxy+bb(:,2)*byy+bb(:,3)*bzy Cz=bb(:,1)*bxz+bb(:,2)*byz+bb(:,3)*bzz ! xxx=bb(:,1)*bb(:,1)*Cx/b4 yyx=bb(:,2)*bb(:,1)*Cy/b4 zzx=bb(:,3)*bb(:,1)*Cz/b4 ! xxy=bb(:,1)*bb(:,2)*Cx/b4 yyy=bb(:,2)*bb(:,2)*Cy/b4 zzy=bb(:,3)*bb(:,2)*Cz/b4 ! xxz=bb(:,1)*bb(:,3)*Cx/b4 yyz=bb(:,2)*bb(:,3)*Cy/b4 zzz=bb(:,3)*bb(:,3)*Cz/b4 ! BnB(:,1)=BnB(:,1)-(xxx+yyx+zzx) BnB(:,2)=BnB(:,2)-(xxy+yyy+zzy) BnB(:,3)=BnB(:,3)-(xxz+yyz+zzz) ! call dot2_mn(BnB,BnB2) normal(:,1)=BnB(:,1)/sqrt(BnB2) normal(:,2)=BnB(:,2)/sqrt(BnB2) normal(:,3)=BnB(:,3)/sqrt(BnB2) if (direction .eq. 'normal') then dir_loc(l1:l2,m_OMP,n_omp,:)=normal else binormal(:,1)=+(bb(:,2)*normal(:,3) & -bb(:,3)*normal(:,2)) binormal(:,2)=-(bb(:,1)*normal(:,3) & -bb(:,3)*normal(:,1)) binormal(:,3)=+(bb(:,1)*normal(:,2) & -bb(:,2)*normal(:,1)) call dot2_mn(binormal,l_d) dir_loc(l1:l2,m_OMP,n_omp,1)=binormal(:,1)/sqrt(l_d) dir_loc(l1:l2,m_OMP,n_omp,2)=binormal(:,2)/sqrt(l_d) dir_loc(l1:l2,m_OMP,n_omp,3)=binormal(:,3)/sqrt(l_d) endif else print*,'No such direction:',direction endif enddo enddo ! ! 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 ! dir_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & dir_loc(l1:l2,m1:m2,n1:n2,:) bb_glob(l1_glob:l2_glob,m1_glob:m2_glob,n1_glob:n2_glob,:)= & bb_loc(l1:l2,m1:m2,n1:n2,:) if (lSF2_v .or. lcorr_v) 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 ! ! Finnishing loop over processors ! enddo enddo ! ! From now on we only consider the 'fluctuations' of B, ! we therefore subtract the imposed field. ! if (lext .and. (my_number .eq. 0)) then print*,'Removing imposed field' bb_glob(:,:,:,1)=bb_glob(:,:,:,1)-B_ext(1) bb_glob(:,:,:,2)=bb_glob(:,:,:,2)-B_ext(2) bb_glob(:,:,:,3)=bb_glob(:,:,:,3)-B_ext(3) endif ! ! --------------------------------------------------- ! Start calculation of the functions of interest ! --------------------------------------------------- ! corr=0 corr_v=0 corr_local=0 corr_v_local=0 SF=0 SF_local=0 SF2=0 SF2_local=0 SF2_v=0 SF2_v_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 ' if (lcorr) print*,'correlation functions' if (lcorr .and. lSF) print*,'and' if (lSF) print*,'structure functions' print*,'------------------------------------' endif ! ! Loop over all grid points in *my* part of the global array. ! do ll=l1_glob,l2_glob if (my_number .eq. 0) print*,'Doing ',ll,' of ',nxgrid SF_local_mm=0 corr_local_mm=0 corr_v_local_mm=0 SF2_local_mm=0 SF2_v_local_mm=0 meanB2_local_mm=0 do mmm=m1_glob,m2_glob SF_local_nn=0 SF2_local_nn=0 SF2_v_local_nn=0 meanB2_local_nn=0 corr_local_nn=0 corr_v_local_nn=0 do nnn=n1_glob,n2_glob B0=bb_glob(ll,mmm,nnn,:) if (lSF2_v .or. lcorr_v) V0=vv_glob(ll,mmm,nnn,:) ! ! 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 average direction of small sub box ! dir_average=dir_glob(ll,mmm,nnn,:) if (subbox_meshpoints .gt. 0) then do sub_l=0,subbox_meshpoints do sub_m=0,subbox_meshpoints do sub_n=0,subbox_meshpoints slp=mod( ll+sub_l+3*nx-1,nx)+1 smp=mod(mmm+sub_m+3*nx-1,nx)+1 snp=mod(nnn+sub_n+3*nx-1,nx)+1 sln=mod( ll-sub_l+3*nx-1,nx)+1 smn=mod(mmm-sub_m+3*nx-1,nx)+1 snn=mod(nnn-sub_n+3*nx-1,nx)+1 dir_average=dir_average+ & dir_glob(slp,smp,snp,:)+ & dir_glob(sln,smn,snn,:) enddo enddo enddo dummy=dir_average(1)*dir_average(1)+ & dir_average(2)*dir_average(2)+ & dir_average(3)*dir_average(3) dir_average(1)=dir_average(1)/sqrt(dummy) dir_average(2)=dir_average(2)/sqrt(dummy) dir_average(3)=dir_average(3)/sqrt(dummy) ! dir_average=dir_average/((2*subbox_meshpoints+1)**3) endif ! ! Find fractional coordinates (for interpolation) ! x1p= ll+i*dir_average(1) y1p=mmm+i*dir_average(2) z1p=nnn+i*dir_average(3) ! ! Find integer coordinates (for interpolation) ! xxp1=mod(int(x1p+3*nx-1),nx)+1 xxp2=mod(xxp1+1+nx-1,nx)+1 yyp1=mod(int(y1p+3*nx-1),nx)+1 yyp2=mod(yyp1+1+nx-1,nx)+1 zzp1=mod(int(z1p+3*nx-1),nx)+1 zzp2=mod(zzp1+1+nx-1,nx)+1 ! ! Determine weight for each of the eight ! interpolation points. ! x2frac=mod(x1p+real(nx)-1,real(nx))+1 - real(xxp1) x1frac=1.-x2frac y2frac=mod(y1p+real(nx)-1,real(nx))+1 - real(yyp1) y1frac=1.-y2frac z2frac=mod(z1p+real(nx)-1,real(nx))+1 - real(zzp1) z1frac=1.-z2frac ! ! Use interpolation to find field at fractional ! grid point. ! B_vec= x1frac*y1frac*z1frac*bb_glob(xxp1,yyp1,zzp1,:)+ & x1frac*y1frac*z2frac*bb_glob(xxp1,yyp1,zzp2,:)+ & x1frac*y2frac*z1frac*bb_glob(xxp1,yyp2,zzp1,:)+ & x1frac*y2frac*z2frac*bb_glob(xxp1,yyp2,zzp2,:)+ & x2frac*y1frac*z1frac*bb_glob(xxp2,yyp1,zzp1,:)+ & x2frac*y1frac*z2frac*bb_glob(xxp2,yyp1,zzp2,:)+ & x2frac*y2frac*z1frac*bb_glob(xxp2,yyp2,zzp1,:)+ & x2frac*y2frac*z2frac*bb_glob(xxp2,yyp2,zzp2,:) if (lSF2_v .or. lcorr_v) then V_vec= x1frac*y1frac*z1frac*vv_glob(xxp1,yyp1,zzp1,:)+& x1frac*y1frac*z2frac*vv_glob(xxp1,yyp1,zzp2,:)+ & x1frac*y2frac*z1frac*vv_glob(xxp1,yyp2,zzp1,:)+ & x1frac*y2frac*z2frac*vv_glob(xxp1,yyp2,zzp2,:)+ & x2frac*y1frac*z1frac*vv_glob(xxp2,yyp1,zzp1,:)+ & x2frac*y1frac*z2frac*vv_glob(xxp2,yyp1,zzp2,:)+ & x2frac*y2frac*z1frac*vv_glob(xxp2,yyp2,zzp1,:)+ & x2frac*y2frac*z2frac*vv_glob(xxp2,yyp2,zzp2,:) endif ! ! Add up the results ! ! Magnetic correlation function if (lcorr) corr_local_nn(counter)=corr_local_nn(counter)+ & B0(1)*B_vec(1)+B0(2)*B_vec(2)+B0(3)*B_vec(3) ! Kinetic correlation function if (lcorr_v) corr_v_local_nn(counter)= & corr_v_local_nn(counter)+ & V0(1)*V_vec(1)+V0(2)*V_vec(2)+V0(3)*V_vec(3) ! Magnetic structure functions from first to sixth ! order. 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) then dir0=dir_average sf1=dir0(1)*B_vec(1)+dir0(2)*B_vec(2)+dir0(3)*B_vec(3) & -(dir0(1)*B0(1)+dir0(2)*B0(2)+dir0(3)*B0(3)) do mom=1,max_moment SF_local_nn(counter,mom)= & SF_local_nn(counter,mom)+abs(sf1)**mom enddo endif ! Magnetic second order structure function. Here we ! consider the whole vector and not just the part along ! the direction of interest. if (lSF2) then tmpB=B0-B_vec SF2_local_nn(counter)=SF2_local_nn(counter)+ & tmpB(1)**2+tmpB(2)**2+tmpB(3)**2 endif ! Kinetic second order structure function. Consider the ! whole vector as for lSF2. if (lSF2_v) then tmpB=V0-V_vec SF2_v_local_nn(counter)=SF2_v_local_nn(counter)+ & tmpB(1)**2+tmpB(2)**2+tmpB(3)**2 endif ! end do ! Calculating the mean of B^2. if (lSF2) meanB2_local_nn=meanB2_local_nn+ & B0(1)**2+B0(2)**2+B0(3)**2 end do meanB2_local_mm=meanB2_local_mm+meanB2_local_nn SF_local_mm=SF_local_mm+SF_local_nn SF2_local_mm=SF2_local_mm+SF2_local_nn SF2_v_local_mm=SF2_v_local_mm+SF2_v_local_nn corr_local_mm=corr_local_mm+corr_local_nn corr_v_local_mm=corr_v_local_mm+corr_v_local_nn end do meanB2_local=meanB2_local+meanB2_local_mm SF_local=SF_local+SF_local_mm SF2_local=SF2_local+SF2_local_mm SF2_v_local=SF2_v_local+SF2_v_local_mm corr_local=corr_local+corr_local_mm corr_v_local=corr_v_local+corr_v_local_mm end do ! ! Summing contribution from all threads ! !$OMP critical if (lcorr) corr=corr+corr_local if (lcorr_v) corr_v=corr_v+corr_v_local if (lSF2) meanB2=meanB2+meanB2_local if (lSF) SF=SF+SF_local if (lSF2) SF2=SF2+SF2_local if (lSF2_v) SF2_v=SF2_v+SF2_v_local !$OMP end critical ! ! Finalizinig parallel session ! !$OMP END PARALLEL ! ! Normalizing ! if (lcorr) corr=corr/corr(1) if (lcorr_v) corr_v=corr_v/corr_v(1) if (lSF) SF=SF/(nw*ncpus) if (lSF2) SF2=SF2/(nw*ncpus) if (lSF2_v) SF2_v=SF2_v/(nw*ncpus) if (lSF2) meanB2=meanB2/(nw*ncpus) ! ! Saving the results ! if (direction .eq. 'parallel') then file_ext='parallel' else if (direction .eq. 'normal') then file_ext='normal' else if (direction .eq. 'binormal') then file_ext='binormal' else print*,'No such direction:',direction endif call chn(subbox_meshpoints,file_subbox) ! ! ! if (lcorr) then if (snap_stop .eq. 0) call safe_character_assign(corr_result_file, & trim(datadir)//'/corr_'//trim(file_ext)//trim(file_subbox)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(corr_result_file, & trim(datadir)//'/corr_'//trim(file_ext)//trim(file_subbox)//'_VAR'// & trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(corr_result_file) open(1,file=corr_result_file) write(1,'(3e12.4)') t write(1,'(3e12.4)') corr close(1) endif if (lcorr_v) then if (snap_stop .eq. 0) call safe_character_assign(corr_v_result_file, & trim(datadir)//'/corr_v_'//trim(file_ext)//trim(file_subbox)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(corr_v_result_file, & trim(datadir)//'/corr_v_'//trim(file_ext)//trim(file_subbox)//'_VAR'// & trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(corr_v_result_file) open(1,file=corr_v_result_file) write(1,'(3e12.4)') t write(1,'(3e12.4)') corr_v close(1) endif if (lSF) then if (snap_stop .eq. 0) call safe_character_assign(SF_result_file, & trim(datadir)//'/SF_'//trim(file_ext)//trim(file_subbox)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(SF_result_file, & trim(datadir)//'/SF_'//trim(file_ext)//trim(file_subbox)//'_VAR'// & trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(SF_result_file) open(1,file=SF_result_file) write(1,'(3e12.4)') t write(1,'(3e12.4)') SF close(1) endif if (lSF2) then if (snap_stop .eq. 0) call safe_character_assign(SF2_result_file, & trim(datadir)//'/SF2_'//trim(file_ext)//trim(file_subbox)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(SF2_result_file, & trim(datadir)//'/SF2_'//trim(file_ext)//trim(file_subbox)//'_VAR'// & trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(SF2_result_file) open(1,file=SF2_result_file) write(1,'(3e12.4)') t,meanB2 write(1,'(3e12.4)') SF2 close(1) endif if (lSF2_v) then if (snap_stop .eq. 0) call safe_character_assign(SF2_v_result_file,& trim(datadir)//'/SF2_v_'//trim(file_ext)//trim(file_subbox)//'.dat') if (snap_stop .ne. 0) call safe_character_assign(SF2_v_result_file,& trim(datadir)//'/SF2_v_'//trim(file_ext)//trim(file_subbox)//'_VAR'// & trim(csnap)//'.dat') close(1) if (ip<8) print*,'Writing '//trim(SF2_v_result_file) open(1,file=SF2_v_result_file) write(1,'(3e12.4)') t write(1,'(3e12.4)') SF2_v close(1) endif ! ! Finnishing loop over snapshots (if snap_stop .ne. 0) ! enddo ! ! Finnishing loop over directions (if direction_ .eq. 'all') ! enddo ! ! Deallocating the allocateable arrays ! deallocate(dir_glob) deallocate(bb_glob) ! end program frenet