! ifort -o qfactor.x qfactor.f90 -fopenmp -mcmodel=medium -O3 -xHost -ipo ! gfortran -o qfactor.x qfactor.f90 -fopenmp -mcmodel=medium -O3 ! -O3, -xHost, -ipo are for a better efficiency ! the efficiency of ifort (version 2021.4.0) is slightly (~1.40 times) faster than gfortran (version 9.3.0) include 'trace_bline.f90' include 'trace_scott.f90' subroutine qfactor0_bridge(i,j) use qfactor_common implicit none real:: bzp, bp(0:2), vp(0:2), rs(0:2), re(0:2), twist0, line_length, q0, q_perp0, incline integer:: i, j, rbs, rbe, s_index, e_index !---------------------------------------------------------------------------- call ij2vp(i, j, vp) call interpolateB(vp, bp) bzp=bp(2) if (scottFlag) then call trace_scott(vp, q0, q_perp0, rs, re, rbs, rbe, line_length, twist0, twistFlag) q(i,j)=q0 q_perp(i,j)=q_perp0 else incline=abs(bzp/norm2(bp)) call trace_bline(vp, rs, re, rbs, rbe, line_length, twist0, twistFlag, incline) endif length(i, j)=line_length if (twistFlag) twist(i, j)=twist0 if( bzp>0.0) then sign2d(i,j)=1.0 reboundary(i, j)=rbe if (vflag) rboundary_tmp(i, j)=1+8*rbe reF(:, i, j)=re if ( (rbe .ne. 0) .and.(rbe .ne. 7)) then call interpolateB(re, bp) e_index=(6-rbe)/2 bnr(i, j)=abs(DBLE(bzp)/bp(e_index)) endif else sign2d(i,j)=-1.0 reboundary(i, j)=rbs if (vflag) rboundary_tmp(i, j)=rbs+8 reF(:, i, j)=rs if ( (rbs .ne. 0) .and.(rbs .ne. 7)) then call interpolateB(rs, bp) s_index=(6-rbs)/2 bnr(i, j)=abs(DBLE(bzp)/bp(s_index)) endif endif if (bzp .eq. 0.0) then sign2d(i,j)=0.0 reboundary(i, j)=0 if (vflag) rboundary_tmp(i, j)=0 endif END subroutine qfactor0_bridge subroutine qfactor0_calculate(i,j) !calculate the Q at the photosphere by Titov (2002), some problematic sites are filled by trace_scott use qfactor_common use trace_common implicit none integer:: i, j, rb, rbs, rbe, m_index, index1, index2, sign_index1, sign_index2 logical:: bkey1, bkey2, bkey11, bkey12, bkey21, bkey22, margin_flag real:: nxx, nxy, nyx, nyy, q0, vp(0:2), q_perp0 real:: rs(0:2), re(0:2), line_length, twist0 !---------------------------------------------------------------------------- rb=reboundary(i,j) if ((rb .eq. 0) .or. (rb .eq. 7) .or. (bnr(i,j) .eq. 0.0)) then q(i,j)=NaN return endif m_index=(6-rb)/2 index1=mod(m_index+1,3) index2=mod(m_index+2,3) margin_flag= (i .eq. 0) .or. (j .eq. 0) .or. (i .eq. q1m1) .or. (j .eq. q2m1) bkey1= (.not. margin_flag) .and. ( rb .eq. reboundary(i+1,j)) .and. ( rb .eq. reboundary(i-1,j)) bkey2= (.not. margin_flag) .and. ( rb .eq. reboundary(i,j+1)) .and. ( rb .eq. reboundary(i,j-1)) bkey11=(i+2 .le. q1m1) .and. ( rb .eq. reboundary(i+1,j)) .and. ( rb .eq. reboundary(i+2,j)) bkey12=(i-2 .ge. 0) .and. ( rb .eq. reboundary(i-1,j)) .and. ( rb .eq. reboundary(i-2,j)) .and. (.not. bkey11) bkey21=(j+2 .le. q2m1) .and. ( rb .eq. reboundary(i,j+1)) .and. ( rb .eq. reboundary(i,j+2)) bkey22=(j-2 .ge. 0) .and. ( rb .eq. reboundary(i,j-1)) .and. ( rb .eq. reboundary(i,j-2)) .and. (.not. bkey21) if (bkey11) sign_index1= 1 if (bkey12) sign_index1=-1 if (bkey21) sign_index2= 1 if (bkey22) sign_index2=-1 if (bkey1) then nxx=reF(index1, i+1,j)-reF(index1, i-1,j) nyx=reF(index2, i+1,j)-reF(index2, i-1,j) else if (bkey11 .or. bkey12) then nxx=(-3*reF(index1, i,j)+4*reF(index1, i+sign_index1,j)-reF(index1, i+2*sign_index1,j))*sign_index1 nyx=(-3*reF(index2, i,j)+4*reF(index2, i+sign_index1,j)-reF(index2, i+2*sign_index1,j))*sign_index1 endif if (bkey2) then nxy=reF(index1, i,j+1)-reF(index1, i,j-1) nyy=reF(index2, i,j+1)-reF(index2, i,j-1) else if (bkey21 .or. bkey22) then nxy=(-3*reF(index1, i,j)+4*reF(index1, i,j+sign_index2)-reF(index1, i,j+2*sign_index2))*sign_index2 nyy=(-3*reF(index2, i,j)+4*reF(index2, i,j+sign_index2)-reF(index2, i,j+2*sign_index2))*sign_index2 endif if ((bkey1 .or. bkey11 .or. bkey12) .and.& (bkey2 .or. bkey21 .or. bkey22)) then q(i,j)=(nxx*nxx+nxy*nxy+nyx*nyx+nyy*nyy) / (bnr(i,j) * (2.0*delta)**2.0) else call ij2vp(i, j, vp) call trace_scott(vp, q0, q_perp0, rs, re, rbs, rbe, line_length, twist0, .false.) q(i,j)=q0 endif end subroutine qfactor0_calculate subroutine qfactor0() use qfactor_common implicit none integer:: i, j !---------------------------------------------------------------------------- allocate(sign2d(0:q1m1, 0:q2m1)) cut_coordinate=0.0 !$OMP PARALLEL DO PRIVATE(i,j), schedule(DYNAMIC) DO j= 0, q2m1 DO i= 0, q1m1 call qfactor0_bridge(i,j) enddo enddo !$OMP END PARALLEL DO if (.not. scottFlag) then !$OMP PARALLEL DO PRIVATE(i,j), schedule(DYNAMIC) DO j= 0, q2m1 DO i= 0, q1m1 call qfactor0_calculate(i,j) enddo enddo !$OMP END PARALLEL DO endif call qmin2() call out_slogq() end subroutine qfactor0 subroutine qcs_bridge(i,j) use qfactor_common use trace_common implicit none real:: bp(0:2), vp(0:2), br(0:2), rs(0:2), re(0:2), & twist0, bn_0, bn_s, bn_e, bnt_tmp, line_length, & vpa1(0:2), rsa1(0:2), rea1(0:2), vpa2(0:2), rsa2(0:2), rea2(0:2), & vpb1(0:2), rsb1(0:2), reb1(0:2), vpb2(0:2), rsb2(0:2), reb2(0:2), & nxx, nxy, nyx, nyy, sxx, sxy, syx, syy, exx, exy, eyx, eyy, & q0, q_perp0, incline, u0(0:2), v0(0:2) integer:: i, j, rbs, rbe, & rbsa1, rbea1, rbsa2, rbea2, rbsb1, rbeb1, rbsb2, rbeb2, & maxdim, index1, index2, pmin_mark(0:2), pmax_mark(0:2), & s_index0, s_index1, s_index2, e_index0, e_index1, e_index2, index_a, index_b logical:: surface_flag, edge_flag, bkey !---------------------------------------------------------------------------- call ij2vp(i, j, vp) if (scottFlag) then call trace_scott(vp, q0, q_perp0, rs, re, rbs, rbe, line_length, twist0, twistFlag) q(i,j)=q0 q_perp(i,j)=q_perp0 else call interpolateB(vp, bp) if (csflag) then incline=abs(dot_product(ev3,bp)/ norm2(bp)) else incline=abs(bp(Normal_index)/ norm2(bp)) endif call trace_bline(vp, rs, re, rbs, rbe, line_length, twist0, twistFlag, incline) endif rsF(:, i, j)=rs reF(:, i, j)=re rsboundary( i, j)=rbs reboundary( i, j)=rbe length(i, j)=line_length if (twistFlag) twist(i, j)=twist0 if (scottFlag) return if ((rbe .eq. 0) .or. (rbs .eq. 0) .or. (rbe .eq. 7) .or. (rbs .eq. 7)) return !---------------------------------------------------------------------------- ! deal with field lines touching the cut plane e_index0=(6-rbe)/2 s_index0=(6-rbs)/2 call interpolateB(re, br) bn_e=br(e_index0) call interpolateB(rs, br) bn_s=br(s_index0) pmin_mark=0 pmax_mark=0 where(vp .lt. (pmin+delta) ) pmin_mark=1 where(vp .gt. (pmax-delta) ) pmax_mark=1 select case( sum(pmin_mark+pmax_mark)) case(0) surface_flag=.false. edge_flag=.false. case(1) surface_flag=.true. edge_flag=.false. case(2:3) surface_flag=.true. edge_flag=.true. tangent_Flag(i, j)=.true.! unreal, for preventing bug q(i,j)=NaN end select if (csflag) then bn_0=sum(DBLE(bp)*ev3) else bn_0=bp(Normal_index) endif bnr(i, j)=abs(DBLE(bn_s)*bn_e/(DBLE(bn_0)**2.)) if (surface_flag) then maxdim=sum(maxloc(pmin_mark+pmax_mark))-1 else maxdim=sum(maxloc(abs(bp)))-1 endif if ( ((incline .le. min_incline) .or. surface_flag) & .and. (maxdim .ne. Normal_index) .and. (.not. edge_flag)) then tangent_Flag(i, j)=.true. index1=mod(maxdim+1,3) index2=mod(maxdim+2,3) if (surface_flag) then u0=[0.,0.,0.] u0(index1)=1. v0=[0.,0.,0.] v0(index2)=1. incline=abs(bp(maxdim)/ norm2(bp)) bn_0=bp(maxdim) else !use the plane perpendicular to the field line bn_0=norm2(bp) v0(maxdim)= bp(index1) v0(index1)=-bp(maxdim) v0(index2)=0. v0 =v0 - dot_product(v0, bp)*bp /(bn_0*bn_0) v0 =v0/norm2(v0) u0(0)=dble(bp(1))*v0(2)-dble(bp(2))*v0(1) u0(1)=dble(bp(2))*v0(0)-dble(bp(0))*v0(2) u0(2)=dble(bp(0))*v0(1)-dble(bp(1))*v0(0) u0 =u0/norm2(u0) incline=1. endif bnt_tmp=abs(bn_s*bn_e/(dble(bn_0)**2.)) vpa1=vp+delta*u0 call trace_bline(vpa1, rsa1, rea1, rbsa1, rbea1, line_length, twist0, .false., incline) vpa2=vp-delta*u0 call trace_bline(vpa2, rsa2, rea2, rbsa2, rbea2, line_length, twist0, .false., incline) vpb1=vp+delta*v0 call trace_bline(vpb1, rsb1, reb1, rbsb1, rbeb1, line_length, twist0, .false., incline) vpb2=vp-delta*v0 call trace_bline(vpb2, rsb2, reb2, rbsb2, rbeb2, line_length, twist0, .false., incline) bkey= ( rbs .eq. rbsa1) .and. ( rbs .eq. rbsa2) .and. & ( rbs .eq. rbsb1) .and. ( rbs .eq. rbsb2) .and. & ( rbe .eq. rbea1) .and. ( rbe .eq. rbea2) .and. & ( rbe .eq. rbeb1) .and. ( rbe .eq. rbeb2) s_index1=mod(s_index0+1,3) s_index2=mod(s_index0+2,3) e_index1=mod(e_index0+1,3) e_index2=mod(e_index0+2,3) if (bkey) then sxx = rsa1(s_index1)-rsa2(s_index1) syx = rsa1(s_index2)-rsa2(s_index2) exx = rea1(e_index1)-rea2(e_index1) eyx = rea1(e_index2)-rea2(e_index2) sxy = rsb1(s_index1)-rsb2(s_index1) syy = rsb1(s_index2)-rsb2(s_index2) exy = reb1(e_index1)-reb2(e_index1) eyy = reb1(e_index2)-reb2(e_index2) nxx = exx*syy - exy*syx nxy = -exx*sxy + exy*sxx nyx = eyx*syy - eyy*syx nyy = -eyx*sxy + eyy*sxx qtmp(i,j) = (nxx*nxx + nxy*nxy + nyx*nyx + nyy*nyy) * bnt_tmp /((2.*delta)**4.) else call trace_scott(vp, q0, q_perp0, rs, re, rbs, rbe, line_length, twist0, .false.) qtmp(i,j)=q0 endif endif END subroutine qcs_bridge subroutine qcs_calculate(i,j) use qfactor_common use trace_common implicit none integer:: i, j, rbs, rbe integer:: s_index0, s_index1, s_index2, sign_s_index1, sign_s_index2, e_index0, e_index1, e_index2, sign_e_index1, sign_e_index2 logical:: bkeys1, bkeys2, bkeys11, bkeys12, bkeys21, bkeys22, bkeye1,bkeye2, bkeye11, bkeye12, bkeye21, bkeye22 logical:: margin_flag1, margin_flag2 real:: sxx, sxy, syx, syy, exx, exy, eyx, eyy, nxx, nxy, nyx, nyy, vp(0:2), q0, q_perp0, rs(0:2), re(0:2), line_length, twist0 !---------------------------------------------------------------------------- rbs=rsboundary(i,j) rbe=reboundary(i,j) if((rbs .eq. 0) .or. (rbe .eq. 0) .or. (rbs .eq. 7) .or. (rbe .eq. 7)) then q(i,j)=NaN return endif if (tangent_Flag(i,j)) then q(i,j)=qtmp(i, j) return endif margin_flag1= (i .eq. 0) .or. (i .eq. q1m1) margin_flag2= (j .eq. 0) .or. (j .eq. q2m1) bkeys1 = (.not. margin_flag1) .and. ( rbs .eq. rsboundary(i+1,j)) .and. ( rbs .eq. rsboundary(i-1,j)) bkeys2 = (.not. margin_flag2) .and. ( rbs .eq. rsboundary(i,j+1)) .and. ( rbs .eq. rsboundary(i,j-1)) bkeye1 = (.not. margin_flag1) .and. ( rbe .eq. reboundary(i+1,j)) .and. ( rbe .eq. reboundary(i-1,j)) bkeye2 = (.not. margin_flag2) .and. ( rbe .eq. reboundary(i,j+1)) .and. ( rbe .eq. reboundary(i,j-1)) bkeye11=(i+2 .le. q1m1) .and. ( rbe .eq. reboundary(i+1,j)) .and. ( rbe .eq. reboundary(i+2,j)) bkeye12=(i-2 .ge. 0) .and. ( rbe .eq. reboundary(i-1,j)) .and. ( rbe .eq. reboundary(i-2,j)) .and. (.not. bkeye11) bkeye21=(j+2 .le. q2m1) .and. ( rbe .eq. reboundary(i,j+1)) .and. ( rbe .eq. reboundary(i,j+2)) bkeye22=(j-2 .ge. 0) .and. ( rbe .eq. reboundary(i,j-1)) .and. ( rbe .eq. reboundary(i,j-2)) .and. (.not. bkeye21) bkeys11=(i+2 .le. q1m1) .and. ( rbs .eq. rsboundary(i+1,j)) .and. ( rbs .eq. rsboundary(i+2,j)) bkeys12=(i-2 .ge. 0) .and. ( rbs .eq. rsboundary(i-1,j)) .and. ( rbs .eq. rsboundary(i-2,j)) .and. (.not. bkeys11) bkeys21=(j+2 .le. q2m1) .and. ( rbs .eq. rsboundary(i,j+1)) .and. ( rbs .eq. rsboundary(i,j+2)) bkeys22=(j-2 .ge. 0) .and. ( rbs .eq. rsboundary(i,j-1)) .and. ( rbs .eq. rsboundary(i,j-2)) .and. (.not. bkeys21) s_index0=(6-rbs)/2 s_index1=mod(s_index0+1,3) s_index2=mod(s_index0+2,3) e_index0=(6-rbe)/2 e_index1=mod(e_index0+1,3) e_index2=mod(e_index0+2,3) if (bkeys11) sign_s_index1= 1 if (bkeys12) sign_s_index1=-1 if (bkeys21) sign_s_index2= 1 if (bkeys22) sign_s_index2=-1 if (bkeye11) sign_e_index1= 1 if (bkeye12) sign_e_index1=-1 if (bkeye21) sign_e_index2= 1 if (bkeye22) sign_e_index2=-1 if (bkeys1) then sxx = rsF(s_index1, i+1, j)-rsF(s_index1, i-1, j) syx = rsF(s_index2, i+1, j)-rsF(s_index2, i-1, j) else if ((bkeys11) .or. (bkeys12)) then sxx = (-3*rsF(s_index1, i, j)+4*rsF(s_index1, i+sign_s_index1, j)-rsF(s_index1, i+2*sign_s_index1, j))*sign_s_index1 syx = (-3*rsF(s_index2, i, j)+4*rsF(s_index2, i+sign_s_index1, j)-rsF(s_index2, i+2*sign_s_index1, j))*sign_s_index1 endif if (bkeys2) then sxy = rsF(s_index1, i, j+1)-rsF(s_index1, i, j-1) syy = rsF(s_index2, i, j+1)-rsF(s_index2, i, j-1) else if ((bkeys21) .or. (bkeys22)) then sxy = (-3*rsF(s_index1, i, j)+4*rsF(s_index1, i, j+sign_s_index2)-rsF(s_index1, i, j+2*sign_s_index2))*sign_s_index2 syy = (-3*rsF(s_index2, i, j)+4*rsF(s_index2, i, j+sign_s_index2)-rsF(s_index2, i, j+2*sign_s_index2))*sign_s_index2 endif if (bkeye1) then exx = reF(e_index1, i+1, j)-reF(e_index1, i-1, j) eyx = reF(e_index2, i+1, j)-reF(e_index2, i-1, j) else if ((bkeye11) .or. (bkeye12)) then exx = (-3*reF(e_index1, i, j)+4*reF(e_index1, i+sign_e_index1, j)-reF(e_index1, i+2*sign_e_index1, j))*sign_e_index1 eyx = (-3*reF(e_index2, i, j)+4*reF(e_index2, i+sign_e_index1, j)-reF(e_index2, i+2*sign_e_index1, j))*sign_e_index1 endif if (bkeye2) then exy = reF(e_index1, i, j+1)-reF(e_index1, i, j-1) eyy = reF(e_index2, i, j+1)-reF(e_index2, i, j-1) else if ((bkeye21) .or. (bkeye22)) then exy = (-3*reF(e_index1, i, j)+4*reF(e_index1, i, j+sign_e_index2)-reF(e_index1, i, j+2*sign_e_index2))*sign_e_index2 eyy = (-3*reF(e_index2, i, j)+4*reF(e_index2, i, j+sign_e_index2)-reF(e_index2, i, j+2*sign_e_index2))*sign_e_index2 endif if ((bkeys1 .or. bkeys11 .or. bkeys12) .and. & (bkeys2 .or. bkeys21 .or. bkeys22) .and. & (bkeye1 .or. bkeye11 .or. bkeye12) .and. & (bkeye2 .or. bkeye21 .or. bkeye22)) then nxx = exx*syy - exy*syx nxy = -exx*sxy + exy*sxx nyx = eyx*syy - eyy*syx nyy = -eyx*sxy + eyy*sxx q(i,j) = (nxx*nxx + nxy*nxy + nyx*nyx + nyy*nyy) * bnr(i,j) / ((2.*delta)**4.) else call ij2vp(i, j, vp) call trace_scott(vp, q0, q_perp0, rs, re, rbs, rbe, line_length, twist0, .false.) q(i,j)=q0 endif end subroutine qcs_calculate subroutine qcs() use qfactor_common implicit none integer:: i,j !---------------------------------------------------------------------------- tangent_Flag=.false. !$OMP PARALLEL DO PRIVATE(i,j), schedule(DYNAMIC) DO j= 0, q2m1 DO i= 0, q1m1 call qcs_bridge(i,j) enddo enddo !$OMP END PARALLEL DO if (.not. scottFlag) then !$OMP PARALLEL DO PRIVATE(i,j), schedule(DYNAMIC) DO j= 0, q2m1 DO i= 0, q1m1 call qcs_calculate(i,j) enddo enddo !$OMP END PARALLEL DO endif call qmin2() end subroutine qcs subroutine ij2vp(i, j, vp) use qfactor_common implicit none integer:: i, j real:: vp(0:2) !---------------------------------------------------------------------------- select case(Normal_index) case(-1) vp=point0+dble(i*delta)*ev1+dble(j*delta)*ev2 case(0) vp=[cut_coordinate , i*delta+yreg(0), j*delta+zreg(0)] case(1) vp=[i*delta+xreg(0), cut_coordinate, j*delta+zreg(0)] case(2) vp=[i*delta+xreg(0), j*delta+yreg(0), cut_coordinate] end select end subroutine ij2vp subroutine qmin2() use qfactor_common implicit none integer:: i, j !---------------------------------------------------------------------------- DO j= 0, q2m1 DO i= 0, q1m1 if (.not. ISNAN(q(i,j)) .and. (q(i,j).lt. 2.0)) q(i,j)=2.0 enddo enddo if (scottFlag) then DO j= 0, q2m1 DO i= 0, q1m1 if (.not. ISNAN(q_perp(i,j)) .and. (q_perp(i,j) .lt. 2.0)) q_perp(i,j)=2.0 enddo enddo endif end subroutine qmin2 subroutine out_slogq() use qfactor_common implicit none integer:: i, j real, allocatable:: slogq(:,:) !---------------------------------------------------------------------------- allocate(slogq(0:q1m1, 0:q2m1)) slogq=sign2d*log10(q) DO j= 0, q2m1 DO i= 0, q1m1 if (ISNAN(q(i,j))) slogq(i,j)=0.0 enddo enddo open(unit=8, file='slogq.bin', access='stream') write(8) slogq close(8) if (scottFlag) then slogq=sign2d*log10(q_perp) DO j= 0, q2m1 DO i= 0, q1m1 if (ISNAN(q_perp(i,j))) slogq(i,j)=0.0 enddo enddo open(unit=8, file='slogq_perp.bin', access='stream') write(8) slogq close(8) endif deallocate(slogq, sign2d) end subroutine out_slogq subroutine show_time(percent) real:: percent integer:: times(8) !---------------------------------------------------------------------------- call date_and_time(VALUES=times) print 600, percent, times(5), times(6), times(7) 600 format( ' ', F6.2, '% ' ,I2.2, ':', I2.2, ':', I2.2) end subroutine show_time program qfactor use qfactor_common use field_common implicit none integer:: k, k0 !---------------------------------------------------------------------------- call initialize() print*, ' _____________________________________' print*, ' schedule time' call show_time(0.0) !---------------------------------------------------------------------------- if (q0flag) then !z=0 call qfactor0() rboundary_tmp=reboundary(0:q1m1, 0:q2m1) ! For avoiding segmentation fault of IO open(unit=8, file='qfactor0.bin', access='stream') write(8) q, length, Bnr, reF, rboundary_tmp close(8) if (twistFlag) then open(unit=8, file='twist.bin', access='stream') write(8) twist close(8) endif if (scottFlag) then open(unit=8, file='q_perp.bin', access='stream') write(8) q_perp close(8) endif endif !---------------------------------------------------------------------------- if (cflag) then !qcs call qcs() open(unit=8, file='qcs.bin', access='stream') write(8) q, length, rsF, reF rboundary_tmp=rsboundary(0:q1m1, 0:q2m1) ! For avoiding segmentation fault of IO write(8) rboundary_tmp rboundary_tmp=reboundary(0:q1m1, 0:q2m1) ! For avoiding segmentation fault of IO write(8) rboundary_tmp close(8) if (twistFlag) then open(unit=8, file='twist.bin', access='stream') write(8) twist close(8) endif if (scottFlag) then open(unit=8, file='q_perp.bin', access='stream') write(8) q_perp close(8) endif endif !---------------------------------------------------------------------------- if (vflag) then !q3d open(unit=1, file='q3d.bin', access='stream') open(unit=2, file='rboundary3d.bin', access='stream') if(scottFlag) open(unit=3, file='q_perp3d.bin', access='stream') if(twistFlag) open(unit=4, file= 'twist3d.bin', access='stream') if (zreg(0) .eq. 0.0) then call qfactor0() write(1) q write(2) rboundary_tmp ! had been transferred to rbs+8*rbe in subroutine qfactor0_bridge if (scottFlag) write(3) q_perp if (twistFlag) write(4) twist k0=1 else k0=0 endif do k=k0, qz-1 cut_coordinate=zreg(0)+k*delta call qcs() rboundary_tmp=rsboundary(0:qx-1, 0:qy-1)+8*reboundary(0:qx-1, 0:qy-1) write(1) q write(2) rboundary_tmp if (scottFlag) write(3) q_perp if (twistFlag) write(4) twist if (mod(k+1, 4) .eq. 0) call show_time(float(k+1)/qz*100.0) enddo close(1) close(2) if (scottFlag) close(3) if (twistFlag) close(4) endif !---------------------------------------------------------------------------- if (.not. (vflag .and. (mod(qz, 4) .eq. 0))) call show_time(100.0) !---------------------------------------------------------------------------- deallocate(Bfield, q, qtmp, reF, length, bnr, reboundary, rboundary_tmp) if (twistFlag) deallocate(twist, curlB) if (scottFlag) deallocate(q_perp, grad_unit_vec_Bfield) if (cflag .or. vflag) deallocate(rsF, rsboundary, tangent_Flag) end program qfactor