! $Id$ ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lspecial = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 0 ! !*************************************************************** ! !--------------------------------------------------------------- ! ! HOW TO USE THIS FILE ! !--------------------------------------------------------------- ! ! The rest of this file may be used as a template for your own ! special module. Lines which are double commented are intended ! as examples of code. Simply fill out the prototypes for the ! features you want to use. ! ! Save the file with a meaningful name, eg. geo_kws.f90 and place ! it in the $PENCIL_HOME/src/special directory. This path has ! been created to allow users ot optionally check their contributions ! in to the Pencil-Code CVS repository. This may be useful if you ! are working on/using the additional physics with somebodyelse or ! may require some assistance from one of the main Pencil-Code team. ! ! To use your additional physics code edit the Makefile.local in ! the src directory under the run directory in which you wish to ! use your additional physics. Add a line with all the module ! selections to say something like: ! ! SPECIAL=special/nstar ! ! Where nstar it replaced by the filename of your new module ! upto and not including the .f90 ! !-------------------------------------------------------------------- ! module Special ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! include '../special.h' ! Global arrays real, dimension (nx,nz) :: uu_average real, dimension (nx,nz) :: phidot_average ! "pencils" real, dimension (nx,3) :: uuadvec_guu,uuadvec_gaa real, dimension (nx) :: uuadvec_grho,uuadvec_glnrho,uuadvec_gss real, dimension (nx) :: uu_residual ! logical :: lno_radial_advection=.false. logical :: lfargoadvection_as_shift=.false. logical :: lkeplerian_gauge=.false. logical :: lremove_volume_average=.false. ! real, dimension (nxgrid) :: xgrid1 real :: nygrid1 ! namelist /special_run_pars/ lno_radial_advection, lfargoadvection_as_shift,& lkeplerian_gauge,lremove_volume_average ! integer :: idiag_nshift=0 ! contains ! !*********************************************************************** subroutine register_special() ! ! Configure pre-initialised (i.e. before parameter read) variables ! which should be know to be able to evaluate ! ! 6-oct-03/tony: coded ! use Cdata ! if (lroot) call svn_id( & "$Id$") ! endsubroutine register_special !*********************************************************************** subroutine initialize_special(f) ! ! called by run.f90 after reading parameters, but before the time loop ! ! 06-oct-03/tony: coded ! use Mpicomm, only: stop_it ! real, dimension (mx,my,mz,mvar+maux) :: f ! ! Make it possible to switch the algorithm off while still ! having this file compiled, for debug purposes. ! if (lrun .and. .not. lfargo_advection) then if (lroot) then print*,"" print*,"Switch" print*," lfargo_advection=T" print*,"in init_pars of start.in if you" print*,"want to use the fargo algorithm" print*,"" endif call warning("initialize_special","") endif ! ! Not implemented for other than cylindrical coordinates ! if (lfargo_advection.and.coord_system/='cylindric') then if (lroot) then print*,"" print*,"Fargo advection is only implemented for" print*,"cylindrical coordinates. Switch" print*," coord_system='cylindric'" print*,"in init_pars of start.in if you" print*,"want to use the fargo algorithm" print*,"" endif call fatal_error("initialize_special","") endif ! ! Not implemented for the energy equation either ! if (lfargo_advection.and.(pretend_lnTT.or.ltemperature)) & call fatal_error("initialize_special","fargo advection not "//& "implemented for the temperature equation") ! ! Stuff that is only calculated once ! xgrid1=1./xgrid nygrid1=1./nygrid ! call keep_compiler_quiet(f) ! endsubroutine initialize_special !*********************************************************************** subroutine init_special(f) ! ! initialise special condition; called from start.f90 ! 06-oct-2003/tony: coded ! use Cdata use Mpicomm ! real, dimension (mx,my,mz,mvar+maux) :: f ! intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine init_special !*********************************************************************** subroutine pencil_criteria_special() ! ! All pencils that this special module depends on are specified here. ! ! 18-07-06/tony: coded ! if (lfargo_advection) then ! lpenc_requested(i_uu)=.true. ! ! For continuity equation ! lpenc_requested(i_divu)=.true. if (ldensity_nolog) then lpenc_requested(i_grho)=.true. else lpenc_requested(i_glnrho)=.true. endif ! ! For velocity advection ! lpenc_requested(i_uij)=.true. ! ! For the induction equation ! if (lmagnetic) then lpenc_requested(i_aa)=.true. lpenc_requested(i_aij)=.true. endif ! ! For the entropy equation ! if (lentropy) lpenc_requested(i_gss)=.true. ! endif ! endsubroutine pencil_criteria_special !*********************************************************************** subroutine calc_pencils_special(f,p) ! use Sub, only: h_dot_grad ! real, dimension(mx,my,mz,mfarray) :: f real, dimension(nx,3) :: uu_advec,tmp2 real, dimension(nx) :: tmp type (pencil_case) :: p integer :: j,nnghost ! if (lfargo_advection) then ! nnghost=n-nghost ! ! Advect by the relative velocity ! uu_residual=p%uu(:,2)-uu_average(:,nnghost) ! ! Advect by the original radial and vertical, but residual azimuthal ! uu_advec(:,1)=p%uu(:,1) uu_advec(:,2)=uu_residual uu_advec(:,3)=p%uu(:,3) ! ! For the continuity equation ! if (ldensity_nolog) then call h_dot_grad(uu_advec,p%grho,uuadvec_grho) else call h_dot_grad(uu_advec,p%glnrho,uuadvec_glnrho) endif ! ! For velocity advection ! ! Note: It is tempting to use ! ! call h_dot_grad(uu_advec,p%uij,p%uu,uuadvec_guu) ! ! instead of the lines coded below, but the line just above ! would introduce the curvature terms with the residual ! velocity. Yet the curvature terms do not enter the ! non-coordinate basis advection that fargo performs. ! These terms have to be added manually. If one uses ! h_dot_grad_vec, then the curvature with residual would ! have to be removed manually and then the full speed ! added again (for the r-component of uuadvec_guu), like ! this: ! ! uuadvec_guu(:,1)=uuadvec_guu(:,1)-& ! rcyl_mn1*((p%uu(:,2)-uu_advec(:,2))*p%uu(:,2)) ! uuadvec_guu(:,2)=uuadvec_guu(:,2)+& ! rcyl_mn1*((p%uu(:,1)-uu_advec(:,1))*p%uu(:,2)) ! ! ! Although working (and more line-economically), the ! piece of code below is more readable in my opinion. ! do j=1,3 call h_dot_grad(uu_advec,p%uij(:,j,:),tmp) tmp2(:,j)=tmp enddo tmp2(:,1)=tmp2(:,1)-rcyl_mn1*p%uu(:,2)*p%uu(:,2) tmp2(:,2)=tmp2(:,2)+rcyl_mn1*p%uu(:,1)*p%uu(:,2) ! uuadvec_guu=tmp2 ! ! Advection of the magnetic potential ! if (lmagnetic) then do j=1,3 call h_dot_grad(uu_advec,p%aij(:,j,:),tmp) tmp2(:,j)=tmp enddo tmp2(:,1)=tmp2(:,1)-rcyl_mn1*p%uu(:,2)*p%aa(:,2) tmp2(:,2)=tmp2(:,2)+rcyl_mn1*p%uu(:,1)*p%aa(:,2) ! uuadvec_gaa=tmp2 endif ! ! Advection of entropy ! if (lentropy) & call h_dot_grad(uu_advec,p%gss,uuadvec_gss) ! call keep_compiler_quiet(f) ! endif ! endsubroutine calc_pencils_special !*********************************************************************** subroutine dspecial_dt(f,df,p) ! ! Calculate right hand side of ONE OR MORE extra coupled PDEs ! along the 'current' Pencil, i.e. f(l1:l2,m,n) where ! m,n are global variables looped over in equ.f90 ! ! Due to the multi-step Runge Kutta timestepping used one MUST always ! add to the present contents of the df array. NEVER reset it to zero. ! ! Several precalculated Pencils of information are passed if for ! efficiency. ! ! 06-oct-03/tony: coded ! use Cdata use Diagnostics use Mpicomm ! real, dimension (mx,my,mz,mvar+maux) :: f real, dimension (mx,my,mz,mvar) :: df type (pencil_case) :: p integer :: nnghost real, dimension (nx) :: nshift,phidot ! intent(in) :: f,p intent(inout) :: df ! ! Estimate the shift by nshift=phidot*dt/dy ! This is just an approximation, since uavg ! changes from one subtimestep to another, and ! the correct cumulative shift is ! ! nshift=0. ! do itsub=1,3 ! nshift=nshift+phidot*dt_sub/dy ! enddo ! ! But it also works fairly well this way, since uavg ! does not change all that much between subtimesteps. ! if (lfargo_advection) then ! if (ldiagnos) then if (idiag_nshift/=0) then nnghost=n-nghost phidot=uu_average(:,nnghost)*rcyl_mn1 nshift=phidot*dt*dy_1(m) call max_mn_name(nshift,idiag_nshift) endif endif ! call keep_compiler_quiet(f) call keep_compiler_quiet(df) call keep_compiler_quiet(p) ! endif ! endsubroutine dspecial_dt !*********************************************************************** subroutine read_special_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_run_pars, IOSTAT=iostat) ! endsubroutine read_special_run_pars !*********************************************************************** subroutine write_special_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_run_pars) ! endsubroutine write_special_run_pars !*********************************************************************** subroutine rprint_special(lreset,lwrite) ! use Diagnostics ! ! Reads and registers print parameters relevant to special ! ! 06-oct-03/tony: coded ! integer :: iname logical :: lreset,lwr logical, optional :: lwrite ! ! Write information to index.pro ! if (lfargo_advection) then ! lwr = .false. if (present(lwrite)) lwr=lwrite ! if (lreset) then idiag_nshift=0 endif ! do iname=1,nname call parse_name(iname,cname(iname),cform(iname),'nshift',idiag_nshift) enddo ! if (lwr) then write(3,*) 'i_nshift=',idiag_nshift endif ! endif ! endsubroutine rprint_special !*********************************************************************** subroutine special_calc_density(f,df,p) ! ! Calculate a additional 'special' term on the right hand side of the ! mass equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! ! 06-oct-03/tony: coded ! use Cdata use EquationOfState ! real, dimension (mx,my,mz,mvar+maux), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! ! Modified continuity equation ! if (lfargo_advection) then if (ldensity_nolog) then df(l1:l2,m,n,ilnrho) = df(l1:l2,m,n,ilnrho) - & uuadvec_grho - p%rho*p%divu else df(l1:l2,m,n,ilnrho) = df(l1:l2,m,n,ilnrho) - & uuadvec_glnrho - p%divu endif endif ! call keep_compiler_quiet(f) ! endsubroutine special_calc_density !*********************************************************************** subroutine special_calc_hydro(f,df,p) ! ! Calculate a additional 'special' term on the right hand side of the ! momentum equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! use Cdata use Diagnostics ! real, dimension (mx,my,mz,mvar+maux), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! ! Modified momentum equation ! if (lfargo_advection) then df(l1:l2,m,n,iux:iuz)=df(l1:l2,m,n,iux:iuz)-uuadvec_guu ! ! The lines below are not symmetric. This is on purpose, to better ! highlight that fargo advects the azimuthal coordinate ONLY!! ! if (lfirst.and.ldt) then advec_uu=abs(p%uu(:,1)) *dx_1(l1:l2)+ & abs(uu_residual)*dy_1( m )*rcyl_mn1+ & abs(p%uu(:,3)) *dz_1( n ) endif endif ! call keep_compiler_quiet(f) ! endsubroutine special_calc_hydro !*********************************************************************** subroutine special_calc_magnetic(f,df,p) ! ! Calculate a additional 'special' term on the right hand side of the ! induction equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! ! 06-oct-03/tony: coded ! use Cdata ! real, dimension (mx,my,mz,mvar+maux), intent(inout) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! if (lfargo_advection) & df(l1:l2,m,n,iax:iaz)=df(l1:l2,m,n,iax:iaz)-uuadvec_gaa ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine special_calc_magnetic !*********************************************************************** subroutine special_calc_energy(f,df,p) ! ! Calculate a additional 'special' term on the right hand side of the ! energy equation. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! ! 06-oct-03/tony: coded ! use Cdata ! real, dimension (mx,my,mz,mvar+maux), intent(in) :: f real, dimension (mx,my,mz,mvar), intent(inout) :: df type (pencil_case), intent(in) :: p ! if (lfargo_advection) & df(l1:l2,m,n,iss)=df(l1:l2,m,n,iss)-uuadvec_gss ! call keep_compiler_quiet(f) call keep_compiler_quiet(p) ! endsubroutine special_calc_energy !*********************************************************************** subroutine special_before_boundary(f) ! ! Possibility to modify the f array before the boundaries are ! communicated. ! ! Some precalculated pencils of data are passed in for efficiency ! others may be calculated directly from the f array ! ! 06-jul-06/tony: coded ! use Cdata use Mpicomm, only: mpiallreduce_sum ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension (nx,nz) :: fsum_tmp real, dimension (nx) :: uphi integer :: nnghost ! ! Just needs to calculate at the first sub-timestep ! if (lfargo_advection.and.lfirst) then ! ! Pre-calculate the average large scale speed of the flow ! fsum_tmp=0. ! do n=n1,n2;do m=m1,m2 nnghost=n-nghost uphi=f(l1:l2,m,n,iuy) fsum_tmp(:,nnghost)=fsum_tmp(:,nnghost)+uphi*nygrid1 enddo;enddo ! ! The sum has to be done processor-wise ! Sum over processors of same ipz, and different ipy ! --only relevant for 3D, but is here for generality ! call mpiallreduce_sum(fsum_tmp,uu_average,& (/nx,nz/),idir=2) !idir=2 is equal to old LSUMY=.true. endif ! endsubroutine special_before_boundary !*********************************************************************** subroutine special_after_timestep(f,df,dt_sub) ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real :: dt_sub ! if (lmagnetic) then if (lkeplerian_gauge) call apply_keplerian_gauge(f) if (lremove_volume_average) call remove_volume_average(f) endif ! if (lfargo_advection) then if (lfargoadvection_as_shift) then call fourier_shift_fargo(f,df,dt_sub) else if (lroot) then print*,'Fargo advection without Fourier shift' print*,'is not functional yet. Advecting only' print*,'at the last subtimestep was leading to' print*,'errors. Rewrite.' call fatal_error("special_after_timestep","") endif call advect_fargo(f) endif endif ! ! Just for test purposes and comparison with the loop advection ! in Stone, J. et al., JCP 250, 509 (2005) ! if (lno_radial_advection) then f(:,:,:,iux) = 0. df(:,:,:,iux) = 0. endif ! call keep_compiler_quiet(dt_sub) ! endsubroutine special_after_timestep !*********************************************************************** subroutine fourier_shift_fargo(f,df,dt_) ! ! Possibility to modify the f array after the evolution equations ! are solved. ! ! In this case, add the fargo shift to the f and df-array, in ! fourier space. ! ! 06-jul-06/tony: coded ! use Sub use Fourier, only: fft_y_parallel use Cdata use Mpicomm ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real, dimension (nx,ny) :: a_re,a_im real, dimension (nx) :: phidot integer :: ivar,ng real :: dt_ ! ! Pencil uses linear velocity. Fargo will shift based on ! angular velocity. Get phidot from uphi. ! do n=n1,n2 ng=n-n1+1 phidot=uu_average(:,ng)*rcyl_mn1 ! do ivar=1,mvar ! a_re=f(l1:l2,m1:m2,n,ivar); a_im=0. ! ! Forward transform. No need for computing the imaginary part. ! The transform is just a shift in y, so no need to compute ! the x-transform either. ! call fft_y_parallel(a_re,a_im,SHIFT_Y=phidot*dt_,lneed_im=.false.) ! ! Inverse transform of the shifted array back into real space. ! No need again for either imaginary part of x-transform. ! call fft_y_parallel(a_re,a_im,linv=.true.) f(l1:l2,m1:m2,n,ivar)=a_re ! ! Also shift df, unless it is the last subtimestep. ! if (.not.llast) then a_re=df(l1:l2,m1:m2,n,ivar); a_im=0. call fft_y_parallel(a_re,a_im,SHIFT_Y=phidot*dt_,lneed_im=.false.) call fft_y_parallel(a_re,a_im,linv=.true.) df(l1:l2,m1:m2,n,ivar)=a_re endif ! enddo enddo ! endsubroutine fourier_shift_fargo !******************************************************************** subroutine advect_fargo(f) ! ! Possibility to modify the f array after the evolution equations ! are solved. ! ! In this case, do the fargo shift to the f and df-array, in ! real space. ! ! 06-jul-06/tony: coded ! use Sub use Cdata use Mpicomm ! real, dimension (mx,my,mz,mfarray) :: f ! real, dimension (nx,nygrid,mvar) :: faux_remap,faux_remap_shift ! real, dimension (nx,ny) :: faux,tmp2 real, dimension (nx,nygrid) :: tmp ! integer :: ivar,ng,ig,mshift,cellshift,i,mserial ! integer, dimension (nx,nz) :: shift_intg real, dimension (nx,nz) :: shift_total,shift_frac ! ! For shift in real space, the shift is done in integer number of ! cells in the azimuthal direction, so, take only the integer part ! of the velocity for fargo advection. ! do n=1,nz phidot_average(:,n) = uu_average(:,n)*rcyl_mn1 enddo ! ! Define the integer circular shift ! shift_total = phidot_average*dt*dy_1(mpoint) shift_intg = nint(shift_total) ! ! Do circular shift of cells ! do n=n1,n2 ! do ivar=1,mvar faux=f(l1:l2,m1:m2,n,ivar) call remap_to_pencil_y(faux,tmp) faux_remap(:,:,ivar)=tmp enddo ! ng=n-n1+1 ! do i=l1,l2 ig=i-l1+1 cellshift=shift_intg(ig,ng) ! do m=1,ny mserial=m+ipy*ny mshift=mserial-cellshift if (mshift .lt. 1 ) mshift = mshift + nygrid if (mshift .gt. nygrid) mshift = mshift - nygrid ! do ivar=1,mvar faux_remap_shift(ig,mserial,ivar) = faux_remap(ig,mshift,ivar) enddo enddo enddo ! do ivar=1,mvar tmp=faux_remap_shift(:,:,ivar) call unmap_from_pencil_y(tmp, tmp2) f(l1:l2,m1:m2,n,ivar)=tmp2 enddo enddo ! ! Fractional step ! shift_frac = shift_total-shift_intg call fractional_shift(f,shift_frac) ! endsubroutine advect_fargo !******************************************************************** subroutine fractional_shift(f,shift_frac) ! use Deriv, only:der use Mpicomm ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df_frac real, dimension (nx,nz) :: shift_frac,uu_frac real, dimension(3) :: dt_sub real, dimension(nx) :: ushift,dfdy,facx integer :: j,itsubstep ! ! Set boundaries ! if (nprocy==1) then f(:,1:m1-1,:,:) = f(:,m2i:m2,:,:) f(:,m2+1:my,:,:) = f(:,m1:m1i,:,:) else call initiate_isendrcv_bdry(f) call finalize_isendrcv_bdry(f) endif ! dt_sub=dt*beta_ts ! facx=1./(dt*dy_1(mpoint)*rcyl_mn1) ! do n=1,nz uu_frac(:,n)=shift_frac(:,n)*facx enddo ! do itsubstep=1,itorder if (itsubstep==1) then df_frac=0. else df_frac=alpha_ts(itsubstep)*df_frac endif ! do n=n1,n2;do m=m1,m2 ! ushift=uu_frac(:,n-n1+1) ! do j=1,mvar call der(f,j,dfdy,2) df_frac(l1:l2,m,n,j)=df_frac(l1:l2,m,n,j)-ushift*dfdy f(l1:l2,m,n,j) = & f(l1:l2,m,n,j)+ dt_sub(itsubstep)*df_frac(l1:l2,m,n,j) enddo ! enddo;enddo ! enddo ! endsubroutine fractional_shift !******************************************************************** subroutine apply_keplerian_gauge(f) ! use Mpicomm , only: mpiallreduce_sum use Deriv, only: der ! ! Substract mean emf from the radial component of the induction ! equation. Activated only when large Bz fields and are present ! keplerian advection. Due to this u_phi x Bz term, the radial ! component of the magnetic potential ! develops a divergence that grows linearly in time. Since it is ! purely divergent, it is okay analytically. But numerically it leads to ! problems if this divergent grows bigger than the curl, which it does ! eventually. ! ! This is a cylindrical version of the rtime_phiavg special file. ! ! 13-sep-07/wlad: adapted from remove_mean_momenta ! real, dimension (mx,my,mz,mfarray), intent (inout) :: f real, dimension (mx,mz) :: fsum_tmp,glambda_rz real, dimension (mx,my,mz) :: lambda real, dimension (nx) :: glambda_z real :: fac integer :: i ! !if (.not.lupdate_bounds_before_special) then ! print*,'The boundaries have not been updated prior ' ! print*,'to calling this subroutine. This may lead ' ! print*,'to troubles since it needs derivatives ' ! print*,'and integrals, thus properly set ghost zones. ' ! print*,'Use lupdate_bounds_before_special=T in ' ! print*,'the run_pars of run.in.' ! call fatal_error("apply_keplerian_gauge","") !endif ! fac = 1.0/nygrid ! ! Set boundaries of iax ! !call update_ghosts(f,iax) ! ! Average over phi - the result is a (nr,nz) array ! fsum_tmp = 0. do m=m1,m2; do n=1,mz fsum_tmp(:,n) = fsum_tmp(:,n) + fac*f(:,m,n,iax) enddo; enddo ! ! The sum has to be done processor-wise ! Sum over processors of same ipz, and different ipy ! call mpiallreduce_sum(fsum_tmp,glambda_rz,(/mx,mz/),idir=2) ! ! Gauge-transform radial A ! do m=m1,m2 f(l1:l2,m,n1:n2,iax) = f(l1:l2,m,n1:n2,iax) - glambda_rz(l1:l2,n1:n2) enddo ! ! Integrate in R to get lambda, using N=6 composite Simpson's rule. ! Ghost zones in r needed for glambda_r. ! do i=l1,l2 ; do n=1,mz lambda(i,:,n) = dx/6.*( glambda_rz(i-3,n)+glambda_rz(i+3,n)+& 4*(glambda_rz(i-2,n)+glambda_rz(i ,n)+glambda_rz(i+2,n))+& 2*(glambda_rz(i-1,n)+glambda_rz(i+1,n))) enddo; enddo ! ! Gauge-transform vertical A. Ghost zones in z needed for lambda. ! do m=m1,m2; do n=n1,n2 call der(lambda,glambda_z,3) f(l1:l2,m,n,iaz) = f(l1:l2,m,n,iaz) - glambda_z enddo; enddo ! endsubroutine apply_keplerian_gauge !******************************************************************** subroutine remove_volume_average(f) ! use Mpicomm , only: mpiallreduce_sum ! ! Substract mean emf from the radial component of the induction ! equation. Activated only when large Bz fields and are present ! keplerian advection. Due to this u_phi x Bz term, the radial ! component of the magnetic potential ! develops a divergence that grows linearly in time. Since it is ! purely divergent, it is okay analytically. But numerically it leads to ! problems if this divergent grows bigger than the curl, which it does ! eventually. ! ! This is a cylindrical version of the rtime_phiavg special file. ! ! 13-sep-07/wlad: adapted from remove_mean_momenta ! real, dimension (mx,my,mz,mfarray), intent (inout) :: f real :: fsum_tmp,mean_ax,fac integer :: i ! fac = 1.0/nwgrid ! ! Set boundaries of iax ! !call update_ghosts(f,iax) ! ! Average over phi - the result is a (nr,nz) array ! fsum_tmp = 0. do m=m1,m2; do n=n1,n2 ; do i=l1,l2 fsum_tmp = fsum_tmp + fac*f(i,m,n,iax) enddo; enddo; enddo ! ! The sum has to be done processor-wise ! Sum over processors of same ipz, and different ipy ! call mpiallreduce_sum(fsum_tmp,mean_ax) ! ! Gauge-transform radial A ! f(l1:l2,m1:m2,n1:n2,iax) = f(l1:l2,m1:m2,n1:n2,iax) - mean_ax ! endsubroutine remove_volume_average !******************************************************************** ! !******************************************************************** !************ DO NOT DELETE THE FOLLOWING ************** !******************************************************************** !** This is an automatically generated include file that creates ** !** copies dummy routines from nospecial.f90 for any Special ** !** routines not implemented in this file ** !** ** include '../special_dummies.inc' !******************************************************************** endmodule Special