! $Id$ ! ! This modules implements viscous heating and diffusion terms ! here for shock viscosity ! nu_total = nu + nu_shock*dx^2*smooth(max5(-(div u)))) ! ! NOTE: this works and has been tested for periodic boundaries. ! With the current version, if your shock fronts reach a non-periodic ! boundary, unexpected things may happen, so you should monitor the ! behavior on the boundaries in this case. ! ! !** 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 :: lshock = .true. ! ! MVAR CONTRIBUTION 0 ! MAUX CONTRIBUTION 1 ! COMMUNICATED AUXILIARIES 1 ! ! PENCILS PROVIDED shock; gshock(3); shock_perp; gshock_perp(3) ! !*************************************************************** ! module Shock ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! include 'shock.h' ! logical :: lshock_first=.true.,lshock_max5=.false., lshock_max3_interp=.false. logical :: lwith_extreme_div=.false. logical :: lgauss_integral=.false. logical :: lgauss_integral_comm_uu=.false. logical :: lcommunicate_uu=.true. logical :: lforce_periodic_shockviscosity=.false. logical :: lshock_modulation_z=.false. real :: div_threshold=0., div_scaling=1. real :: z_shock=0., width_shock=.1 real, dimension (3,3,3) :: smooth_factor ! namelist /shock_run_pars/ & lshock_first, lshock_max5, div_threshold, div_scaling, & lgauss_integral, lcommunicate_uu, lgauss_integral_comm_uu, & lshock_max3_interp, lforce_periodic_shockviscosity, & lshock_modulation_z, z_shock, width_shock ! integer :: idiag_shockmax=0 ! DIAG_DOC: Max shock number ! interface shock_max3 module procedure shock_max3_farray module procedure shock_max3_pencil endinterface ! interface shock_max3_interp module procedure shock_max3_pencil_interp endinterface ! interface shock_smooth module procedure shock_smooth_farray module procedure shock_smooth_pencil endinterface ! interface shock_divu module procedure shock_divu_farray module procedure shock_divu_pencil endinterface ! interface shock_divu_perp module procedure shock_divu_perp_pencil endinterface ! contains !*********************************************************************** subroutine register_shock() ! ! 19-nov-02/tony: coded ! 24-jan-05/tony: modified from visc_shock.f90 ! use FArrayManager ! call farray_register_auxiliary('shock',ishock,communicated=.true.) ! ! Writing files for use with IDL ! if (naux+naux_com < maux+maux_com) aux_var(aux_count)=',shock $' if (naux+naux_com == maux+maux_com) aux_var(aux_count)=',shock' aux_count=aux_count+1 if (lroot) write(15,*) 'shock = fltarr(mx,my,mz)*one' ! endsubroutine register_shock !*********************************************************************** subroutine initialize_shock(f) ! ! 20-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! ! Initialize shock profile to zero ! f(:,:,:,ishock)=0.0 ! ! Calculate factors for polynomial smoothing ! smooth_factor=1. if (nxgrid/=1) then smooth_factor(1,:,:)=smooth_factor(1,:,:)*0.25 smooth_factor(2,:,:)=smooth_factor(2,:,:) *0.5 smooth_factor(3,:,:)=smooth_factor(3,:,:)*0.25 else smooth_factor(1,:,:)=0. smooth_factor(3,:,:)=0. endif ! if (nygrid/=1) then smooth_factor(:,1,:)=smooth_factor(:,1,:)*0.25 smooth_factor(:,2,:)=smooth_factor(:,2,:) *0.5 smooth_factor(:,3,:)=smooth_factor(:,3,:)*0.25 else smooth_factor(:,1,:)=0. smooth_factor(:,3,:)=0. endif ! if (nzgrid/=1) then smooth_factor(:,:,1)=smooth_factor(:,:,1)*0.25 smooth_factor(:,:,2)=smooth_factor(:,:,2) *0.5 smooth_factor(:,:,3)=smooth_factor(:,:,3)*0.25 else smooth_factor(:,:,1)=0. smooth_factor(:,:,3)=0. endif if (lroot) print*,"initialize_shock: prenormalised shock_factor sum=",sum(smooth_factor) smooth_factor=smooth_factor/sum(smooth_factor) ! if (div_threshold/=0.) lwith_extreme_div=.true. ! ! Die if periodic boundary condition for shock viscosity, but not for ! velocity field. It can lead to subtle numerical errors near non- ! periodic boundaries if the shock viscosity is assumed periodic. ! if (lrun .and. .not. lforce_periodic_shockviscosity) then if (bcx(ishock)=='p' .and. .not. all(bcx(iux:iuz)=='p' )) then if (lroot) then print*, 'initialize_shock: shock viscosity has bcx=''p'', but the velocity field is not' print*, ' periodic! (you must set a proper boundary condition for the' print*, ' shock viscosity)' print*, 'initialize_shock: bcx=', bcx print*, 'initialize_shock: to suppress this error,' print*, ' set lforce_periodic_shockviscosity=T in &shock_run_pars' endif call fatal_error('initialize_shock','') endif if (bcy(ishock)=='p' .and. .not. all(bcy(iux:iuz)=='p')) then if (lroot) then print*, 'initialize_shock: shock viscosity has bcy=''p'', but the velocity field is not' print*, ' periodic! (you must set a proper boundary condition for the' print*, ' shock viscosity)' print*, 'initialize_shock: bcy=', bcy print*, 'initialize_shock: to suppress this error,' print*, ' set lforce_periodic_shockviscosity=T in &shock_run_pars' endif call fatal_error('initialize_shock','') endif if (bcz(ishock)=='p' .and. .not. all(bcz(iux:iuz)=='p')) then if (lroot) then print*, 'initialize_shock: shock viscosity has bcz=''p'', but the velocity field is not' print*, ' periodic! (you must set a proper boundary condition for the' print*, ' shock viscosity)' print*, 'initialize_shock: bcz=', bcz print*, 'initialize_shock: to suppress this error,' print*, ' set lforce_periodic_shockviscosity=T in &shock_run_pars' endif call fatal_error('initialize_shock','') endif endif ! endsubroutine initialize_shock !*********************************************************************** subroutine read_shock_run_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=shock_run_pars, IOSTAT=iostat) ! endsubroutine read_shock_run_pars !*********************************************************************** subroutine write_shock_run_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=shock_run_pars) ! endsubroutine write_shock_run_pars !*********************************************************************** subroutine rprint_shock(lreset,lwrite) ! ! Writes ishock to index.pro file ! ! 24-nov-03/tony: adapted from rprint_ionization ! use Diagnostics, only: parse_name ! logical :: lreset logical, intent(in), optional :: lwrite integer :: iname ! ! reset everything in case of reset ! (this needs to be consistent with what is defined above!) ! if (lreset) then idiag_shockmax=0 lwith_extreme_div=.false. lgauss_integral=.false. lgauss_integral_comm_uu=.false. endif ! ! iname runs through all possible names that may be listed in print.in ! if (lroot.and.ip<14) print*,'rprint_shock: run through parse list' do iname=1,nname call parse_name(iname,cname(iname),cform(iname),& 'shockmax',idiag_shockmax) enddo ! ! check for those quantities for which we want video slices ! if (lwrite_slices) then where(cnamev=='shock') cformv='DEFINED' endif ! endsubroutine rprint_shock !*********************************************************************** subroutine get_slices_shock(f,slices) ! ! Write slices for animation of shock variable. ! ! 26-jul-06/tony: coded ! use Slices_methods, only: assign_slices_scal ! real, dimension (mx,my,mz,mfarray), intent(in) :: f type (slice_data) :: slices ! ! Loop over slices ! select case (trim(slices%name)) ! ! Shock profile ! case ('shock'); call assign_slices_scal(slices,f,ishock) ! endselect ! endsubroutine get_slices_shock !*********************************************************************** subroutine pencil_criteria_shock() ! ! All pencils that the Viscosity module depends on are specified here. ! ! 20-11-04/anders: coded ! if (idiag_shockmax/=0) then lpenc_diagnos(i_shock)=.true. endif ! endsubroutine pencil_criteria_shock !*********************************************************************** subroutine pencil_interdep_shock(lpencil_in) ! ! Interdependency among pencils from the Viscosity module is specified here. ! ! 20-11-04/anders: coded ! logical, dimension (npencils) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_shock !*********************************************************************** subroutine calc_pencils_shock(f,p) ! ! Calculate Viscosity pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 20-11-04/anders: coded ! use Diagnostics, only: max_mn_name use Sub, only: grad ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(in) :: f intent(inout) :: p ! shock if (lpencil(i_shock)) p%shock=f(l1:l2,m,n,ishock) ! gshock if (lpencil(i_gshock)) call grad(f,ishock,p%gshock) ! if (ldiagnos.and.idiag_shockmax/=0) call max_mn_name(p%shock,idiag_shockmax) ! endsubroutine calc_pencils_shock !*********************************************************************** subroutine calc_shock_profile_simple(f) ! ! Calculate divu based shock profile to be used in viscosity and ! diffusion type terms. ! ! 23-nov-02/tony: coded ! !for debug use Debug_IO ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: tmp ! ! Exit if were're not using the "simple" shock profile code or it's the wrong time. ! if (lgauss_integral.or.lgauss_integral_comm_uu.or. & lcommunicate_uu.or.(lshock_first.and.(.not.lfirst))) return ! ! calculate shock viscosity only when nu_shock /=0 ! ! calculate (-divu)_+ ! call shock_divu(f,f(:,:,:,ishock)) ! if (lwith_extreme_div) then where ((f(:,:,:,ishock) > 0.) .and. (f(:,:,:,ishock) < div_threshold)) & f(:,:,:,ishock)=0. where (f(:,:,:,ishock) > 0.) & f(:,:,:,ishock)=f(:,:,:,ishock)*div_scaling f(:,:,:,ishock)=abs(f(:,:,:,ishock)) else f(:,:,:,ishock)=max(0., -f(:,:,:,ishock)) endif ! ! Take the max over 5 neighboring points and smooth. ! Note that this means that we'd need 4 ghost zones, so ! we use one-sided formulae on processor boundaries. ! For that use lshock_max5=.true. ! ! To get the same result with and without MPI ! we take only the max over 3 points. ! if (lshock_max5) then call shock_max5(f(:,:,:,ishock),tmp) else call shock_max3(f(:,:,:,ishock),tmp) endif ! ! Save test data and scale to match the maximum expected result of smoothing ! call shock_smooth(tmp,f(:,:,:,ishock)) ! ! scale with min(dx,dy,dz)**2 ! f(:,:,:,ishock) = f(:,:,:,ishock) * dxmin**2 ! !debug only line:- ! if (ip==0) call output(trim(directory_dist)//'/shockvisc.dat',f(:,:,:,ishock),1) endsubroutine calc_shock_profile_simple !*********************************************************************** subroutine calc_shock_profile(f) ! ! Calculate divu based shock profile to be used in viscosity and ! diffusion type terms. ! ! 12-apr-05/tony: coded ! use Boundcond use Mpicomm use Sub use Magnetic, only: bb_unitvec_shock ! logical :: early_finalize real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: tmp real, dimension (mz) :: modulation_z real, dimension (mx) :: penc,penc_perp real, dimension (mx,3) :: bb_hat integer :: jj,kk ! if ((.not.lshock_first).or.lfirst) then if (lgauss_integral) then ! ! Calculate contributions to other processors ! ! call shock_calc_externalboundary(f) ! ! Initiate asyncronous communication of contributions to other processors ! if (nxgrid/=1) call bcshock_per_x(f) call initiate_isendrcv_bdry(f,ishock,ishock) ! ! Calculate all local shock profile contributions ! ! call shock_calc_internalboundary(f) ! call shock_calc_body(f) ! ! Finalize all shock profile communications ! call finalize_isendrcv_bdry(f,ishock,ishock) if (nygrid/=1) call bcshock_per_y(f) if (nzgrid/=1) call bcshock_per_z(f) !FIX ME !! CALCULATE shock in REAL boundary locations... and neighbours! !!! ! do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,ishock)=max(f(l1:l2,m,n,ishock),0.) enddo; enddo ! call scale_and_chop_internalboundary(f) !f(:,:,:,ishock) = tmp * dxmin**2 elseif (lgauss_integral_comm_uu) then ! ! need to finalize communication early either for test purposes ! early_finalize=test_nonblocking ! ! Initiate (non-blocking) communication of uu and do boundary conditions. ! Required order: ! 1. x-boundaries (x-ghost zones will be communicated) - done above ! 2. communication ! 3. y- and z-boundaries ! call initiate_isendrcv_bdry(f,iux,iuz) if (early_finalize) then call finalize_isendrcv_bdry(f,iux,iuz) call boundconds_y(f,iux,iuz) call boundconds_z(f,iux,iuz) endif ! ! do loop over y and z ! set indices and check whether communication must now be completed ! if test_nonblocking=.true., we communicate immediately as a test. ! do imn=1,ny*nz n=nn(imn) m=mm(imn) lfirstpoint=(imn==1) ! true for very first m-n loop llastpoint=(imn==(ny*nz)) ! true for very last m-n loop ! ! make sure all ghost points are set ! if (.not.early_finalize.and.necessary(imn)) then call finalize_isendrcv_bdry(f,iux,iuz) call boundconds_y(f,iux,iuz) call boundconds_z(f,iux,iuz) endif ! ! Calculate all local shock profile contributions ! ! call shock_calc_body(f) enddo ! ! Scale and chop the shock profile, as needed ! do n=n1,n2; do m=m1,m2 f(l1:l2,m,n,ishock)=max(f(l1:l2,m,n,ishock),0.) enddo; enddo ! elseif (lcommunicate_uu) then ! Communicate uu ghost zones call boundconds_x(f,iux,iuz) call initiate_isendrcv_bdry(f,iux,iuz) f(:,:,:,ishock)=impossible ! ! Divu over internal region ! do n=n1+1,n2-1; do m=m1+1,m2-1 call shock_divu(f,penc) f(:,m,n,ishock)=max(0.,-penc) if (ldivu_perp) then call bb_unitvec_shock(f,bb_hat) call shock_divu_perp(f,bb_hat,penc,penc_perp) f(:,m,n,ishock_perp)=max(0.,-penc_perp) endif enddo; enddo ! ! Max3 over internal region ! if (lshock_max3_interp) then do n=n1+2,n2-2; do m=m1+2,m2-2 call shock_max3_interp(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3_interp(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif enddo; enddo else do n=n1+2,n2-2; do m=m1+2,m2-2 call shock_max3(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif enddo; enddo endif ! ! Smooth over internal region ! do n=n1+3,n2-3; do m=m1+3,m2-3 call shock_smooth(tmp,penc) f(:,m,n,ishock)=penc if (ldivu_perp) then !tobi: this needs to be commented out !call shock_smooth(tmp_perp,penc_perp) f(:,m,n,ishock_perp)=penc_perp endif enddo; enddo ! ! End communication of uu ghost zones and set global boundary conditions. ! call finalize_isendrcv_bdry(f,iux,iuz) call boundconds_y(f,iux,iuz) call boundconds_z(f,iux,iuz) ! ! Divu over external region ! do n=2,mz-1; do jj=1,3 m=1+jj call shock_divu(f,penc) f(:,m,n,ishock)=max(-penc,0.) if (ldivu_perp) then call bb_unitvec_shock(f,bb_hat) call shock_divu_perp(f,bb_hat,penc,penc_perp) f(:,m,n,ishock_perp)=max(0.,-penc_perp) endif m=my-jj call shock_divu(f,penc) f(:,m,n,ishock)=max(-penc,0.) if (ldivu_perp) then call bb_unitvec_shock(f,bb_hat) call shock_divu_perp(f,bb_hat,penc,penc_perp) f(:,m,n,ishock_perp)=max(0.,-penc_perp) endif enddo; enddo do kk=1,3; do m=5,my-4 n=1+kk call shock_divu(f,penc) f(:,m,n,ishock)=max(-penc,0.) if (ldivu_perp) then call bb_unitvec_shock(f,bb_hat) call shock_divu_perp(f,bb_hat,penc,penc_perp) f(:,m,n,ishock_perp)=max(0.,-penc_perp) endif n=mz-kk call shock_divu(f,penc) f(:,m,n,ishock)=max(-penc,0.) if (ldivu_perp) then call bb_unitvec_shock(f,bb_hat) call shock_divu_perp(f,bb_hat,penc,penc_perp) f(:,m,n,ishock_perp)=max(0.,-penc_perp) endif enddo; enddo ! ! Max over external region ! if (lshock_max3_interp) then do n=3,mz-2; do jj=2,4 m=1+jj call shock_max3_interp(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3_interp(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif m=my-jj call shock_max3_interp(f,ishock,penc) tmp(:,m,n)=penc enddo; enddo do kk=2,4; do m=6,my-5 n=1+kk call shock_max3_interp(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3_interp(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif n=mz-kk call shock_max3_interp(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3_interp(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif enddo; enddo else do n=3,mz-2; do jj=2,4 m=1+jj call shock_max3(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif m=my-jj call shock_max3(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif enddo; enddo do kk=2,4; do m=6,my-5 n=1+kk call shock_max3(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif n=mz-kk call shock_max3(f,ishock,penc) tmp(:,m,n)=penc if (ldivu_perp) then call shock_max3(f,ishock_perp,penc_perp) !tobi: this needs to be commented out !tmp_perp(:,m,n)=penc_perp endif enddo; enddo endif ! ! Smooth over external region ! do n=4,mz-3; do jj=3,5 m=1+jj call shock_smooth(tmp,penc) f(:,m,n,ishock)=penc if (ldivu_perp) then !tobi: this needs to be commented out !call shock_smooth(tmp_perp,penc_perp) f(:,m,n,ishock_perp)=penc_perp endif m=my-jj call shock_smooth(tmp,penc) f(:,m,n,ishock)=penc if (ldivu_perp) then !tobi: this needs to be commented out !call shock_smooth(tmp_perp,penc_perp) f(:,m,n,ishock_perp)=penc_perp endif enddo; enddo do kk=3,5; do m=7,my-6 n=1+kk call shock_smooth(tmp,penc) f(:,m,n,ishock)=penc if (ldivu_perp) then !tobi: this needs to be commented out !call shock_smooth(tmp_perp,penc_perp) f(:,m,n,ishock_perp)=penc_perp endif n=mz-kk call shock_smooth(tmp,penc) f(:,m,n,ishock)=penc if (ldivu_perp) then !tobi: this needs to be commented out !call shock_smooth(tmp_perp,penc_perp) f(:,m,n,ishock_perp)=penc_perp endif enddo; enddo ! ! Non-MPI version: ! ! do n=2,mz-1; do m=2,my-1 ! call shock_divu(f,penc) ! f(:,m,n,ishock)=max(0.,-penc) ! enddo; enddo ! do n=3,mz-2; do m=3,my-2 ! call shock_max3(f,ishock,penc) ! tmp(:,m,n)=penc ! enddo; enddo ! do n=4,mz-3; do m=4,my-3 ! call shock_smooth(tmp,penc) ! f(:,m,n,ishock)=penc ! enddo; enddo !! !! Set boundary conditions on shock viscosity in the x-direction. The two other !! directions are taken care of later by pde (early_finalize does not work !! properly with shock). !! ! call boundconds_x(f,ishock,ishock) ! ! ! Scale with min(dx,dy,dz)**2 ! f(:,:,:,ishock) = f(:,:,:,ishock) * dxmin**2 if (ldivu_perp) then f(:,:,:,ishock_perp) = f(:,:,:,ishock_perp) * dxmin**2 endif ! ! allow for overall profile to turn off shock viscosity below z_shock with width_shock ! if (lshock_modulation_z) then modulation_z=.5*(1.+erfunc((z-z_shock)/width_shock)) f(:,:,:,ishock)=f(:,:,:,ishock)*spread(spread(modulation_z,1,mx),2,my) endif ! endif endif ! !debug only line:- ! if (ip==0) call output(trim(directory_dist)//'/shockvisc.dat',f(:,:,:,ishock),1) endsubroutine calc_shock_profile !*********************************************************************** !Utility routines - poss need moving elsewhere subroutine shock_max5(f,maxf) ! ! return array maxed with by 2 points either way ! skipping 1 data point all round ! ! 23-nov-02/tony: coded - from sub.f90 - nearmax ! 26-may-03/axel: maxf and f where interchanged in y-chunk ! real, dimension (mx,my,mz) :: f real, dimension (mx,my,mz) :: maxf, tmp ! ! x-direction, f -> maxf ! check for degeneracy ! if (nxgrid/=1) then if (mx>=5) then maxf(1 ,:,:) = max(f(1 ,:,:), & f(2 ,:,:), & f(3 ,:,:)) maxf(2 ,:,:) = max(f(1 ,:,:), & f(2 ,:,:), & f(3 ,:,:), & f(4 ,:,:)) maxf(3:mx-2,:,:) = max(f(1:mx-4,:,:), & f(2:mx-3,:,:), & f(3:mx-2,:,:), & f(4:mx-1,:,:), & f(5:mx ,:,:)) maxf(mx-1,:,:) = max(f(mx-3 ,:,:), & f(mx-2 ,:,:), & f(mx-1 ,:,:), & f(mx ,:,:)) maxf(mx ,:,:) = max(f(mx-2 ,:,:), & f(mx-1 ,:,:), & f(mx ,:,:)) elseif (mx==4) then maxf(1,:,:)=max(f(1,:,:),f(2,:,:),f(3,:,:)) maxf(2,:,:)=max(f(1,:,:),f(2,:,:),f(3,:,:),f(4,:,:)) maxf(3,:,:)=maxf(2,:,:) maxf(4,:,:)=max(f(2,:,:),f(3,:,:),f(4,:,:)) elseif (mx==3) then maxf(1,:,:)=max(f(1,:,:),f(2,:,:),f(3,:,:)) maxf(2,:,:)=maxf(1,:,:) maxf(3,:,:)=maxf(1,:,:) elseif (mx==2) then maxf(1,:,:)=max(f(1,:,:),f(2,:,:)) maxf(2,:,:)=maxf(1,:,:) else maxf=f endif else maxf=f endif ! ! y-direction, maxf -> f (swap back again) ! check for degeneracy ! if (nygrid/=1) then if (my>=5) then tmp(:,1 ,:) = max(maxf(:,1 ,:), & maxf(:,2 ,:), & maxf(:,3 ,:)) tmp(:,2 ,:) = max(maxf(:,1 ,:), & maxf(:,2 ,:), & maxf(:,3 ,:), & maxf(:,4 ,:)) tmp(:,3:my-2,:) = max(maxf(:,1:my-4,:), & maxf(:,2:my-3,:), & maxf(:,3:my-2,:), & maxf(:,4:my-1,:), & maxf(:,5:my ,:)) tmp(:,my-1,:) = max(maxf(:,my-3 ,:), & maxf(:,my-2 ,:), & maxf(:,my-1 ,:), & maxf(:,my ,:)) tmp(:,my ,:) = max(maxf(:,my-2 ,:), & maxf(:,my-1 ,:), & maxf(:,my ,:)) elseif (my==4) then tmp(:,1,:)=max(maxf(:,1,:),maxf(:,2,:),maxf(:,3,:)) tmp(:,2,:)=max(maxf(:,1,:),maxf(:,2,:),maxf(:,3,:),maxf(:,4,:)) tmp(:,3,:)=tmp(:,2,:) tmp(:,4,:)=max(maxf(:,2,:),maxf(:,3,:),maxf(:,4,:)) elseif (my==3) then tmp(:,1,:)=max(maxf(:,1,:),maxf(:,2,:),maxf(:,3,:)) tmp(:,2,:)=f(:,1,:) tmp(:,3,:)=f(:,1,:) elseif (my==2) then tmp(:,1,:)=max(maxf(:,1,:),maxf(:,2,:)) tmp(:,2,:)=tmp(:,1,:) else tmp=maxf endif else tmp=maxf endif ! ! z-direction, f -> maxf ! check for degeneracy ! if (nzgrid/=1) then if (mz>=5) then maxf(:,:,1 ) = max(tmp(:,:,1 ), & tmp(:,:,2 ), & tmp(:,:,3 )) maxf(:,:,2 ) = max(tmp(:,:,1 ), & tmp(:,:,2 ), & tmp(:,:,3 ), & tmp(:,:,4 )) maxf(:,:,3:mz-2) = max(tmp(:,:,1:mz-4), & tmp(:,:,2:mz-3), & tmp(:,:,3:mz-2), & tmp(:,:,4:mz-1), & tmp(:,:,5:mz )) maxf(:,:,mz-1 ) = max(tmp(:,:,mz-3 ), & tmp(:,:,mz-2 ), & tmp(:,:,mz-1 ), & tmp(:,:,mz )) maxf(:,:,mz ) = max(tmp(:,:,mz-2 ), & tmp(:,:,mz-1 ), & tmp(:,:,mz )) elseif (mz==4) then maxf(:,:,1)=max(tmp(:,:,1),tmp(:,:,2),tmp(:,:,3)) maxf(:,:,2)=max(tmp(:,:,1),tmp(:,:,2),tmp(:,:,3),tmp(:,:,4)) maxf(:,:,3)=maxf(:,:,2) maxf(:,:,4)=max(tmp(:,:,2),tmp(:,:,3),tmp(:,:,4)) elseif (mz==3) then maxf(:,:,1)=max(tmp(:,:,1),tmp(:,:,2),tmp(:,:,3)) maxf(:,:,2)=maxf(:,:,1) maxf(:,:,3)=maxf(:,:,1) elseif (mz==2) then maxf(:,:,1)=max(tmp(:,:,1),tmp(:,:,2)) maxf(:,:,2)=maxf(:,:,1) else maxf=tmp endif else maxf=tmp endif ! endsubroutine shock_max5 !*********************************************************************** !Utility routines - poss need moving elsewhere subroutine shock_max3_farray(f,maxf) ! ! return array maxed with by 2 points either way ! skipping 1 data point all round ! ! 29-nov-03/axel: adapted from shock_max5 ! real, dimension (mx,my,mz) :: f real, dimension (mx,my,mz) :: maxf, tmp ! ! x-direction, f -> maxf ! check for degeneracy ! if (nxgrid/=1) then if (mx>=3) then maxf(1 ,:,:) = max(f(1 ,:,:), & f(2 ,:,:)) maxf(2:mx-1,:,:) = max(f(1:mx-2,:,:), & f(2:mx-1,:,:), & f(3:mx ,:,:)) maxf( mx ,:,:) = max(f( mx-1,:,:), & f( mx ,:,:)) else maxf=f endif else maxf=f endif ! ! y-direction, maxf -> f (swap back again) ! check for degeneracy ! if (nygrid/=1) then if (my>=3) then tmp(:,1 ,:) = max(maxf(:,1 ,:), & maxf(:,2 ,:)) tmp(:,2:my-1,:) = max(maxf(:,1:my-2,:), & maxf(:,2:my-1,:), & maxf(:,3:my ,:)) tmp(:, my ,:) = max(maxf(:, my-1,:), & maxf(:, my ,:)) else tmp=maxf endif else tmp=maxf endif ! ! z-direction, f -> maxf ! check for degeneracy ! if (nzgrid/=1) then if (mz>=3) then maxf(:,:,1 ) = max(tmp(:,:,1 ), & tmp(:,:,2 )) maxf(:,:,2:mz-1) = max(tmp(:,:,1:mz-2), & tmp(:,:,2:mz-1), & tmp(:,:,3:mz )) maxf(:,:,mz ) = max(tmp(:,:, mz-1), & tmp(:,:, mz )) endif else maxf=tmp endif ! endsubroutine shock_max3_farray !*********************************************************************** !Utility routines - poss need moving elsewhere subroutine shock_max3_pencil(f,j,maxf) ! ! return array maxed with by 2 points either way ! skipping 1 data point all round ! ! 27-apr-03/tony: adapted from shock_max3 ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx) :: maxf integer :: j integer :: ii,jj,kk ! maxf=f(:,m,n,j) if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do jj=-1,1 do ii=-1,1 if (kk/=0.or.jj/=0.or.ii/=0) then maxf(l1-1:l2+1)=max(maxf(l1-1:l2+1),f(l1-1+ii:l2+1+ii,m+jj,n+kk,j)) endif enddo enddo enddo elseif ((nxgrid/=1).and.(nygrid/=1)) then do jj=-1,1 do ii=-1,1 if (jj/=0.or.ii/=0) then maxf(l1-1:l2+1)=max(maxf(l1-1:l2+1),f(l1-1+ii:l2+1+ii,m+jj,n1 ,j)) endif enddo enddo elseif ((nxgrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do ii=-1,1 if (kk/=0.or.ii/=0) then maxf(l1-1:l2+1)=max(maxf(l1-1:l2+1),f(l1-1+ii:l2+1+ii,m1 ,n+kk,j)) endif enddo enddo elseif ((nygrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do jj=-1,1 if (kk/=0.or.jj/=0) then maxf(l1 :l2 )=max(maxf(l1 :l2 ),f(l1 :l2 ,m+jj,n+kk,j)) endif enddo enddo elseif (nxgrid/=1) then do ii=-1,1 if (ii/=0) then maxf(l1-1:l2+1)=max(maxf(l1-1:l2+1),f(l1-1+ii:l2+1+ii,m1 ,n1 ,j)) endif enddo elseif (nygrid/=1) then do jj=-1,1 if (jj/=0) then maxf(l1 :l2 )=max(maxf(l1 :l2 ),f(l1 :l2 ,m+jj,n1 ,j)) endif enddo elseif (nzgrid/=1) then do kk=-1,1 if (kk/=0) then maxf(l1 :l2 )=max(maxf(l1 :l2 ),f(l1 :l2 ,m1 ,n+kk,j)) endif enddo endif ! ! endsubroutine shock_max3_pencil !*********************************************************************** !Utility routines - poss need moving elsewhere subroutine shock_max3_pencil_interp(f,j,maxf) ! ! return array maxed with by 2 points either way ! skipping 1 data point all round ! ! 14-aug-06/tony: adapted from shock_max3 ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx) :: tmp_penc real, dimension (mx) :: maxf real, parameter :: t1 = 1. real, parameter :: t2 = 0.70710678 real, parameter :: t3 = 0.57735027 real, parameter, dimension (-1:1,-1:1,-1:1) :: interp = reshape(& (/ t3, t2, t3, & t2, t1, t2, & t3, t2, t3, & t2, t1, t2, & t1, 0., t1, & t2, t1, t2, & t3, t2, t3, & t2, t1, t2, & t3, t2, t3 /),& (/ 3,3,3 /)) integer :: j integer :: ii,jj,kk ! maxf=0. if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then tmp_penc=f(:,m,n,j) do kk=-1,1 do jj=-1,1 do ii=-1,1 maxf(3:mx-2)=max(maxf(3:mx-2),interp(ii,jj,kk)*(f(3+ii:mx-2+ii,m+jj,n+kk,j)-tmp_penc(3:mx-2))) enddo enddo enddo else call fatal_error("shock_max3_pencil_interp", & "Tony got lazy and only implemented the 3D case") endif maxf(3:mx-2)=maxf(3:mx-2)+f(3:mx-2,m,n,j) ! ! endsubroutine shock_max3_pencil_interp !*********************************************************************** subroutine shock_max5_pencil(f,j,maxf) ! ! return array maxed with by 2 points either way ! skipping 1 data point all round ! ! 27-apr-03/tony: adapted from shock_max3 ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx) :: maxf integer :: j integer :: ii,jj,kk ! maxf=f(:,m,n,j) if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do kk=-2,2 if (kk==0) cycle do jj=-2,2 if (jj==0) cycle do ii=-2,2 if (ii==0) cycle maxf(2+ii:mx-1+ii)=max(maxf(2+ii:mx-1+ii),f(2+ii:mx-1+ii,m+jj,n+kk,j)) enddo enddo enddo elseif ((nxgrid/=1).and.(nygrid/=1)) then do jj=-2,2 if (jj==0) cycle do ii=-2,2 if (ii==0) cycle maxf(2+ii:mx-1+ii)=max(maxf(2+ii:mx-1+ii),f(2+ii:mx-1+ii,m+jj,n1,j)) enddo enddo elseif ((nxgrid/=1).and.(nzgrid/=1)) then do kk=-2,2 if (kk==0) cycle do ii=-2,2 if (ii==0) cycle maxf(2+ii:mx-1+ii)=max(maxf(2+ii:mx-1+ii),f(2+ii:mx-1+ii,m1,n+kk,j)) enddo enddo elseif ((nygrid/=1).and.(nzgrid/=1)) then do kk=-2,2 if (kk==0) cycle do jj=-2,2 if (jj==0) cycle maxf(l1:l2)=max(maxf(l1:l2),f(l1:l2,m+jj,n+kk,j)) enddo enddo elseif (nxgrid/=1) then do ii=-2,2 if (ii==0) cycle maxf(2+ii:mx-1+ii)=max(maxf(2+ii:mx-1+ii),f(2+ii:mx-1+ii,m1,n1,j)) enddo elseif (nygrid/=1) then do jj=-2,2 if (jj==0) cycle maxf(l1:l2)=max(maxf(l1:l2),f(l1:l2,m+jj,n1,j)) enddo elseif (nzgrid/=1) then do kk=-2,2 if (kk==0) cycle maxf(l1:l2)=max(maxf(l1:l2),f(l1:l2,m1,n+kk,j)) enddo endif ! endsubroutine shock_max5_pencil !*********************************************************************** subroutine shock_smooth_farray(f,smoothf) ! ! return array smoothed with by 2 points either way ! skipping 3 data points all round ! i.e. result valid () ! ! 23-nov-02/tony: coded ! real, dimension (mx,my,mz) :: f real, dimension (mx,my,mz) :: smoothf, tmp ! ! check for degeneracy ! if (nxgrid/=1) then if (mx>=3) then smoothf(1 ,:,:) = 0.25 * (3.*f(1 ,:,:) + & f(2 ,:,:)) ! smoothf(2:mx-1,:,:) = 0.25 * ( f(1:mx-2,:,:) + & 2.*f(2:mx-1,:,:) + & f(3:mx ,:,:)) ! smoothf(mx ,:,:) = 0.25 * ( f(mx-1 ,:,:) + & 3.*f(mx ,:,:)) else smoothf=f endif else smoothf=f endif ! ! check for degeneracy ! if (nygrid/=1) then if (my>=3) then tmp(:,1 ,:) = 0.25 * (3.*smoothf(:,1 ,:) + & smoothf(:,2 ,:)) ! tmp(:,2:my-1,:) = 0.25 * ( smoothf(:,1:my-2,:) + & 2.*smoothf(:,2:my-1,:) + & smoothf(:,3:my ,:)) ! tmp(:,my ,:) = 0.25 * ( smoothf(:,my-1 ,:) + & 3.*smoothf(:,my ,:)) else tmp=smoothf endif else tmp=smoothf endif ! ! check for degeneracy ! if (nzgrid/=1) then if (mz>=3) then smoothf(:,:,1 ) = 0.25 * (3.*tmp(:,:,1 ) + & tmp(:,:,2 )) ! smoothf(:,:,2:mz-1) = 0.25 * ( tmp(:,:,1:mz-2) + & 2.*tmp(:,:,2:mz-1) + & tmp(:,:,3:mz )) ! smoothf(:,:,mz ) = 0.25 * ( tmp(:,:,mz-1 ) + & 3.*tmp(:,:,mz )) else smoothf=tmp endif else smoothf=tmp endif ! endsubroutine shock_smooth_farray !*********************************************************************** subroutine shock_smooth_pencil(f,smoothf) ! ! return array smoothed with by 2 points either way ! skipping 3 data point all round ! i.e. result valid () ! ! 23-nov-02/tony: coded ! real, dimension (mx,my,mz) :: f real, dimension (mx) :: smoothf integer :: ii,jj,kk ! smoothf=0. if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do jj=-1,1 do ii=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2+ii,2+jj,2+kk)*f(l1+ii:l2+ii,m+jj,n+kk) enddo enddo enddo elseif ((nxgrid/=1).and.(nygrid/=1)) then do jj=-1,1 do ii=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2+ii,2+jj,2 )*f(l1+ii:l2+ii,m+jj,n1 ) enddo enddo elseif ((nxgrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do ii=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2+ii,2 ,2+kk)*f(l1+ii:l2+ii,m1 ,n+kk) enddo enddo elseif ((nygrid/=1).and.(nzgrid/=1)) then do kk=-1,1 do jj=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2 ,2+jj,2+kk)*f(l1 :l2 ,m+jj,n+kk) enddo enddo elseif (nxgrid/=1) then do ii=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2+ii,2 ,2 )*f(l1+ii:l2+ii,m1 ,n1 ) enddo elseif (nygrid/=1) then do jj=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2 ,2+jj,2 )*f(l1 :l2 ,m+jj,n1 ) enddo elseif (nzgrid/=1) then do kk=-1,1 smoothf(l1:l2) = smoothf(l1:l2) & + smooth_factor(2 ,2 ,2+kk)*f(l1 :l2 ,m1 ,n+kk) enddo endif ! endsubroutine shock_smooth_pencil !*********************************************************************** subroutine shock_divu_farray(f,df) ! ! calculate divergence of a vector U, get scalar ! accurate to 2nd order, explicit, centred and left and right biased ! ! 23-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray), intent(in) :: f real, dimension (mx,my,mz), intent(out) :: df real :: fac ! df=0. ! if (nxgrid/=1) then fac=1./(2.*dx) df(1 ,:,:) = df(1 ,:,:) & + ( 4.*f(2,:,:,iux) & - 3.*f(1,:,:,iux) & - f(3,:,:,iux))*fac df(2:mx-1,:,:) = df(2:mx-1,:,:) & + ( f(3:mx,:,:,iux)-f(1:mx-2,:,:,iux) ) * fac df(mx ,:,:) = df(mx ,:,:) & + ( 3.*f(mx ,:,:,iux) & - 4.*f(mx-1,:,:,iux) & + f(mx-2,:,:,iux))*fac endif ! if (nygrid/=1) then fac=1./(2.*dy) df(:,1 ,:) = df(:,1 ,:) & + ( 4.*f(:,2,:,iuy) & - 3.*f(:,1,:,iuy) & - f(:,3,:,iuy))*fac df(:,2:my-1,:) = df(:,2:my-1,:) & + (f(:,3:my,:,iuy)-f(:,1:my-2,:,iuy))*fac df(:,my ,:) = df(:,my ,:) & + ( 3.*f(:,my ,:,iuy) & - 4.*f(:,my-1,:,iuy) & + f(:,my-2,:,iuy))*fac endif ! if (nzgrid/=1) then fac=1./(2.*dz) df(:,:,1 ) = df(:,:,1 ) & + ( 4.*f(:,:,2,iuz) & - 3.*f(:,:,1,iuz) & - f(:,:,3,iuz))*fac df(:,:,2:mz-1) = df(:,:,2:mz-1) & + (f(:,:,3:mz,iuz)-f(:,:,1:mz-2,iuz))*fac df(:,:,mz ) = df(:,:,mz ) & + ( 3.*f(:,:,mz ,iuz) & - 4.*f(:,:,mz-1,iuz) & + f(:,:,mz-2,iuz))*fac endif ! endsubroutine shock_divu_farray !*********************************************************************** subroutine shock_divu_pencil(f,df) ! ! calculate divergence of a vector U, get scalar ! accurate to 2nd order, explicit, centred an left and right biased ! ! 23-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx) :: df real :: fac ! df=0. if (nxgrid/=1) then fac=1./(2.*dx) df(2:mx-1) = df(2:mx-1) & + (f(3:mx ,m ,n ,iux) & - f(1:mx-2,m ,n ,iux) ) & * fac endif ! if (nygrid/=1) then fac=1./(2.*dy) df = df & + (f(: ,m+1 ,n ,iuy) & - f(: ,m-1 ,n ,iuy) ) & * fac endif ! if (nzgrid/=1) then fac=1./(2.*dz) df = df & + (f(: ,m ,n+1 ,iuz) & - f(: ,m ,n-1 ,iuz) ) & * fac endif endsubroutine shock_divu_pencil !*********************************************************************** subroutine shock_divu_perp_pencil(f,bb_hat,divu,divu_perp) ! ! Calculate `perpendicular divergence' of u. ! nabla_perp.uu = nabla.uu - (1/b2)*bb.(bb.nabla)*uu ! ! 16-aug-06/tobi: coded ! real, dimension (mx,my,mz,mfarray), intent (in) :: f real, dimension (mx,3), intent (in) :: bb_hat real, dimension (mx), intent (in) :: divu real, dimension (mx), intent (out) :: divu_perp ! real, dimension (mx) :: fac ! divu_perp = divu ! if (nxgrid/=1) then fac=bb_hat(:,1)/(2*dx) divu_perp(l1-2:l2+2) = divu_perp(l1-2:l2+2) & - fac(l1-2:l2+2)*(bb_hat(l1-2:l2+2,1)*( f(l1-1:l2+3,m ,n ,iux) & - f(l1-3:l2+1,m ,n ,iux) ) & + bb_hat(l1-2:l2+2,2)*( f(l1-1:l2+3,m ,n ,iuy) & - f(l1-3:l2+1,m ,n ,iuy) ) & + bb_hat(l1-2:l2+2,3)*( f(l1-1:l2+3,m ,n ,iuz) & - f(l1-3:l2+1,m ,n ,iuz) ) ) endif ! if (nygrid/=1) then fac=bb_hat(:,2)/(2*dy) divu_perp(l1-2:l2+2) = divu_perp(l1-2:l2+2) & - fac(l1-2:l2+2)*(bb_hat(l1-2:l2+2,1)*( f(l1-2:l2+2,m+1,n ,iux) & - f(l1-2:l2+2,m-1,n ,iux) ) & + bb_hat(l1-2:l2+2,2)*( f(l1-2:l2+2,m+1,n ,iuy) & - f(l1-2:l2+2,m-1,n ,iuy) ) & + bb_hat(l1-2:l2+2,3)*( f(l1-2:l2+2,m+1,n ,iuz) & - f(l1-2:l2+2,m-1,n ,iuz) ) ) endif ! if (nzgrid/=1) then fac=bb_hat(:,3)/(2*dz) divu_perp(l1-2:l2+2) = divu_perp(l1-2:l2+2) & - fac(l1-2:l2+2)*(bb_hat(l1-2:l2+2,1)*( f(l1-2:l2+2,m ,n+1,iux) & - f(l1-2:l2+2,m ,n-1,iux) ) & + bb_hat(l1-2:l2+2,2)*( f(l1-2:l2+2,m ,n+1,iuy) & - f(l1-2:l2+2,m ,n-1,iuy) ) & + bb_hat(l1-2:l2+2,3)*( f(l1-2:l2+2,m ,n+1,iuz) & - f(l1-2:l2+2,m ,n-1,iuz) ) ) endif ! endsubroutine shock_divu_perp_pencil !*********************************************************************** subroutine shock_smooth_cube_diamond7(f,df) ! ! calculate divergence of a vector U smoothed over a 7x7x7 cube ! by using avergaging by a volume integral and transforming to ! a surface flux integral using Gauss' Theorm ! ! 01-apr-05/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: df real :: cube, diamond, fac integer :: i,j,k ! df=0. ! if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do k=n1,n2 do j=m1,m2 do i=l1,l2 ! df(i,j,k)=sum(f(i-3:i+3, j-3:j+3, k+3 , iuz)) - sum(f(i-3:i+3, j-3:j+3, k-3 , iuz)) + & ! sum(f(i-3:i+3, j+3 , k-3:k+3, iuy)) - sum(f(i-3:i+3, j-3 , k-3:k+3, iuy)) + & ! sum(f(i+3 , j-3:j+3, k-3:k+3, iux)) - sum(f(i-3 , j-3:j+3, k-3:k+3, iux)) enddo enddo enddo df=df*(dx*dy) elseif ((nxgrid/=1).and.(nzgrid/=1)) then do k=n1,n2 do j=m1,m2 do i=l1,l2 fac = (-1./18.) / dx diamond = ( & f(i , j, k+3 , iuz) & + f(i+1, j, k+2 , iuz) + f(i+1, j, k+2 , iux) & + f(i+2, j, k+1 , iuz) + f(i+2, j, k+1 , iux) & + f(i+3, j, k , iux) & - f(i+2, j, k-1 , iuz) + f(i+2, j, k-1 , iux) & - f(i+1, j, k-2 , iuz) + f(i+1, j, k-2 , iux) & - f(i , j, k-3 , iuz) & - f(i-1, j, k-2 , iuz) - f(i-1, j, k-2 , iux) & - f(i-2, j, k-1 , iuz) - f(i-2, j, k-1 , iux) & - f(i-3, j, k , iux) & + f(i-2, j, k+1 , iuz) - f(i-2, j, k+1 , iux) & + f(i-1, j, k+2 , iuz) - f(i-1, j, k+2 , iux) & ) * fac ! fac = (1./16.) / dx cube = ( sum(f(i-2:i+2, j, k-2 , iuz)) & - sum(f(i-2:i+2, j, k+2 , iuz)) & + sum(f(i-2 , j, k-2:k+2, iux)) & - sum(f(i+2 , j, k-2:k+2, iux)) & ) * fac ! if (lwith_extreme_div) then if ((diamond < 0.) .and. (diamond > -div_threshold)) then diamond=0. elseif (diamond<-div_threshold) then diamond=diamond*div_scaling endif if ((cube < 0.) .and. (cube > -div_threshold)) then cube=0. elseif (cube<-div_threshold) then cube=cube*div_scaling endif df(i,j,k)=0.5*(diamond + cube) !df(i,j,k)=sqrt(0.5*(diamond**2 + cube**2)) else df(i,j,k)=0.5*(max(diamond,0.) + max(cube,0.)) !df(i,j,k)=sqrt(0.5*(max(diamond,0.)**2 + max(cube,0.)**2)) endif enddo enddo enddo else call fatal_error('shock_smooth_cube_diamond7','CASE NOT IMPLEMENTED') endif ! endsubroutine shock_smooth_cube_diamond7 !*********************************************************************** function scale_and_chop(value) ! real :: scale_and_chop real :: value ! if (lwith_extreme_div) then if ((value < 0.) .and. (value > -div_threshold)) then scale_and_chop=0. elseif (value<-div_threshold) then scale_and_chop=-value*div_scaling else scale_and_chop=value endif else scale_and_chop=max(value,0.) endif ! endfunction scale_and_chop !*********************************************************************** subroutine scale_and_chop_internalboundary(f) ! real, dimension (mx,my,mz,mfarray) :: f integer :: i,j,k ! if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do j=m1i,m2i do i=0,nghost-1 f(l1+i,j,n1,ishock)=scale_and_chop(f(l1+i,j,n1,ishock)) f(l2-i,j,n1,ishock)=scale_and_chop(f(l2-i,j,n1,ishock)) enddo enddo do j=0,nghost-1 do i=l1,l2 f(i,m1+j,n1,ishock)=scale_and_chop(f(i,m1+j,n1,ishock)) f(i,m2-j,n1,ishock)=scale_and_chop(f(i,m2-j,n1,ishock)) enddo enddo elseif ((nxgrid/=1).and.(nygrid/=1)) then do j=m1i,m2i do i=0,nghost-1 f(l1+i,j,n1,ishock)=scale_and_chop(f(l1+i,j,n1,ishock)) f(l2-i,j,n1,ishock)=scale_and_chop(f(l2-i,j,n1,ishock)) enddo enddo do j=0,nghost-1 do i=l1,l2 f(i,m1+j,n1,ishock)=scale_and_chop(f(i,m1+j,n1,ishock)) f(i,m2-j,n1,ishock)=scale_and_chop(f(i,m2-j,n1,ishock)) enddo enddo elseif ((nxgrid/=1).and.(nzgrid/=1)) then do k=n1i,n2i do i=0,nghost-1 f(l1+i,m1,k,ishock)=scale_and_chop(f(l1+i,m1,k,ishock)) f(l2-i,m1,k,ishock)=scale_and_chop(f(l2-i,m1,k,ishock)) enddo enddo do k=0,nghost-1 do i=l1,l2 f(i,m1,n1+k,ishock)=scale_and_chop(f(i,m1,n1+k,ishock)) f(i,m1,n2-k,ishock)=scale_and_chop(f(i,m1,n2-k,ishock)) enddo enddo elseif ((nygrid/=1).and.(nzgrid/=1)) then do k=0,nghost-1 do j=m1i,m2i f(l1,j,n1+k,ishock)=scale_and_chop(f(l1,j,n1+k,ishock)) f(l1,j,n2-k,ishock)=scale_and_chop(f(l1,j,n2-k,ishock)) enddo enddo do k=n1,n2 do j=0,nghost-1 f(l1,m1+j,k,ishock)=scale_and_chop(f(l1,m1+j,k,ishock)) f(l1,m2-j,k,ishock)=scale_and_chop(f(l1,m2-j,k,ishock)) enddo enddo elseif (nxgrid/=1) then do i=0,nghost-1 f(l1+i,m1,n1,ishock)=scale_and_chop(f(l1+i,m1,n1,ishock)) f(l2-i,m1,n1,ishock)=scale_and_chop(f(l2-i,m1,n1,ishock)) enddo elseif (nygrid/=1) then do j=0,nghost-1 f(l1,m1+j,n1,ishock)=scale_and_chop(f(l1,m1+j,n1,ishock)) f(l1,m2-j,n1,ishock)=scale_and_chop(f(l1,m2-j,n1,ishock)) enddo elseif (nzgrid/=1) then do k=0,nghost-1 f(n1,m1,n1+k,ishock)=scale_and_chop(f(l1,m1,n1+k,ishock)) f(n1,m1,n2-k,ishock)=scale_and_chop(f(l1,m1,n2-k,ishock)) enddo endif ! endsubroutine scale_and_chop_internalboundary !*********************************************************************** subroutine shock_smooth_octagon7(f,df) ! ! calculate divergence of a vector U smoothed over a 7x7x7 cube ! by using avergaging by a volume integral and transforming to ! a surface flux integral using Gauss' Theorm ! ! 01-apr-05/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz) :: df real :: octagon, fac_diag, fac_straight integer :: i,j,k ! df=0. ! if ((nxgrid/=1).and.(nygrid/=1).and.(nzgrid/=1)) then do k=n1,n2 do j=m1,m2 do i=l1,l2 ! df(i,j,k)=sum(f(i-3:i+3, j-3:j+3, k+3 , iuz)) - sum(f(i-3:i+3, j-3:j+3, k-3 , iuz)) + & ! sum(f(i-3:i+3, j+3 , k-3:k+3, iuy)) - sum(f(i-3:i+3, j-3 , k-3:k+3, iuy)) + & ! sum(f(i+3 , j-3:j+3, k-3:k+3, iux)) - sum(f(i-3 , j-3:j+3, k-3:k+3, iux)) enddo enddo enddo df=df*(dx*dy) elseif ((nxgrid/=1).and.(nzgrid/=1)) then ! ! Top ! fac_diag = -sqrt(dx**2+dz**2)/(25.*dx*dz) fac_straight = -(1./(25*dx)) k=n2+1 j=m1 do i=l1,l2 octagon = ( & + f(i+2 , j, k+2 , iuz) + f(i+2 , j, k+2 , iux) & + f(i+3 , j, k+1 , iuz) + f(i+3 , j, k+1 , iux) & - f(i+1 , j, k-3 , iuz) + f(i+1 , j, k-3 , iux) & ! Bottom right / - f(i+2 , j, k-2 , iuz) + f(i+2 , j, k-2 , iux) & - f(i+3 , j, k-1 , iuz) + f(i+3 , j, k-1 , iux) & - f(i-1 , j, k-3 , iuz) - f(i-1 , j, k-3 , iux) & ! Bottom left \ - f(i-2 , j, k-2 , iuz) - f(i-2 , j, k-2 , iux) & - f(i-3 , j, k-1 , iuz) - f(i-3 , j, k-1 , iux) & + f(i-2 , j, k+2 , iuz) - f(i-2 , j, k+2 , iux) & + f(i-3 , j, k+1 , iuz) - f(i-3 , j, k+1 , iux) & ) * fac_diag & + ( & f(i+3 , j, k-1 , iux) - f(i-3 , j, k-1 , iux) & ! left and right | + f(i+3 , j, k , iux) - f(i-3 , j, k , iux) & + f(i+3 , j, k+1 , iux) - f(i-3 , j, k+1 , iux) & + f(i-1 , j, k+2 , iuz) - f(i-1 , j, k-3 , iuz) & ! top and bottom - + f(i , j, k+2 , iuz) - f(i , j, k-3 , iuz) & + f(i+1 , j, k+2 , iuz) - f(i+1 , j, k-3 , iuz) & + f(i-2 , j, k+2 , iuz) & + f(i+2 , j, k+2 , iuz) & ) * fac_straight ! df(i,j,k)=scale_and_chop(octagon) enddo ! ! Bottom ! fac_diag = -sqrt(dx**2+dz**2)/(25.*dx*dz) fac_straight = -(1./(25*dx)) k=n2+1 j=m1 do i=l1,l2 octagon = ( & + f(i+1 , j, k+3 , iuz) + f(i+1 , j, k+3 , iux) & + f(i+2 , j, k+2 , iuz) + f(i+2 , j, k+2 , iux) & + f(i+3 , j, k+1 , iuz) + f(i+3 , j, k+1 , iux) & - f(i+2 , j, k-2 , iuz) + f(i+2 , j, k-2 , iux) & - f(i+3 , j, k-1 , iuz) + f(i+3 , j, k-1 , iux) & - f(i-2 , j, k-2 , iuz) - f(i-2 , j, k-2 , iux) & - f(i-3 , j, k-1 , iuz) - f(i-3 , j, k-1 , iux) & + f(i-2 , j, k+2 , iuz) - f(i-2 , j, k+2 , iux) & + f(i-3 , j, k+1 , iuz) - f(i-3 , j, k+1 , iux) & ) * fac_diag & + ( & f(i+3 , j, k-1 , iux) - f(i-3 , j, k-1 , iux) & ! left and right | + f(i+3 , j, k , iux) - f(i-3 , j, k , iux) & + f(i+3 , j, k+1 , iux) - f(i-3 , j, k+1 , iux) & + f(i-1 , j, k+3 , iuz) - f(i-1 , j, k-2 , iuz) & ! top and bottom - + f(i , j, k+3 , iuz) - f(i , j, k-2 , iuz) & + f(i+1 , j, k+3 , iuz) - f(i+1 , j, k-2 , iuz) & - f(i-2 , j, k-2 , iuz) & - f(i+2 , j, k-2 , iuz) & ) * fac_straight ! df(i,j,k)=scale_and_chop(octagon) enddo ! ! Bulk ! fac_diag = -sqrt(dx**2+dz**2)/(28.*dx*dz) fac_straight = -(1./(28.*dx)) j=m1 do k=n1i,n2i do i=l1i,l2i octagon = ( & f(i+1 , j, k+3 , iuz) + f(i+1 , j, k+3 , iux) & ! Top right \ + f(i+2 , j, k+2 , iuz) + f(i+2 , j, k+2 , iux) & + f(i+3 , j, k+1 , iuz) + f(i+3 , j, k+1 , iux) & - f(i+1 , j, k-3 , iuz) + f(i+1 , j, k-3 , iux) & ! Bottom right / - f(i+2 , j, k-2 , iuz) + f(i+2 , j, k-2 , iux) & - f(i+3 , j, k-1 , iuz) + f(i+3 , j, k-1 , iux) & - f(i-1 , j, k-3 , iuz) - f(i-1 , j, k-3 , iux) & ! Bottom left \ - f(i-2 , j, k-2 , iuz) - f(i-2 , j, k-2 , iux) & - f(i-3 , j, k-1 , iuz) - f(i-3 , j, k-1 , iux) & + f(i-1 , j, k+3 , iuz) - f(i-1 , j, k+3 , iux) & ! Top left / + f(i-2 , j, k+2 , iuz) - f(i-2 , j, k+2 , iux) & + f(i-3 , j, k+1 , iuz) - f(i-3 , j, k+1 , iux) & ) * fac_diag & + ( & f(i+3 , j, k-1 , iux) - f(i-3 , j, k-1 , iux) & ! left and right | + f(i+3 , j, k , iux) - f(i-3 , j, k , iux) & + f(i+3 , j, k+1 , iux) - f(i-3 , j, k+1 , iux) & + f(i-1 , j, k+3 , iuz) - f(i-1 , j, k-3 , iuz) & ! top and bottom - + f(i , j, k+3 , iuz) - f(i , j, k-3 , iuz) & + f(i+1 , j, k+3 , iuz) - f(i+1 , j, k-3 , iuz) & ) * fac_straight ! df(i,j,k)=scale_and_chop(octagon) enddo enddo else call fatal_error('shock_smooth_octagon7','CASE NOT IMPLEMENTED') endif ! endsubroutine shock_smooth_octagon7 !*********************************************************************** subroutine bcshock_per_x(f) ! ! periodic boundary condition ! 11-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f ! if (nprocx==1.and.lperi(1)) then f(l1:l1i ,:,:,ishock) = f(l1:l1i,:,:,ishock) + f(l2+1:mx,:,:,ishock) ! f(l2+1:mx,:,:,ishock) = f(l1:l1i,:,:,ishock) f(l2i:l2 ,:,:,ishock) = f(l2i:l2,:,:,ishock) + f(1:l1-1 ,:,:,ishock) ! f(1:l1-1 ,:,:,ishock) = f(l2i:l2,:,:,ishock) endif ! endsubroutine bcshock_per_x !*********************************************************************** subroutine bcshock_per_y(f) ! ! periodic boundary condition ! 11-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f ! if (nprocy==1.and.lperi(2)) then f(:,m1:m1i,:,ishock) = f(:,m1:m1i,:,ishock) + f(:,m2+1:my,:,ishock) ! f(:,m2+1:my,:,ishock) = f(:,m1:m1i,:,ishock) f(:,m2i:m2,:,ishock) = f(:,m2i:m2,:,ishock) + f(:,1:m1-1 ,:,ishock) ! f(:,1:m1-1,:,ishock) = f(:,m2i:m2,:,ishock) endif ! endsubroutine bcshock_per_y !*********************************************************************** subroutine bcshock_per_z(f) ! ! periodic boundary condition ! 11-nov-02/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f ! if (nprocz==1.and.lperi(3)) then f(:,:,n1:n1i ,ishock) = f(:,:,n1:n1i,ishock) + f(:,:,n2+1:mz,ishock) ! f(:,:,n2+1:mz,ishock) = f(:,:,n1:n1i,ishock) f(:,:,n2i:n2 ,ishock) = f(:,:,n2i:n2,ishock) + f(:,:,1:n1-1 ,ishock) ! f(:,:,1:n1-1 ,ishock) = f(:,:,n2i:n2,ishock) endif ! endsubroutine bcshock_per_z !*********************************************************************** ! ! include 'shock_profile.inc' ! endmodule Shock