MODULE INITIALIZE implicit none ! -------READING IN BASIC PARAMETER FILES-------! ! --params.i : contains the simulation parameters ! --fftw3.f : The fftw parameter file ! --mpif.h : The MPI parameter file INCLUDE 'params.i' INCLUDE 'fftw3.f' ! INCLUDE '/opt/openmpi-ifort/1.3.3-x86_64/include/mpif.h' INCLUDE '/opt/openmpi-ifort/1.4.3-x86_64/include/mpif.h' ! INCLUDE '/usr/lpp/ppe.poe/include/thread64/mpif.h' ! INCLUDE THE LOCATION OF THE mpif.h FILE ! INCLUDE '/usr/lpp/ppe.poe/include/thread64/mpif.h' ! ---------------------------- integer time,step_rk,dimfive,option integer MAXTIME,e_rad,num_steps, steps, o_rad, freq_filtering integer*8 fftw_plan_fwd_x,fftw_plan_inv_x integer*8 fftw_plan_fwd_y,fftw_plan_inv_y real*8 pi,deltat,rsun,z(nz), x(nx), y(ny) real*8 dimc, diml, spongex(nx), dx, dy, p00(nz),& dimrho, normx, normy,nu, spongey(ny), spongez(nz),height(nz), visc(nz), cq(nz) parameter(normx=1.0/DBLE(nx), normy = 1.0/DBLE(ny)) parameter(rsun = 69598946770.0, pi = 3.14159265358979) real*8 stretch(nz), unstretch(nz), stretchx, stretchy, decay_coeff_x(nx/2+1) real*8 decay_coeff_y(ny/2+1) ! MPI Related definitions integer numtasks, rank, status(MPI_STATUS_SIZE) ! Domain distribution ! integer, dimension(:), allocatable :: dim2, ystart, fwd_recvtypes, fwd_sendtypes integer, dimension(:), allocatable :: dim1, inv_recvtypes, inv_sendtypes, dimz integer, dimension(:), allocatable :: fwd_sdispls, inv_sdispls, fwd_z_displs integer, dimension(:), allocatable :: fwd_rdispls, inv_rdispls integer, dimension(:), allocatable :: fwd_z_send, inv_z_send, fwd_z_recv, inv_z_recv ! Definitions of rank and number of processes ! integer datatype, transp ! external datatype ! external transp ! END OF MPI DEFINITIONS ! Common to all calculations ! div - divergence of the velocity vector ! vr - forcing applied on the RHS of the radial momentum equation ! gamma - the adiabatic index of the background solar state ! p0, rho0, T0, c_speed, and g - pressure, density, temperature, sound speed ! and gravity of the background solar state ! c2 = c_speed^2, rhoinv = 1/rho0 ! gradp, gradvr - radial velocity and pressure gradients in 3-space ! omega - vorticities vector ! a - array that contains the 5 variables in the calculation: density, radial ! velocity, latitudinal velocity, longitudinal velocity and pressure in that order ! temp_step, scr - scratch space arrays used in the time evolution algorithm ! DOMBROSKI: "pr", "dc2" FOR Birch et al. 2009 SOUND PERTURBATION ! ! real*8 pr ! real*8, allocatable, dimension(:,:,:) :: dc2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real*8, allocatable, dimension(:,:,:) :: div, vr, tempin,c2, c2a, omega_r,spongexyz real*8, allocatable, dimension(:,:) :: forcing real*8, allocatable, dimension(:,:,:) :: p0, gradp0_z, c_speed, rho0, gradrho0_z, rhoinv, c2rho0, c2div_v0 ,reduction, gamma, g real*8, dimension(nz) :: source_dep real*8, allocatable, target, dimension(:,:,:,:) :: gradvz,gradp real*8, allocatable, target, dimension(:,:,:,:) :: a,temp_step,scr real*8, pointer, dimension(:,:,:) :: rho,RHSv_z,dvzdz, & & RHSv_y,RHScont,p,RHSp,gradp_z,gradp_y,gradp_x,RHSv_x real*8, allocatable, target, dimension(:,:,:) :: dvxdx, dvydy ! For periodic horizontal boundaries complex*16, allocatable, dimension(:) :: eyekx, eyeky complex*16, parameter :: eye = (0.0d0,1.0d0) !-----Quadratic and Lagrange interpolation stuff------! integer time1, time_old real*8 x0, x1, x2, x3, x4, x5, x6 real*8, dimension(:,:), allocatable :: z_i, z_iplus1, LC0, LC1, LC2, LC3, LC4, LC5, LC6 !-----Flow stuf------! real*8, allocatable, target, dimension(:,:,:,:) :: omega,v0,gradrho,omega0,advect0 real*8, allocatable, dimension(:,:,:) :: flow_x, flow_y, flow_z, div_v0 real*8, pointer, dimension(:,:,:) :: v0_x, v0_y, v0_z, omega0_x, omega0_y, omega0_z, gradrho_x, & gradrho_y, gradrho_z, advect0_x, advect0_y, advect0_z, & omega_x, omega_y, omega_z !-----Magnetic field stuff-----! real*8 dimb!, reduction(nz) real*8, allocatable, dimension(:,:,:) :: box, boy, boz, gradp0_x, gradp0_y, dzbx, dzflux2, & curlbox, curlboy, curlboz, curlbx, curlby, curlbz, gradrho0_x, gradrho0_y,& dzby, dzflux1, flux1, flux2, flux3 real*8, pointer, dimension(:,:,:) :: bx, by, bz, RHSb_x, RHSb_y, RHSb_z, v_x, v_y, v_z !-----DISPLACEMENT--------! real*8, allocatable, target, dimension(:,:,:,:) :: scratch real*8, pointer, dimension(:,:,:) :: dxixdx, dxiydy, dxizdz, RHSxi_x, RHSxi_y, RHSxi_z, & xi_x, xi_y, xi_z !------PML STUFF------! real*8, allocatable, dimension(:,:,:) :: az, bzpml real*8, target, allocatable, dimension(:,:,:,:) :: scrpml, psipml, pmlvars real*8, pointer, dimension(:,:,:) :: psivz, psip, psidzbx, psidzby, RHSpsiinductionbx,& psiinductionbx, psiinductionby, RHSpsivz, RHSpsip, & RHSpsidzbx, RHSpsidzby, RHSpsiinductionby Contains !========================================================================== Subroutine Initialize_all() implicit none integer i, j, k,fwd_sprev, inv_sprev, rem1, fwd_prev_z integer fwd_rprev, inv_rprev, ierr, rem2, rem3 integer (KIND=MPI_ADDRESS_KIND) dum,sizeofdouble real*8, allocatable, dimension(:,:,:) :: in complex*16, allocatable, dimension(:,:,:) :: co1 real*8 param1, param2, value1, value2, tempxyz, const, themax call MPI_INIT(ierr) call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr) call MPI_COMM_SIZE(MPI_COMM_WORLD, numtasks, ierr) call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_PRECISION,dum,sizeofdouble,ierr) call PARSE_OPTIONS () if (.NOT. (periodic)) then do i=1,nx x(i) = DBLE(i-1.0)/DBLE(nx - 1.0) enddo do i=1,ny y(i) = DBLE(i-1.0)/DBLE(ny - 1.0) enddo else do i=1,nx x(i) = DBLE(i-1.0)/DBLE(nx) enddo do i=1,ny y(i) = DBLE(i-1.0)/DBLE(ny) enddo allocate(eyekx(nx/2+1), eyeky(ny/2+1)) do i=1,nx/2+1 eyekx(i) = DBLE(i-1.0)*eye*normx*2.0*pi enddo do i=1,ny/2+1 eyeky(i) = DBLE(i-1.0)*eye*normy*2.0*pi enddo endif ! DECAY_COEFF'S CHANGED TO A SHARP CUTOFF FUNCTION HERE do i=1,nx/2+1 !decay_coeff_x(i) = 1.0/(1.0 + exp((i-nx/3.0)*1.))*normx if (i.le.DBLE(nx)/3.) then decay_coeff_x(i) = normx else decay_coeff_x(i) = 0.0 endif enddo do i=1,ny/2+1 !decay_coeff_y(i) = 1.0/(1.0 + exp((i-ny/3.0)*1.))*normy if (i.le.DBLE(ny)/3.) then decay_coeff_y(i) = normy else decay_coeff_y(i) = 0.0 endif enddo stretchx = rsun/xlength stretchy = rsun/ylength allocate(dim1(0:numtasks-1), dim2(0:numtasks-1), dimz(0:numtasks-1),& fwd_sendtypes(0:numtasks-1), inv_sendtypes(0:numtasks-1), & fwd_recvtypes(0:numtasks-1), inv_recvtypes(0:numtasks-1), & ystart(0:numtasks-1), fwd_sdispls(0:numtasks-1), fwd_rdispls(0:numtasks-1),& inv_sdispls(0:numtasks-1), inv_rdispls(0:numtasks-1), fwd_z_send(0:numtasks-1),& inv_z_send(0:numtasks-1), fwd_z_recv(0:numtasks-1), inv_z_recv(0:numtasks-1),& fwd_z_displs(0:numtasks-1)) dim1 = nx/numtasks dim2 = ny/numtasks dimz = nz/numtasks rem1 = nx - dim1(0)*numtasks rem2 = ny - dim2(0)*numtasks rem3 = nz - dimz(0)*numtasks do i=0,rem1-1 dim1(i) = dim1(i) +1 enddo do i=0,rem2-1 dim2(i) = dim2(i) +1 enddo do i=0,rem3-1 dimz(i) = dimz(i) +1 enddo ! Fwd transpose ! Sent packets are of dimensions dim1(receiver_rank), dim2(rank), nz ! Received packets are of dimensions dim2(sender_rank), dim1(rank), nz fwd_sprev = 1 fwd_rprev = 1 fwd_prev_z = 1 ystart(0) = 1 do j=0,numtasks-1 fwd_sendtypes(j) = transp(nx,dim1(j),dim2(rank),nz) fwd_recvtypes(j) = datatype(ny,dim2(j),dim1(rank),nz) fwd_z_send(j) = datatype_z(nx,nx,dim2(rank),dimz(j)) fwd_z_recv(j) = datatype_z(nx,nx,dim2(j),dimz(rank)) fwd_sdispls(j) = fwd_sprev fwd_sprev = fwd_sprev + dim1(j) fwd_rdispls(j) = fwd_rprev fwd_rprev = fwd_rprev + dim2(j) fwd_z_displs(j) = fwd_prev_z fwd_prev_z = fwd_prev_z + dimz(j) if (j .NE. 0) ystart(j) = ystart(j-1) + dim2(j-1) enddo ! Inv transpose ! Sent packets are of dimensions dim2(receiver_rank), dim1(rank), nz ! Received packets are of dimensions dim1(sender_rank), dim2(rank), nz inv_sprev = 1 inv_rprev = 1 do j=0,numtasks-1 inv_sendtypes(j) = transp(ny,dim2(j),dim1(rank),nz) inv_recvtypes(j) = datatype(nx,dim1(j),dim2(rank),nz) inv_z_send(j) = datatype_z(nx,nx,dim2(j),dimz(rank)) inv_z_recv(j) = datatype_z(nx,nx,dim2(rank),dimz(j)) inv_sdispls(j) = inv_sprev inv_sprev = inv_sprev + dim2(j) inv_rdispls(j) = inv_rprev inv_rprev = inv_rprev + dim1(j) enddo ! The size of the fifth dimension dimfive = 5 if (magnetic) dimfive = 8 if (displ) dimfive = 6 !----COMPUTING FFT PLANS allocate(in(nx,dim2(rank),nz), co1(nx/2+1,dim2(rank),nz)) call dfftw_plan_guru_dft_r2c(fftw_plan_fwd_x,1,nx,1,1,1,nz*dim2(rank),& & nx,(nx/2+1),in(1,1,1),co1(1,1,1),FFTW_ESTIMATE) call dfftw_plan_guru_dft_c2r(fftw_plan_inv_x,1,nx,1,1,1,nz*dim2(rank),& & (nx/2+1),nx,co1(1,1,1),in(1,1,1),FFTW_ESTIMATE) deallocate(in,co1) allocate(in(ny,dim1(rank),nz), co1(ny/2+1,dim1(rank),nz)) call dfftw_plan_guru_dft_r2c(fftw_plan_fwd_y,1,ny,1,1,1,nz*dim1(rank),& & ny,(ny/2+1),in(1,1,1),co1(1,1,1),FFTW_ESTIMATE) call dfftw_plan_guru_dft_c2r(fftw_plan_inv_y,1,ny,1,1,1,nz*dim1(rank),& & (ny/2+1),ny,co1(1,1,1),in(1,1,1),FFTW_ESTIMATE) deallocate(in,co1) !----- if (magnetic) then if (.not. displ) allocate(a(nx,dim2(rank),nz,8), temp_step(nx,dim2(rank),nz,8),& scr(nx,dim2(rank),nz,8)) if (displ) allocate(a(nx,dim2(rank),nz,6), temp_step(nx,dim2(rank),nz,6),& scr(nx,dim2(rank),nz,6),scratch(nx,dim2(rank),nz,5)) allocate(box(nx,dim2(rank),nz), boy(nx, dim2(rank), nz), boz(nx, dim2(rank), nz),& curlbox(nx,dim2(rank),nz), curlboy(nx,dim2(rank),nz), curlboz(nx,dim2(rank),nz),& curlbx(nx,dim2(rank),nz), curlby(nx,dim2(rank),nz), curlbz(nx,dim2(rank),nz),& gradrho0_x(nx, dim2(rank), nz), gradrho0_y(nx, dim2(rank), nz), gradp0_x(nx, dim2(rank), nz),& gradp0_y(nx, dim2(rank), nz),reduction(nx,dim2(rank),nz),dzbx(nx,dim2(rank),nz), & dzflux2(nx,dim2(rank),nz), dzflux1(nx,dim2(rank),nz), dzby(nx,dim2(rank),nz),& flux1(nx,dim2(rank),nz), flux2(nx,dim2(rank),nz), flux3(nx,dim2(rank),nz)) if (USE_PML) & allocate(scrpml(nx,dim2(rank),nzpml,6),pmlvars(nx,dim2(rank),nzpml,6),& psipml(nx,dim2(rank),nzpml,6)) else if (.not. displ) allocate(a(nx,dim2(rank),nz,5), temp_step(nx,dim2(rank),nz,5),& scr(nx,dim2(rank),nz,5), scrpml(nx,dim2(rank),nzpml,2), pmlvars(nx,dim2(rank),nzpml,2),& psipml(nx,dim2(rank),nzpml,2)) if (displ) allocate(a(nx,dim2(rank),nz,6), temp_step(nx,dim2(rank),nz,6),& scr(nx,dim2(rank),nz,6), scrpml(nx,dim2(rank),nzpml,2), pmlvars(nx,dim2(rank),nzpml,2),& psipml(nx,dim2(rank),nzpml,2), scratch(nx, dim2(rank),nz,2)) endif ! DOMBROSKI: "r", "dc2" FOR Birch et al. 2009 SOUND PERTURBATION ! !allocate(dc2(nx,dim2(rank),nz)) ! gamma, g for three-dimensional gamma1 distribution allocate(gamma(nx,dim2(rank),nz),g(nx,dim2(rank),nz)) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! allocate(gradp(nx,dim2(rank),nz,3),c2a(nx,dim2(rank),nz),& gradvz(nx,dim2(rank),nz,1),c2(nx,dim2(rank),nz), div(nx,dim2(rank),nz),& dvxdx(nx,dim2(rank),nz), forcing(nx,dim2(rank)), dvydy(nx,dim2(rank),nz),& spongexyz(nx, dim2(rank), nz), gradrho0_z(nx, dim2(rank), nz),& gradp0_z(nx, dim2(rank), nz), rhoinv(nx, dim2(rank), nz), rho0(nx, dim2(rank), nz), & p0(nx, dim2(rank), nz), c_speed(nx, dim2(rank), nz), c2rho0(nx,dim2(rank),nz),& z_i(nx,dim2(rank)), z_iplus1(nx,dim2(rank)), LC0(nx,dim2(rank)), LC1(nx,dim2(rank)), LC2(nx,dim2(rank)),& LC3(nx,dim2(rank)), LC4(nx,dim2(rank)), LC5(nx,dim2(rank)), LC6(nx,dim2(rank)),& az(nx,dim2(rank),nz),bzpml(nx,dim2(rank),nz)) if (flows) then allocate(v0(nx,dim2(rank),nz,3)) v0 = 0.0 v0_x => v0(:,:,:,1) v0_y => v0(:,:,:,2) v0_z => v0(:,:,:,3) if (.not. displ) then allocate(omega(nx,dim2(rank),nz,3),omega0(nx,dim2(rank),nz,3),& div_v0(nx,dim2(rank),nz), advect0(nx,dim2(rank),nz,3), gradrho(nx,dim2(rank),nz,3),& flow_x(nx,dim2(rank),nz), flow_y(nx,dim2(rank),nz), flow_z(nx,dim2(rank),nz), c2div_v0(nx,dim2(rank),nz) ) gradrho_x => gradrho(:,:,:,1) gradrho_y => gradrho(:,:,:,2) gradrho_z => gradrho(:,:,:,3) omega0_x => omega0(:,:,:,1) omega0_y => omega0(:,:,:,2) omega0_z => omega0(:,:,:,3) omega_x => omega(:,:,:,1) omega_y => omega(:,:,:,2) omega_z => omega(:,:,:,3) advect0_x => advect0(:,:,:,1) advect0_y => advect0(:,:,:,2) advect0_z => advect0(:,:,:,3) endif endif call solar_data () ! Setup Radial derivative stuff ! evaluate dz/di call dbyd2(unstretch(1),z(1),1,nz,1) !do k=1,nz ! stretch(k) = 1.0/unstretch(k) !enddo ! Parameters for writefits () ! X0 = (DBLE(ny) - 1.0)*0.5 ! scaleb = 10.0**(-9.0) * 2.06265/1.4959787 ! rarcsec = scaleb*rsun*10.0**(-2.0) ! ngood = nx*ny ! ngood1 = nx*ny/2 ! timestep = deltat*diml/dimc ! rin = z(1) ! rout = z(nz) ! Boundary condition parameters for the radial derivative. ! ibc(1) = 1 ! ibc(2) = 1 ! ibc(3) = 1 ! ibc(4) = 1 ! ibc(5) = 4 a = 0.0 pmlvars = 0.0 psivz => psipml(:,:,:,1) psip => psipml(:,:,:,2) RHSpsivz => scrpml(:,:,:,1) RHSpsip => scrpml(:,:,:,2) gradp_x => gradp(:,:,:,1) gradp_y => gradp(:,:,:,2) gradp_z => gradp(:,:,:,3) if (.not. displ) then rho => temp_step(:,:,:,1) v_x => temp_step(:,:,:,2) v_y => temp_step(:,:,:,3) v_z => temp_step(:,:,:,4) p => temp_step(:,:,:,5) dvzdz => gradvz(:,:,:,1) RHScont => scr(:,:,:,1) ! Continuity eqn. RHSv_x => scr(:,:,:,2) ! Radial momentum RHSv_y => scr(:,:,:,3) ! Latitudinal momentum RHSv_z => scr(:,:,:,4) ! Longitudinal momentum RHSp => scr(:,:,:,5) ! Pressure eqn. elseif (displ) then xi_x => temp_step(:,:,:,1) ! Continuity eqn. xi_y => temp_Step(:,:,:,2) ! Radial momentum xi_z => temp_step(:,:,:,3) ! Latitudinal momentum v_x => temp_step(:,:,:,4) ! Longitudinal momentum v_y => temp_step(:,:,:,5) ! Pressure eqn. v_z => temp_step(:,:,:,6) ! Pressure eqn. RHSxi_x => scr(:,:,:,1) ! Continuity eqn. RHSxi_y => scr(:,:,:,2) ! Radial momentum RHSxi_z => scr(:,:,:,3) ! Latitudinal momentum RHSv_x => scr(:,:,:,4) ! Longitudinal momentum RHSv_y => scr(:,:,:,5) ! Pressure eqn. RHSv_z => scr(:,:,:,6) ! Pressure eqn. dxizdz => gradvz(:,:,:,1) dxixdx => dvxdx dxiydy => dvydy p => scratch(:,:,:,2) rho => scratch(:,:,:,1) endif if (magnetic) then box = 0.0 boy = 0.0 boz = 0.0 if (.not. displ) then bx => temp_step(:,:,:,6) by => temp_step(:,:,:,7) bz => temp_step(:,:,:,8) RHSb_x => scr(:,:,:,6) ! equation for b_x RHSb_y => scr(:,:,:,7) ! equation for b_y RHSb_z => scr(:,:,:,8) ! equation for b_z endif if (USE_PML) then RHSpsidzbx => scrpml(:,:,:,3) RHSpsidzby => scrpml(:,:,:,4) RHSpsiinductionbx => scrpml(:,:,:,5) RHSpsiinductionby => scrpml(:,:,:,6) psidzbx => psipml(:,:,:,3) psidzby => psipml(:,:,:,4) psiinductionbx => psipml(:,:,:,5) psiinductionby => psipml(:,:,:,6) endif if (displ) then bx => scratch(:,:,:,3) by => scratch(:,:,:,4) bz => scratch(:,:,:,5) endif endif call MPI_BARRIER(MPI_COMM_WORLD, ierr) end Subroutine Initialize_all !============================================================================== Subroutine solar_data implicit none integer k!,q !real*8 data(nz,6),temp, zstep, z1 ! Pressure, Density, Sound Speed, Gamma, Gravity read in via fits file in ! driver.f90. Nondimensional solar radius (z) and 'stretch' (di/dz) come ! from text file below. open(7,file = file_data,form = 'formatted',status = 'old') ! read first z value into z1 do k=1,nz read(7,*) z(k), stretch(k) ! scale by 'nz' for consistency with 'dbyd2.f' definition: stretch(k) = stretch(k)/nz enddo close(7) ! OLD ALGORITHM: ! if (magnetic) then ! open(17,file = file_reduction,form = 'formatted',status = 'old') ! do k = 1,nz ! read(17,*) temp ! reduction(:,:,k) = temp ! enddo ! close(17) ! endif ! Data in the text file is arranged as ! Non-dimensional Solar radius, sound speed, density, pressure, gravity, gamma_1 ! In CGS Units ! solar mass: 1.9891d33 g ! read the full file into data(nz,6) ! do k = 1,nz ! read(7,*) data(k,:) ! enddo ! q = 1 ! dimc = data(q,2) ! dimrho = data(q,3) ! diml = rsun ! dimb = (4.0*pi*dimc**2.0*dimrho)**0.5 ! do k =1,nz ! z(k) = data(k,1) ! c_speed(:,:,k) = data(k,2)/dimc !*0.01 ! c2(:,:,k) = (data(k,2)/dimc)**2.0 !c_speed(4,5,k)**2.0 ! rho0(:,:,k) = data(k,3)/dimrho ! p0(:,:,k) = data(k,4)/(dimrho*dimc**2.0) ! g(k) = data(k,5)*diml/dimc**2.0 ! gamma(k) = data(k,6) ! enddo ! p00 = p0(1,1,:) height = (z-1.)*695.98994 end Subroutine solar_data !========================================================================== SUBROUTINE PARSE_OPTIONS implicit none ! USE_PML OPTIONS ARE ALL EVEN ! if ((.not. magnetic) .and. (.not. flows)) then if ((.not. TEST_IN_2D) .and. (.not. USE_PML) .and. (.not. displ)) option = 1 if ((.not. TEST_IN_2D) .and. (USE_PML) .and. (.not. displ)) option = 2 if ((.not. TEST_IN_2D) .and. (.not. USE_PML) .and. displ) option = 3 if ((.not. TEST_IN_2D) .and. (USE_PML) .and. displ) option = 4 if ((TEST_IN_2D) .and. (.not. USE_PML)) option = 5 if ((TEST_IN_2D) .and. (USE_PML)) option = 6 endif if (magnetic) then if ((.not. TEST_IN_2D) .and. (.not. USE_PML) .and. (.not. displ)) option = 7 if ((.not. TEST_IN_2D) .and. (USE_PML) .and. (.not. displ)) option = 8 if ((.not. TEST_IN_2D) .and. (.not. USE_PML) .and. displ) option = 9 if ((.not. TEST_IN_2D) .and. (USE_PML) .and. displ) option = 10 if ((TEST_IN_2D) .and. (.not. USE_PML)) option = 11 if ((TEST_IN_2D) .and. (USE_PML)) option = 12 endif if (flows) then if (.not. USE_PML .and. (.not. displ)) option = 13 if (USE_PML .and. (.not. displ)) option = 14 if (.not. USE_PML .and. displ) option = 15 if (USE_PML .and. displ) option = 16 endif if (rank==0) print *, 'Option = ', option END SUBROUTINE PARSE_OPTIONS !========================================================================== FUNCTION TRANSP_COMP(sizeofcol,m,n,p) implicit none integer transp_comp,ierr,row,m,n,sizeofcol,p,transp2d integer(KIND=MPI_ADDRESS_KIND) extent2d, dum, extent ! Datatype for transpose of an m X n X p matrix call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_COMPLEX,dum, extent,ierr) extent2d = sizeofcol*n*extent call MPI_TYPE_VECTOR(n,1,sizeofcol,MPI_DOUBLE_COMPLEX,row,ierr) call MPI_TYPE_CREATE_HVECTOR(m,1,extent,row,transp2d,ierr) call MPI_TYPE_CREATE_HVECTOR(p,1,extent2d,transp2d,transp_comp,ierr) call MPI_TYPE_COMMIT(transp_comp,ierr) END FUNCTION TRANSP_COMP !========================================================================== FUNCTION DATATYPE_COMP(sizeofcol,m,n,p) implicit none integer datatype_comp,ierr,newtype,m,n,sizeofcol,p,datatype2d integer (KIND=MPI_ADDRESS_KIND) extentofcol, dum, extent2d, extent ! Data type for a submatrix m X n X p in a supermatrix M X n X p call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_COMPLEX,dum,extent,ierr) extentofcol = sizeofcol*extent extent2d = sizeofcol*n*extent call MPI_TYPE_VECTOR(m,1,1,MPI_DOUBLE_COMPLEX,newtype,ierr) call MPI_TYPE_CREATE_HVECTOR(n,1,extentofcol,newtype,datatype2d,ierr) call MPI_TYPE_CREATE_HVECTOR(p,1,extent2d,datatype2d,datatype_comp,ierr) call MPI_TYPE_COMMIT(datatype_comp,ierr) END FUNCTION DATATYPE_COMP !========================================================================== FUNCTION DATATYPE_Z(sizeofcol,m,n,p) implicit none integer datatype_z,ierr,newtype,m,n,sizeofcol,p,datatype2d integer (KIND=MPI_ADDRESS_KIND) extentofcol, dum, extent2d, extent ! Data type for a submatrix m X n X p in a supermatrix m X n X P call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_PRECISION,dum,extent,ierr) extentofcol = sizeofcol*extent extent2d = sizeofcol*n*extent call MPI_TYPE_VECTOR(m,1,1,MPI_DOUBLE_PRECISION,newtype,ierr) call MPI_TYPE_CREATE_HVECTOR(n,1,extentofcol,newtype,datatype2d,ierr) call MPI_TYPE_CREATE_HVECTOR(p,1,extent2d,datatype2d,datatype_z,ierr) call MPI_TYPE_COMMIT(datatype_z,ierr) END FUNCTION DATATYPE_Z !========================================================================== FUNCTION TRANSP(sizeofcol,m,n,p) implicit none integer transp,ierr,row,m,n,sizeofcol,p,transp2d integer(KIND=MPI_ADDRESS_KIND) extent2d, dum, extent ! Datatype for transpose of an m X n X p matrix call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_PRECISION, dum, extent,ierr) extent2d = sizeofcol*n*extent call MPI_TYPE_VECTOR(n,1,sizeofcol,MPI_DOUBLE_PRECISION,row,ierr) call MPI_TYPE_CREATE_HVECTOR(m,1,extent,row,transp2d,ierr) call MPI_TYPE_CREATE_HVECTOR(p,1,extent2d,transp2d,transp,ierr) call MPI_TYPE_COMMIT(transp,ierr) END FUNCTION TRANSP !========================================================================== FUNCTION DATATYPE(sizeofcol,m,n,p) implicit none integer datatype,ierr,newtype,m,n,sizeofcol,p,datatype2d integer (KIND=MPI_ADDRESS_KIND) extentofcol, dum, extent2d,extent ! Data type for a submatrix m X n X p in a supermatrix M X n X p call MPI_TYPE_GET_EXTENT(MPI_DOUBLE_PRECISION,dum, extent,ierr) extentofcol = sizeofcol*extent extent2d = sizeofcol*n*extent call MPI_TYPE_VECTOR(m,1,1,MPI_DOUBLE_PRECISION,newtype,ierr) call MPI_TYPE_CREATE_HVECTOR(n,1,extentofcol,newtype,datatype2d,ierr) call MPI_TYPE_CREATE_HVECTOR(p,1,extent2d,datatype2d,datatype,ierr) call MPI_TYPE_COMMIT(datatype,ierr) END FUNCTION DATATYPE !========================================================================== END MODULE INITIALIZE