! $Id$ ! ! This module provide a way for users to specify custom ! (i.e. not in the standard Pencil Code) physics, diagnostics etc. ! ! The module provides a set of standard hooks into the Pencil-Code and ! currently allows the following customizations: ! ! Description | Relevant function call ! --------------------------------------------------------------------------- ! Special variable registration | register_special ! (pre parameter read) | ! Special variable initialization | initialize_special ! (post parameter read) | ! Special variable finalization | finalize_special ! (deallocation, etc.) | ! | ! Special initial condition | init_special ! this is called last so may be used to modify | ! the mvar variables declared by this module | ! or optionally modify any of the other f array | ! variables. The latter, however, should be | ! avoided where ever possible. | ! | ! Special term in the mass (density) equation | special_calc_density ! Special term in the momentum (hydro) equation | special_calc_hydro ! Special term in the energy equation | special_calc_energy ! Special term in the induction (magnetic) | special_calc_magnetic ! equation | ! | ! Special equation | dspecial_dt ! NOT IMPLEMENTED FULLY YET - HOOKS NOT PLACED INTO THE PENCIL-CODE ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of special_dummies.inc) the number of f array ! variables and auxiliary variables added by this module ! ! CPARAM logical, parameter :: lspecial = .true. ! ! MAUX CONTRIBUTION 18 ! ! PENCILS PROVIDED stress_ij(6) !*************************************************************** ! ! HOW TO USE THIS FILE ! -------------------- ! ! Change the line above to ! lspecial = .true. ! to enable use of special hooks. ! ! 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 SVN 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/geo_kws ! ! Where geo_kws it replaced by the filename of your new module ! upto and not including the .f90 ! module Special ! use Cparam use Cdata use Initcond use General, only: keep_compiler_quiet use Messages, only: svn_id, fatal_error ! implicit none ! include '../special.h' ! ! Declare index of new variables in f array (if any). ! character (len=labellen) :: inithij='nothing' character (len=labellen) :: initgij='nothing' character (len=labellen) :: ctrace_factor='1/3' character (len=labellen) :: cstress_prefactor='6' character (len=labellen) :: fourthird_in_stress='4/3' character (len=labellen) :: cc_light='1' character (len=labellen) :: aux_stress='stress' real :: amplhij=0., amplgij=0., dummy=0. real :: trace_factor=0., stress_prefactor, fourthird_factor, EGWpref real :: nscale_factor_conformal=1., tshift=0. logical :: lno_transverse_part=.false. logical :: lswitch_sign_e_X=.true., ldebug_print=.false., lkinGW=.true. logical :: lStress_as_aux=.true., lreynolds=.false. logical :: lggTX_as_aux=.true., lhhTX_as_aux=.true. logical :: lremove_mean_hij=.false., lremove_mean_gij=.false. real, dimension(3,3) :: ij_table real :: c_light2=1. ! ! Do this here because shared variables for this array doesn't work on Beskow. ! integer, parameter :: nk=nxgrid/2 real, dimension(nk) :: specGWs ,specGWh ,specStr real, dimension(nk) :: specGWshel,specGWhhel,specStrhel public :: specGWs, specGWshel, specGWh, specGWhhel, specStr, specStrhel ! ! input parameters namelist /special_init_pars/ & ctrace_factor, cstress_prefactor, fourthird_in_stress, lno_transverse_part, & inithij, initgij, amplhij, amplgij, lStress_as_aux, & lggTX_as_aux, lhhTX_as_aux ! ! run parameters namelist /special_run_pars/ & ctrace_factor, cstress_prefactor, fourthird_in_stress, lno_transverse_part, & ldebug_print, lswitch_sign_e_X, & nscale_factor_conformal, tshift, cc_light, & lStress_as_aux, lkinGW, aux_stress, & lggTX_as_aux, lhhTX_as_aux, lremove_mean_hij, lremove_mean_gij ! ! Diagnostic variables (needs to be consistent with reset list below). ! integer :: idiag_g11pt=0 ! DIAG_DOC: $g_{11}(x_1,y_1,z_1,t)$ integer :: idiag_g22pt=0 ! DIAG_DOC: $g_{22}(x_1,y_1,z_1,t)$ integer :: idiag_g33pt=0 ! DIAG_DOC: $g_{33}(x_1,y_1,z_1,t)$ integer :: idiag_g12pt=0 ! DIAG_DOC: $g_{12}(x_1,y_1,z_1,t)$ integer :: idiag_g23pt=0 ! DIAG_DOC: $g_{23}(x_1,y_1,z_1,t)$ integer :: idiag_g31pt=0 ! DIAG_DOC: $g_{31}(x_1,y_1,z_1,t)$ integer :: idiag_hhTpt=0 ! DIAG_DOC: $h_{T}(x_1,y_1,z_1,t)$ integer :: idiag_hhXpt=0 ! DIAG_DOC: $h_{X}(x_1,y_1,z_1,t)$ integer :: idiag_ggTpt=0 ! DIAG_DOC: $\dot{h}_{T}(x_1,y_1,z_1,t)$ integer :: idiag_ggXpt=0 ! DIAG_DOC: $\dot{h}_{X}(x_1,y_1,z_1,t)$ integer :: idiag_hhTp2=0 ! DIAG_DOC: $h_{T}(x_1,y_1,z_1,t)$ integer :: idiag_hhXp2=0 ! DIAG_DOC: $h_{X}(x_1,y_1,z_1,t)$ integer :: idiag_ggTp2=0 ! DIAG_DOC: $\dot{h}_{T}(x_1,y_1,z_1,t)$ integer :: idiag_ggXp2=0 ! DIAG_DOC: $\dot{h}_{X}(x_1,y_1,z_1,t)$ integer :: idiag_hrms=0 ! DIAG_DOC: $\bra{h_T^2+h_X^2}^{1/2}$ integer :: idiag_EEGW=0 ! DIAG_DOC: $\bra{g_T^2+g_X^2}\,c^2/(32\pi G)$ integer :: idiag_gg2m=0 ! DIAG_DOC: $\bra{g_T^2+g_X^2}$ integer :: idiag_hhT2m=0 ! DIAG_DOC: $\bra{h_T^2}$ integer :: idiag_hhX2m=0 ! DIAG_DOC: $\bra{h_X^2}$ integer :: idiag_hhTXm=0 ! DIAG_DOC: $\bra{h_T h_X}$ integer :: idiag_ggT2m=0 ! DIAG_DOC: $\bra{g_T^2}$ integer :: idiag_ggX2m=0 ! DIAG_DOC: $\bra{g_X^2}$ integer :: idiag_ggTXm=0 ! DIAG_DOC: $\bra{g_T g_X}$ ! contains !*********************************************************************** subroutine register_special() ! ! Set up indices for variables in special modules. ! ! 3-aug-17/axel: coded ! use Sub, only: register_report_aux use FArrayManager ! if (lroot) call svn_id( & "$Id$") ! ! Register ggT and ggX as auxiliary arrays ! May want to do this only when Fourier transform is enabled. ! if (lggTX_as_aux) then call farray_register_auxiliary('ggT',iggT) call farray_register_auxiliary('ggX',iggX) call farray_register_auxiliary('ggTim',iggTim) call farray_register_auxiliary('ggXim',iggXim) endif ! if (lhhTX_as_aux) then call farray_register_auxiliary('hhT',ihhT) call farray_register_auxiliary('hhX',ihhX) call farray_register_auxiliary('hhTim',ihhTim) call farray_register_auxiliary('hhXim',ihhXim) endif ! if (lStress_as_aux) then call farray_register_auxiliary('StT',iStressT) call farray_register_auxiliary('StX',iStressX) call farray_register_auxiliary('StTim',iStressTim) call farray_register_auxiliary('StXim',iStressXim) call farray_register_auxiliary('Str',iStress_ij,vector=6) endif ! if (lroot) call svn_id( & "$Id$") ! endsubroutine register_special !*********************************************************************** subroutine initialize_special(f) ! ! Called after reading parameters, but before the time loop. ! ! 06-oct-03/tony: coded ! use EquationOfState, only: cs0 ! real, dimension (mx,my,mz,mfarray) :: f ! ! set index table ! ij_table(1,1)=1 ij_table(2,2)=2 ij_table(3,3)=3 ij_table(1,2)=4 ij_table(2,3)=5 ij_table(3,1)=6 ij_table(2,1)=4 ij_table(3,2)=5 ij_table(1,3)=6 ! ! determine trace factor ! select case (ctrace_factor) case ('0'); trace_factor=.0 case ('1/2'); trace_factor=.5 case ('1/3'); trace_factor=onethird case default call fatal_error("initialize_special: No such value for ctrace_factor:" & ,trim(ctrace_factor)) endselect ! ! Determine fourthird_in_stress. This factor is normally 4/3, but it can be ! set to unity in case we want to pretend that the kinematic Beltrami field ! has the same prefactor as a magnetic one. ! select case (fourthird_in_stress) case ('1'); fourthird_factor=1. case ('4/3'); fourthird_factor=fourthird case default call fatal_error("initialize_special: No such value for fourthird_in_stress:" & ,trim(fourthird_in_stress)) endselect ! ! determine stress_prefactor and GW energy prefactor, ! which is EGWpref=.5*16.*pi/stress_prefactor**2 ! At the moment, only the case stress_prefactor=6 is to be used. ! The other cases are kept for backward compatibility. ! select case (cstress_prefactor) case ('1'); stress_prefactor=1.; EGWpref=8.*pi case ('6'); stress_prefactor=6.; EGWpref=1./(32.*pi) case ('16pi'); stress_prefactor=16.*pi; EGWpref=1./(32.*pi) case ('16piG/c^2'); stress_prefactor=16.*pi*G_Newton_cgs/c_light_cgs**2; EGWpref=c_light_cgs**2/(32.*pi*G_Newton_cgs) case default call fatal_error("initialize_special: No such value for ctrace_factor:" & ,trim(ctrace_factor)) endselect if (headt) print*,'stress_prefactor=',stress_prefactor if (headt) print*,'EGWpref=',EGWpref ! ! set speed of light ! select case (cc_light) case ('1'); c_light2=1. case ('cgs'); c_light2=c_light_cgs**2 case default call fatal_error("initialize_special: No such value for cc_light:" & ,trim(ctrace_factor)) endselect if (headt) print*,'c_light2=',c_light2 ! ! Determine whether Reynolds stress needs to be computed: ! Compute Reynolds stress when lhydro or lhydro_kinematic, ! provided lkinGW=T ! lreynolds=(lhydro.or.lhydro_kinematic).and.lkinGW ! ! give a warning if cs0**2=1 ! ! if (cs0==1.) call fatal_error('gravitational_waves_hij6', & ! 'cs0 should probably not be unity') ! call keep_compiler_quiet(f) ! endsubroutine initialize_special !*********************************************************************** subroutine finalize_special(f) ! ! Called right before exiting. ! ! 14-aug-2011/Bourdin.KIS: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! call keep_compiler_quiet(f) ! endsubroutine finalize_special !*********************************************************************** subroutine init_special(f) ! ! initialise special condition; called from start.f90 ! 06-oct-2003/tony: coded ! real, dimension (mx,my,mz,mfarray) :: f ! intent(inout) :: f ! ! initial condition for hij ! select case (inithij) case ('nothing') if (lroot) print*,'init_special: nothing' f(:,:,:,ihhT:ihhXim)=0. case default call fatal_error("init_special: No such value for inithij:" & ,trim(inithij)) endselect ! ! initial condition for gij ! select case (initgij) case ('nothing') if (lroot) print*,'init_special: nothing' f(:,:,:,iggT:iggXim)=0. case default call fatal_error("init_special: No such value for initgij:" & ,trim(initgij)) endselect ! endsubroutine init_special !*********************************************************************** subroutine pencil_criteria_special() ! ! All pencils that this special module depends on are specified here. ! ! 1-apr-06/axel: coded ! ! Velocity field needed for Reynolds stress ! if (lreynolds) then lpenc_requested(i_uu)=.true. lpenc_requested(i_rho)=.true. if (trace_factor/=0.) lpenc_requested(i_u2)=.true. endif ! ! Magnetic field needed for Maxwell stress ! if (lmagnetic) then lpenc_requested(i_bb)=.true. if (trace_factor/=0.) lpenc_requested(i_b2)=.true. endif ! endsubroutine pencil_criteria_special !*********************************************************************** subroutine pencil_interdep_special(lpencil_in) ! ! Interdependency among pencils provided by this module are specified here. ! ! 18-jul-06/tony: coded ! logical, dimension(npencils), intent(inout) :: lpencil_in ! call keep_compiler_quiet(lpencil_in) ! endsubroutine pencil_interdep_special !*********************************************************************** subroutine calc_pencils_special(f,p) ! ! Calculate Special pencils. ! Most basic pencils should come first, as others may depend on them. ! ! 24-aug-17/axel: coded ! 7-jun-18/axel: included 4/3*rho factor ! use Deriv, only: derij ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (nx) :: tmp type (pencil_case) :: p ! intent(in) :: f intent(inout) :: p integer :: i, j, ij ! ! Construct stress tensor; notice opposite signs for u and b. ! p%stress_ij=0.0 do j=1,3 do i=1,j ij=ij_table(i,j) if (lreynolds) p%stress_ij(:,ij)=p%stress_ij(:,ij)+p%uu(:,i)*p%uu(:,j)*fourthird_factor*p%rho if (lmagnetic) p%stress_ij(:,ij)=p%stress_ij(:,ij)-p%bb(:,i)*p%bb(:,j) ! ! Remove trace. ! if (i==j) then if (lreynolds) p%stress_ij(:,ij)=p%stress_ij(:,ij)-trace_factor*p%u2*fourthird_factor*p%rho if (lmagnetic) p%stress_ij(:,ij)=p%stress_ij(:,ij)+trace_factor*p%b2 endif enddo enddo ! 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 for ! efficiency. ! ! This routine computes various diagnostic quantities, but those ! would be more easily done further below when sign_switch is known. ! But this only affects the helicities. The values of EEGW and hrms are ! however correct and agree with those of gravitational_waves_hij6. ! ! 06-oct-03/tony: coded ! 07-feb-18/axel: added nscale_factor=0 (no expansion), =.5 (radiation era) ! use Diagnostics ! real, dimension (mx,my,mz,mfarray) :: f real, dimension (mx,my,mz,mvar) :: df real :: scale_factor, stress_prefactor2, sign_switch=0 type (pencil_case) :: p ! integer :: ij ! intent(in) :: p intent(inout) :: f,df ! ! Identify module and boundary conditions. ! if (headtt.or.ldebug) print*,'dspecial_dt: SOLVE dspecial_dt' ! ! Compute scale factor. ! Note: to prevent division by zero, it is best to put tstart=1. in start.in. ! If that is not done, one can put here tshift=1., for example. ! If that is not the case either, we put scale_factor=1. ! At the next timestep, this will no longer be a problem. ! if (t+tshift==0.) then scale_factor=1. else scale_factor=(t+tshift)**nscale_factor_conformal endif stress_prefactor2=stress_prefactor/scale_factor ! ! Assemble rhs of GW equations. ! do ij=1,6 f(l1:l2,m,n,iStress_ij+ij-1)=stress_prefactor2*p%stress_ij(:,ij) enddo ! ! diagnostics ! if (ldiagnos) then if (lggTX_as_aux) then if (idiag_EEGW/=0) call sum_mn_name((f(l1:l2,m,n,iggT)**2+f(l1:l2,m,n,iggTim)**2 & +f(l1:l2,m,n,iggX)**2+f(l1:l2,m,n,iggXim)**2 & )*nwgrid*EGWpref,idiag_EEGW) if (idiag_gg2m/=0) call sum_mn_name((f(l1:l2,m,n,iggT)**2+f(l1:l2,m,n,iggTim)**2 & +f(l1:l2,m,n,iggX)**2+f(l1:l2,m,n,iggXim)**2 & )*nwgrid,idiag_gg2m) if (idiag_ggT2m/=0) call sum_mn_name((f(l1:l2,m,n,iggT)**2+f(l1:l2,m,n,iggTim)**2 & )*nwgrid,idiag_ggT2m) if (idiag_ggX2m/=0) call sum_mn_name((f(l1:l2,m,n,iggX)**2+f(l1:l2,m,n,iggXim)**2 & )*nwgrid,idiag_ggX2m) if (idiag_ggTXm/=0) call sum_mn_name((f(l1:l2,m,n,iggT )*f(l1:l2,m,n,iggXim) & -f(l1:l2,m,n,iggTim)*f(l1:l2,m,n,iggX ) & )*nwgrid*sign_switch,idiag_ggTXm) endif if (lhhTX_as_aux) then if (idiag_hrms/=0) call sum_mn_name((f(l1:l2,m,n,ihhT)**2+f(l1:l2,m,n,ihhTim)**2 & +f(l1:l2,m,n,ihhX)**2+f(l1:l2,m,n,ihhXim)**2 & )*nwgrid,idiag_hrms,lsqrt=.true.) if (idiag_hhT2m/=0) call sum_mn_name((f(l1:l2,m,n,ihhT)**2+f(l1:l2,m,n,ihhTim)**2 & )*nwgrid,idiag_hhT2m) if (idiag_hhX2m/=0) call sum_mn_name((f(l1:l2,m,n,ihhX)**2+f(l1:l2,m,n,ihhXim)**2 & )*nwgrid,idiag_hhX2m) if (idiag_hhTXm/=0) call sum_mn_name((f(l1:l2,m,n,ihhT )*f(l1:l2,m,n,ihhXim) & -f(l1:l2,m,n,ihhTim)*f(l1:l2,m,n,ihhX ) & )*nwgrid*sign_switch,idiag_hhTXm) endif ! if (lroot.and.m==mpoint.and.n==npoint) then if (idiag_g11pt/=0) call save_name(f(lpoint,m,n,igij+1-1),idiag_g11pt) if (idiag_g22pt/=0) call save_name(f(lpoint,m,n,igij+2-1),idiag_g22pt) if (idiag_g33pt/=0) call save_name(f(lpoint,m,n,igij+3-1),idiag_g33pt) if (idiag_g12pt/=0) call save_name(f(lpoint,m,n,igij+4-1),idiag_g12pt) if (idiag_g23pt/=0) call save_name(f(lpoint,m,n,igij+5-1),idiag_g23pt) if (idiag_g31pt/=0) call save_name(f(lpoint,m,n,igij+6-1),idiag_g31pt) if (lhhTX_as_aux) then if (idiag_hhTpt/=0) call save_name(f(lpoint,m,n,ihhT),idiag_hhTpt) if (idiag_hhXpt/=0) call save_name(f(lpoint,m,n,ihhX),idiag_hhXpt) endif if (lggTX_as_aux) then if (idiag_ggTpt/=0) call save_name(f(lpoint,m,n,iggT),idiag_ggTpt) if (idiag_ggXpt/=0) call save_name(f(lpoint,m,n,iggX),idiag_ggXpt) endif endif ! if (lroot.and.m==mpoint2.and.n==npoint2) then if (lhhTX_as_aux) then if (idiag_hhTp2/=0) call save_name(f(lpoint2,m,n,ihhT),idiag_hhTp2) if (idiag_hhXp2/=0) call save_name(f(lpoint2,m,n,ihhX),idiag_hhXp2) endif if (lggTX_as_aux) then if (idiag_ggTp2/=0) call save_name(f(lpoint2,m,n,iggT),idiag_ggTp2) if (idiag_ggXp2/=0) call save_name(f(lpoint2,m,n,iggX),idiag_ggXp2) endif endif endif ! endsubroutine dspecial_dt !*********************************************************************** subroutine read_special_init_pars(iostat) ! use File_io, only: parallel_unit ! integer, intent(out) :: iostat ! read(parallel_unit, NML=special_init_pars, IOSTAT=iostat) ! endsubroutine read_special_init_pars !*********************************************************************** subroutine write_special_init_pars(unit) ! integer, intent(in) :: unit ! write(unit, NML=special_init_pars) ! endsubroutine write_special_init_pars !*********************************************************************** 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 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 ! ! 13-may-18/axel: added remove_mean_value for hij and gij ! use Sub, only: remove_mean_value ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! endsubroutine special_before_boundary !*********************************************************************** subroutine special_after_boundary(f) ! ! Possibility to modify the f array after the boundaries are ! communicated. ! ! 07-aug-17/axel: coded ! real, dimension (mx,my,mz,mfarray), intent(inout) :: f ! endsubroutine special_after_boundary !*********************************************************************** subroutine special_after_timestep(f,df,dt_,llast) ! ! Possibility to modify the f and df after df is updated. ! Used for the Fargo shift, for instance. ! ! 27-nov-08/wlad: coded ! logical, intent(in) :: llast real, dimension(mx,my,mz,mfarray), intent(inout) :: f real, dimension(mx,my,mz,mvar), intent(inout) :: df real, intent(in) :: dt_ ! ! Compute the transverse part of the stress tensor by going into Fourier space. ! if (llast) call compute_gT_and_gX_from_gij(f,'St') ! call keep_compiler_quiet(df) call keep_compiler_quiet(dt_) ! endsubroutine special_after_timestep !*********************************************************************** subroutine compute_gT_and_gX_from_gij(f,label) ! ! Compute the transverse part of the stress tensor by going into Fourier space. ! ! 07-aug-17/axel: coded ! use Fourier, only: fourier_transform, fft_xyz_parallel use SharedVariables, only: put_shared_variable ! real, dimension (:,:,:,:), allocatable :: Tpq_re, Tpq_im real, dimension (:,:,:), allocatable :: S_T_re, S_T_im, S_X_re, S_X_im real, dimension (:), allocatable :: kx, ky, kz real, dimension (mx,my,mz,mfarray) :: f real, dimension (6) :: Pij, e_T, e_X, Sij_re, Sij_im real, dimension (3) :: e1, e2 ! real, dimension(nk) :: specGWs ,specGWh ,specStr ! real, dimension(nk) :: specGWshel,specGWhhel,specStrhel integer :: i,j,p,q,ik,ikx,iky,ikz,stat,ij,pq,ip,jq,jStress_ij logical :: lscale_tobox1=.true. real :: kscale_factor, fact, sign_switch real :: ksqr, one_over_k2, k1, k2, k3, k1sqr, k2sqr, k3sqr real :: hhTre, hhTim, hhXre, hhXim, coefAre, coefAim real :: ggTre, ggTim, ggXre, ggXim, coefBre, coefBim real :: cosot, sinot, om12, om1, om, omt1 intent(inout) :: f character (len=2) :: label ! ! Check that the relevant arrays are registered ! if (.not.lStress_as_aux.and.label=='St') call fatal_error('compute_gT_and_gX_from_gij','lStress_as_aux must be true') ! ! For testing purposes, if lno_transverse_part=T, we would not need to ! compute the Fourier transform, so we would skip the rest. ! ! Allocate memory for arrays. ! allocate(S_T_re(nx,ny,nz),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for S_T_re') ! allocate(S_T_im(nx,ny,nz),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for S_T_im') ! allocate(S_X_re(nx,ny,nz),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for S_X_re') ! allocate(S_X_im(nx,ny,nz),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for S_X_im') ! allocate(Tpq_re(nx,ny,nz,6),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for Tpq_re') ! allocate(Tpq_im(nx,ny,nz,6),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij','Could not allocate memory for Tpq_im') ! allocate(kx(nxgrid),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij', & 'Could not allocate memory for kx') allocate(ky(nygrid),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij', & 'Could not allocate memory for ky') allocate(kz(nzgrid),stat=stat) if (stat>0) call fatal_error('compute_gT_and_gX_from_gij', & 'Could not allocate memory for kz') ! ! calculate k^2 ! kscale_factor=1 if (lscale_tobox1) kscale_factor=2*pi/Lx kx=cshift((/(i-(nxgrid+1)/2,i=0,nxgrid-1)/),+(nxgrid+1)/2)*kscale_factor ! kscale_factor=1 if (lscale_tobox1) kscale_factor=2*pi/Ly ky=cshift((/(i-(nygrid+1)/2,i=0,nygrid-1)/),+(nygrid+1)/2)*kscale_factor ! kscale_factor=1 if (lscale_tobox1) kscale_factor=2*pi/Lz kz=cshift((/(i-(nzgrid+1)/2,i=0,nzgrid-1)/),+(nzgrid+1)/2)*kscale_factor ! ! Set k^2 array. Note that in Fourier space, kz is the fastest index and has ! the full nx extent (which, currently, must be equal to nxgrid). ! But call it one_over_k2. ! ! Assemble stress, Tpq ! Tpq_re=0.0 Tpq_im=0.0 do ij=1,6 if (label=='St') then jStress_ij=iStress_ij-1+ij Tpq_re(:,:,:,ij)=f(l1:l2,m1:m2,n1:n2,jStress_ij) endif enddo ! ! Fourier transform all 6 components ! do ij=1,6 call fourier_transform(Tpq_re(:,:,:,ij),Tpq_im(:,:,:,ij)) !-- call fft_xyz_parallel(Tpq_re(:,:,:,ij),Tpq_im(:,:,:,ij)) enddo ! ! Set ST=SX=0 and reset all spectra. ! S_T_re=0. ; S_T_im=0. S_X_re=0. ; S_X_im=0. specGWs=0.; specGWshel=0. specGWh=0.; specGWhhel=0. specStr=0.; specStrhel=0. ! ! P11, P22, P33, P12, P23, P31 ! do iky=1,nz do ikx=1,ny do ikz=1,nx ! ! compute e_T and e_X; determine first preferred direction, ! which is a component with the smallest component by modulus. ! k1=kx(ikx+ipy*ny) k2=ky(iky+ipz*nz) k3=kz(ikz+ipx*nx) k1sqr=k1**2 k2sqr=k2**2 k3sqr=k3**2 ksqr=k1sqr+k2sqr+k3sqr ! ! find two vectors e1 and e2 to compute e_T and e_X ! if (lroot.and.ikx==1.and.iky==1.and.ikz==1) then e1=0. e2=0. Pij(1)=0. Pij(2)=0. Pij(3)=0. Pij(4)=0. Pij(5)=0. Pij(6)=0. one_over_k2=0. else one_over_k2=1./ksqr if(abs(k1)