Commit 81f5fed5 authored by Eric Buchlin's avatar Eric Buchlin

Added parallelized version

parent 6208ced0
.PHONY: all clean outclean run mpirun mpirunlsf
# Platform-dependent variables
include Makefile.platform
# Module alternatives for current model
include Makefile.model
MODULES = diagnostics.o \
boundary-$(BOUNDARY).o \
initialize.o \
ioparams.o iosimstate-$(IODATA).o iodata-$(IODATA).o \
nonlin-$(NONLIN).o \
propagation-$(PROPAGATION).o \
shell-atm.o \
types.o \
unlimitstack.o
# default number of processors for interactive MPI runs
NP=4
all: shell-atm
params: params.txt
shell-atm: $(MODULES)
$(MPIFC) $(MPIFFLAGS) -o shell-atm $(MODULES)
clean: outclean
rm -f *~ shell-atm *.o *.mod
outclean:
rm -f out* nohup.out
# default run
run: mpirun
# Parallel run with no queuing system
mpirun: shell-atm outclean
@echo "Parallel run on $(NP) processors (can be overriden by make NP=n mpirun)"
@echo "Press Enter if OK, or press Ctrl-C"
@read dummy
$(MPIRUN) -np $(NP) ./shell-atm | tee output.log
# LSF queuing system (cosmos.jpl.nasa.gov, clx.cineca.it...)
mpirunlsf: shell-atm outclean
@echo "--------------------- mpi_lsf.txt -----------------------"
@cat mpi_lsf.txt
@echo "---------------------------------------------------------"
@echo "Press Enter if mpi_lsf.txt is OK, or press Ctrl-C"
@read dummy
bsub < mpi_lsf.txt
# PBS queuing system (grid1.arcetri.astro.it...)
mpirunpbs: shell-atm outclean
@echo "--------------------- mpi_pbs.txt -----------------------"
@cat mpi_pbs.txt
@echo "---------------------------------------------------------"
@echo "Press Enter if mpi_pbs.txt is OK, or press Ctrl-C"
@read dummy
qsub < mpi_pbs.txt
# LoadLeveler queuing system
mpirunll: shell-atm outclean
@echo "---------------------- mpi_ll.txt -----------------------"
@cat mpi_ll.txt
@echo "---------------------------------------------------------"
@echo "Press Enter if mpi_ll.txt is OK, or press Ctrl-C"
@read dummy
llsubmit mpi_ll.txt
# Interactive run (without MPI)
interactiverun: shell-atm outclean
./shell-atm | tee output.log
# Module dependencies
diagnostics.mod: diagnostics.o
boundary.mod: boundary-$(BOUNDARY).o
initialize.mod: initialize.o
iodata.mod: iodata-$(IODATA).o
ioparams.mod: ioparams.o
iosimstate.mod: iosimstate-$(IODATA).o
nonlin.mod: nonlin-$(NONLIN).o
propagation.mod: propagation-$(PROPAGATION).o
types.mod: types.o
boundary-photoopen.o boundary-photosphere.o boundary-photospheres.o boundary-zpopen.o: types.mod diagnostics.mod
boundary-closed.o boundary-periodic.o: types.mod
diagnostics.o: types.mod nonlin.mod
initialize.o: types.mod diagnostics.mod
iodata-formatted.o iodata-unformatted.o: types.mod diagnostics.mod iosimstate.mod
ioparams.o: types.mod diagnostics.mod boundary.mod \
initialize.mod iodata.mod iosimstate.mod nonlin.mod propagation.mod
iosimstate-formatted.o iosimstate-unformatted.o: types.mod
nonlin-none.o nonlin-goy.o nonlin-nlshell.o: types.mod
propagation-fromm.o propagation-laxw: types.mod
shell-atm.o: types.mod nonlin.mod unlimitstack.o \
initialize.mod diagnostics.mod iodata.mod ioparams.mod iosimstate.mod
# The generic rules for these targets (modules with no alternative)
# are not enough for GNU-make because it thinks that .mod
# are Modula-2 programs!
diagnostics.o: diagnostics.f90
$(MPIFC) $(MPIFFLAGS) -c $<
initialize.o: initialize.f90
$(MPIFC) $(MPIFFLAGS) -c $<
ioparams.o: ioparams.f90
$(MPIFC) $(MPIFFLAGS) -c $<
types.o: types.f90
$(MPIFC) $(MPIFFLAGS) -c $<
# Generic rules
.SUFFIXES: .f90 .c .o
.f90.o:
$(MPIFC) $(MPIFFLAGS) -c $<
.c.o:
$(CC) -c $<
## Alternative for atmosphere
ATMOSPHERE = unif
## Alternative for forcing or boundary conditions
BOUNDARY = photospheres
## Alternative for over-expansion
FEXPAND = none
## Alternative for non-linear terms
NONLIN = goy
## Alternative for propagation
PROPAGATION = laxw
## Alternative for units of variables
UNIT = CGS
## Formatted or unformatted (binary) output
IODATA = formatted
CC = gcc
FC = mpif90
FFLAGS = -O2
MPIFC = $(FC)
MPIFFLAGS = $(FFLAGS)
MPIRUN = mpirun
# for debug:
#MPIFOPTS = -check all -warn all -warn nounused -stand f95
!+
!***********************************************************
! Module boundary-photospheres
!***********************************************************
! Variables and routines for boundary conditions and forcing
! For a photosphere on both sides of the box
!
! Authors:
! AV Andrea Verdini
! EB Eric Buchlin
!
! Modifications history:
! 07 Nov 05 EB Created
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module boundary
use types
!-
implicit none
! Identification of module alternative
character (len=*), parameter :: module_boundary = "photospheres"
contains
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Boundary conditions (reflection and forcing)
! A velocity is imposed, which leads to partial reflection
! MPI: OK (done)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine force (zp, zm, k, p, t)
use diagnostics
!-
implicit none
type(params), intent(inout) :: p
complex(kind=dfloat), intent(inout), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
complex(kind=dfloat), intent(in), dimension(p%nmin:p%nmax) :: k
real(kind=dfloat), intent(in) :: t
! boundary conditions on velocity (forcing)
complex(kind=dfloat), dimension(:), save, allocatable :: u1, u2
! amplitudes of forcing velocity for the 3 forced modes
real(kind=dfloat), dimension(0:nkforce-1) :: tmp1
real(kind=dfloat) :: tmp
! are u1 and u2 already allocated?
logical, save :: isallocated = .FALSE.
! MPI variables
integer :: status(mpi_status_size)
if (.NOT. isallocated) then
allocate (u1(p%nmin:p%nmax), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'
allocate (u2(p%nmin:p%nmax), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'
u1 = 0; u2 = 0
isallocated = .TRUE.
end if
! p%tstar is the time after which the phase of one component
! of forcing should be changed
if (modulo (t, p%tstar) .lt. p%delta_t) then
! change complex phase of 1st component, for both boundaries
call random_number (harvest=tmp1)
! indices must match the size of tmp1
if (p%mpi%me .eq. 0) p%fc%a1 = exp (cmplx (0._dfloat, tmp1 * pi * 2._dfloat))
call random_number (harvest=tmp1)
if (p%mpi%me .eq. p%mpi%np-1) p%fc%a2 = exp (cmplx (0._dfloat, tmp1 * pi * 2._dfloat))
else if (modulo (t + p%tstar/2._dfloat, p%tstar) .lt. p%delta_t) then
! change complex phase of 2nd component, for both boundaries
call random_number (harvest=tmp1)
if (p%mpi%me .eq. 0) p%fc%b1 = exp (cmplx (0._dfloat, tmp1 * pi * 2._dfloat))
call random_number (harvest=tmp1)
if (p%mpi%me .eq. p%mpi%np-1) p%fc%b2 = exp (cmplx (0._dfloat, tmp1 * pi * 2._dfloat))
end if
! set forcing (velocity field at boundaries)
u1(2:2+nkforce-1) = (sin (pi * t / p%tstar) ** 2 * p%fc%a1 + &
sin (pi * (t / p%tstar + .5_dfloat)) ** 2 * p%fc%b1) * p%forcamp
! same kind of forcing on 2nd boundary
u2(2:2+nkforce-1) = (sin (pi * t / p%tstar) ** 2 * p%fc%a2 + &
sin (pi * (t / p%tstar + .5_dfloat)) ** 2 * p%fc%b2) * p%forcamp
! apply this forcing by imposing the velocity field,
! and compute power of forcing
if (p%mpi%me .eq. 0) then
zp(:, 0) = -zm(:, 0) + 2._dfloat * u1
p%forcp = energy_tot1d (zp(:,0), p) - &
energy_tot1d (zm(:,0), p)
p%forcp = p%forcp * p%b0
end if
if (p%mpi%me .eq. p%mpi%np - 1) then
zm(:, p%nz-1) = -zp(:, p%nz-1) + 2._dfloat * u2
tmp = energy_tot1d (zm(:,p%nz-1), p) - &
energy_tot1d (zp(:,p%nz-1), p)
tmp = tmp * p%b0
if (p%mpi%np .ne. 1) then
call mpi_send (tmp, 1, mpi_double_precision, &
0, 11, mpi_comm_world, ierr)
end if
end if
if (p%mpi%me .eq. 0 .and. p%mpi%np .ne. 1) then
call mpi_recv (tmp, 1, mpi_double_precision, &
p%mpi%np-1, 11, mpi_comm_world, status, ierr)
end if
if (p%mpi%me .eq. 0) then
p%forcp = (p%forcp + tmp) * pi ** 3 / p%k0 ** 2
end if
end subroutine force
end module boundary
!+
!***********************************************************
! Module diagnostics
!***********************************************************
! Computation of diagnostics from fields, integrated quantities, conversions
!
! Authors:
! AV Andrea Verdini
! EB Eric Buchlin
!
! Modifications history:
! 22 Oct 05: EB Created (moved from shell-atm)
! 05 Nov 05: EB Take density into account
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module diagnostics
use mpi
use types
!-
implicit none
! Identification of module alternative
character (len=*), parameter :: module_diagnostics = "default"
contains
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Return the total energy (massic density) contained in shells
! of a loop cross-section
! As we don't know the density in this section or its position,
! we can't take it into account right here
! MPI: OK (no need)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function energy_tot1d (z, p)
!-
implicit none
type (params), intent(in) :: p
real (kind=dfloat) :: energy_tot1d
complex(kind=dfloat), intent(in), dimension(:) :: z
energy_tot1d = real (dot_product (z, z), kind=dfloat) * 0.25_dfloat
end function energy_tot1d
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Return the total energy contained in shells of a field (whole
! box, all processors)
! Take into acount mass density, and separation between planes
! MPI: OK (done)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function energy_tot (z, p)
!-
implicit none
type (params), intent(in) :: p
real (kind=dfloat) :: energy_tot, tmpe
complex(kind=dfloat), intent(in), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: z
integer :: i
! compute total energy contained in shells of a field
! (all planes in one processor)
tmpe = 0._dfloat
do i = 0, p%nz - 1
! take into account effect of boundary conditions (reflection):
! the boundary planes should be counted half
if (((i .eq. 0 ) .and. (p%mpi%me .eq. 0 )) .or. &
&((i .eq. p%nz-1) .and. (p%mpi%me .eq. p%mpi%np-1))) then
tmpe = tmpe + energy_tot1d (z(:, i), p) * p%rho0 * .5_dfloat
else
tmpe = tmpe + energy_tot1d (z(:, i), p) * p%rho0
end if
end do
! mpi_allreduce is overkill if we need to know the total energy only
! in the root processor, but we don't know if it is the case
! (unless we write another procedure for this case)
call mpi_allreduce (tmpe, energy_tot, 1, mpi_double_precision, &
mpi_sum, mpi_comm_world, ierr)
! convert this volumic energy to an energy
energy_tot = energy_tot * p%dz * pi ** 3 / p%k0 ** 2
end function energy_tot
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Return the total dissipation for a field u or b (whole
! box, all processors)
! Argument field is u or b
! Take into acount mass density, and separation between planes
! MPI: OK (done)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function diss_tot (z, k, p, disscoeff)
!-
implicit none
type (params), intent(in) :: p
complex(kind=dfloat), intent(in), dimension(p%nmin:p%nmax) :: k
real (kind=dfloat), intent(in) :: disscoeff
real (kind=dfloat) :: diss_tot, tmpd
complex(kind=dfloat), intent(in), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: z
integer :: i
! compute total dissipation contained in shells of a field
! (all planes of one processor)
tmpd = 0._dfloat
do i = 0, p%nz - 1
! take into account effect of boundary conditions (reflexion)
! see also energy_tot
if (((i .eq. 0 ) .and. (p%mpi%me .eq. 0 )) .or. &
&((i .eq. p%nz-1) .and. (p%mpi%me .eq. p%mpi%np-1))) then
tmpd = tmpd + energy_tot1d (z(:, i) * k, p) * p%rho0 * 2._dfloat
else
tmpd = tmpd + energy_tot1d (z(:, i) * k, p) * p%rho0 * 4._dfloat
end if
end do
! mpi_allreduce is overkill if we need to know the total energy only
! in the root processor, but we don't know if it is the case
! (unless we write another procedure for this case)
call mpi_allreduce (tmpd, diss_tot, 1, mpi_double_precision, &
mpi_sum, mpi_comm_world, ierr)
! convert this volumic energy to an energy
diss_tot = disscoeff * diss_tot * p%dz * pi ** 3 / p%k0 ** 2
end function diss_tot
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Transform (u, b) to (zp, zm) for the planes of one processor
! MPI: OK (no need)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine ub_to_zpm (u, b, zp, zm, p)
!-
implicit none
type(params), intent(in) :: p
complex(kind=dfloat), intent(in), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: u, b
complex(kind=dfloat), intent(out), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
zp = u + b
zm = u - b
end subroutine ub_to_zpm
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Transform (zp, zm) to (u, b) for the planes of one processor
! MPI: OK (no need)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine zpm_to_ub (zp, zm, u, b, p)
!-
implicit none
type(params), intent(in) :: p
complex(kind=dfloat), intent(in), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
complex(kind=dfloat), intent(out), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: u, b
u = (zp + zm) / 2._dfloat
b = (zp - zm) / 2._dfloat
end subroutine zpm_to_ub
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Check whether time scales are shorter than time step
! MPI: OK (done)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine check_time_scales (zp, zm, k, p, ti, t, tscales)
!-
implicit none
type(params), intent(inout) :: p
complex(kind=dfloat), intent(in), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
complex(kind=dfloat), intent(in), dimension(p%nmin:p%nmax) :: k
integer , intent(in) :: ti
real(kind=dfloat) , intent(in) :: t
type(timescales) , intent(inout) :: tscales
real(kind=dfloat) :: tmp
integer :: iz
integer, save :: ic = 0
! compute smallest non-linear time scale
tmp = 0
do iz = 0, p%nz - 1
tmp = max (tmp, &
maxval (abs (k * zp(:,iz) )), &
maxval (abs (k * zm(:,iz) )))
end do
tscales%taunl = 1._dfloat / tmp
call mpi_allreduce (tscales%taunl, tmp, 1, mpi_double_precision, &
mpi_min, mpi_comm_world, ierr)
tscales%taunl = tmp ! not useful on other processors than root
if (ic .eq. 0) then ! do this only at the first call
! compute smallest dissipation time scale
tmp = abs (k(p%nmax) ** 2)
tscales%taunu = min (1._dfloat / p%nu , 1._dfloat / p%eta) / tmp
! compute smallest Alfven time scale
tscales%taual = p%dz / p%b0
! compute crossing time
tscales%taucr = p%tnz * p%dz / p%b0
end if
if (p%mpi%me .eq. 0) then
! Update time step (adaptative time step)
! CFL condition with security factor
! No constraint from taunu as dissipation is implicit
p%delta_t = min (tscales%taunl, tscales%taual) / 5._dfloat
! Time step and scales information
100 format ('ti:', i10, ' t:', ES12.5, ' dt:', ES10.3, &
' tNL:', ES10.3)
101 format (' tcr:', ES10.3, ' tnu:', ES10.3, ' tAl:', ES10.3)
if (mod (ic, 100) == 0) then
print 100, ti, t, p%delta_t, tscales%taunl
if (mod (ic, 1000) == 0) print 101, &
tscales%taucr, tscales%taunu, tscales%taual
call flush (6)
end if
endif
call mpi_bcast (p%delta_t, 1, mpi_double_precision, &
0, mpi_comm_world, ierr)
ic = ic + 1
end subroutine check_time_scales
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialize or print wallclock time needed for simulation
! Prints information on efficiency
! MPI: OK (done)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine wallclock (p, ti, init)
!-
type(params), intent(in) :: p
integer, intent(in) :: ti
logical, intent(in), optional :: init
real(kind=dfloat), save :: clock0, clock1
real(kind=dfloat) :: cd
if (p%mpi%me .eq. 0) then
if (present (init) .and. init) then
call cpu_time (clock0)
else
call cpu_time (clock1)
cd = clock1 - clock0
print '(A,ES10.3,A,ES10.3,A)', 'Simulation duration: ', cd, &
' s (', cd / ti / p%tnz * p%mpi%np, ' s/nt/nz*np )'
! (don't write / (ti * p%nz) to avoid integer overflows)
end if
end if
end subroutine wallclock
end module diagnostics
!+
!***********************************************************
! Module initialize
!***********************************************************
! Initialization of fields, of atmosphere, of wavenumbers and
! of random numbers
!
! Authors:
! AV Andrea Verdini
! EB Eric Buchlin
!
! Modifications history:
! 22 Oct 05: EB Created (moved from shell-atm)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module initialize
use types
!-
implicit none
! Identification of module alternative
character (len=*), parameter :: module_initialize = "default"
contains
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Initialization of shells values
! MPI: OK (no need, but be careful when changing initial field)
! We can also write the initial fields in simstate.dat
! by an external program (or IDL or Python script) and start from there
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine init_shells (zp, zm, k, p)
use diagnostics
!-
implicit none
type(params), intent(in) :: p
complex(kind=dfloat), intent(out), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
complex(kind=dfloat), intent(in), dimension(p%nmin:p%nmax) :: k
zp = 0._dfloat
zm = 0._dfloat
end subroutine init_shells
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Return in array k the values of the wavenumbers for each
! index value in the perpendicular direction
! Initialize random number generator
! MPI: OK (no need, if random numbers are generated for both boundaries
! when used, on each processor, even if boundaries are on different
! processors)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine get_wavenumbers (k, p)
!-
implicit none
type(params), intent(in) :: p
complex(kind=dfloat), intent(out), dimension(p%nmin:p%nmax) :: k
integer :: i, n
integer, dimension(:), allocatable :: randseed
! get size of seed and current seed and initialize to a fixed value
call random_seed (size=n)
allocate (randseed (n))
do i = 1, n
randseed(i) = i
end do
call random_seed (put=randseed)
! Initialize array of wavenumbers corresponding to each shell
k(p%nmin) = (p%k0 * (p%lambda ** p%nmin)) * (1._dfloat, 0._dfloat)
do i = p%nmin + 1, p%nmax
k(i) = k(i - 1) * p%lambda
end do
end subroutine get_wavenumbers
end module initialize
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#!/bin/bash
#@ class = clallmds
#@ job_name = shell-atm
#@ total_tasks = 64
#@ node = 4
#@ wall_clock_limit = 20:00:00
#@ output = output.$(jobid).log
#@ error = error.$(jobid).log
#@ environment = COPY_ALL
#@ job_type = mpich
#@ queue
module load gnu-env/4.7.2 openmpi/1.8.4_gnu47
mpirun -bycore -bind-to-core -report-bindings ./shell-atm
!+
!***********************************************************
! Module nonlin-goy
!***********************************************************
! Integration of non-linear terms through shell-models
!
! Authors:
! AV Andrea Verdini
! EB Eric Buchlin
!
! Modifications history:
! 24 Oct 05: EB Created (moved from shell-atm)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module nonlin
use types
!-
implicit none
! Identification of module alternative
character (len=*), parameter :: module_nonlin = "goy"
contains
!+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Perform one step of 3rd order Runge-Kutta method integration
! Tesi di laurea, Marco Ghiglione, p 38
! Version for shell-loop: deals with 2D fields, so that Alfven-wave
! propagation can be directly handled by the routine computing the
! right handside of the differential equation (derivs2d)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine integrate_step_nonlin (zp, zm, k, p)
!-
implicit none
type(params), intent(in) :: p
complex(kind=dfloat), intent(inout), &
dimension(p%nmin:p%nmax, 0:p%nz-1) :: zp, zm
complex(kind=dfloat), intent(in), dimension(p%nmin:p%nmax) :: k
! temporary arrays for numerical scheme
complex(kind=dfloat), dimension(:,:), allocatable, save :: vp, vm, &
gp, gm, tp, tm
! are these temporary arrays already allocated?
logical, save :: isallocated = .FALSE.
integer :: aerr
if (.NOT. isallocated) then
allocate (vp(p%nmin:p%nmax, 0:p%nz-1), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'
allocate (vm(p%nmin:p%nmax, 0:p%nz-1), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'
allocate (gp(p%nmin:p%nmax, 0:p%nz-1), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'
allocate (gm(p%nmin:p%nmax, 0:p%nz-1), stat=aerr)
if (aerr .ne. 0) stop 'Allocation error'