1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8!!@LICENSE
9!
10! ==================================================================
11! Allocation, reallocation, and deallocation utility routines
12! Written by J.M.Soler. May 2000.
13! ==================================================================
14! SUBROUTINE alloc_default( old, new, restore, &
15!                           copy, shrink, imin, routine )
16!   Sets defaults for allocation
17! INPUT (optional):
18!   type(allocDefaults) restore : default settings to be restored
19!   logical             copy    : Copy old array to new array?
20!   logical             shrink  : Reduce array size?
21!   integer             imin    : First index (typically 1 in Fortan,
22!                                              0 in C)
23!   character(len=*)    routine : Name of calling routine
24! OUTPUT (optional):
25!   type(allocDefaults) old     : default settings before the call
26!   type(allocDefaults) new     : default settings after the call
27! BEHAVIOR:
28!   All these defaults can be superseeded by optional arguments in
29!     each call to re_alloc.
30!   Initial default values: copy    = .true.
31!                           shrink  = .true.
32!                           imin    = 1
33!                           routine = 'unknown'
34!   If restore is present together with any of copy, shrink, imin, or
35!   routine, these are applied AFTER resetting the restore defaults.
36! USAGE:
37!   In order to restore the allocation defaults possibly set by the
38! calling routine, the suggested construction is:
39!   use alloc_module
40!   type(allocDefaults) oldDefaults
41!   call alloc_default( old=oldDefaults, routine=..., &
42!                       copy=..., shrink=... )
43!   call re_alloc(...)
44!   call alloc_default( restore=oldDefaults )
45! Notice that, if the restore call is skipped, the new defaults will
46! stay in effect until a new call to alloc_dafault is made.
47! ==================================================================
48! SUBROUTINE alloc_report( level, unit, file, printNow, threshold, shutdown )
49!   Sets the output file for the allocation report
50! INPUT (optional):
51!   integer      :: level     : Level (detail) of report
52!   integer      :: unit      : Output file unit
53!   character*(*):: file      : Output file name
54!   logical      :: printNow  : If present & true => print report now
55!   real(dp)     :: threshold : Memory threshold (in bytes) to print
56!                               the memory use of any given array
57!   logical      :: shutdown  : If present & true then close output unit.
58! BEHAVIOR:
59!   The detail/extent of the report increses with the value of level:
60! level=0 : no report at all (the default)
61! level=1 : only total memory peak and where it occurred
62! level=2 : detailed report created but printed only upon request
63! level=3 : detailed report printed at every new memory peak
64! level=4 : print every individual reallocation or deallocation
65!   If unit is present, alloc_report merely takes note of it for
66! future use, assuming that it has been already open outside.
67! In this case, file is not used.
68!   If unit is absent, and file is present, a file with that
69! name is open for future use.
70!   If both arguments are absent, a file named 'alloc_report'
71! is open for future use.
72!   If alloc_report is called with printNow=.true. several times in
73! a program, with the same unit or file argument, the subsequent
74! reports are written consecutively in the same file, each with a
75! time stamp header.
76!   If threshold is not present, threshold=0 is assumed.
77!   When performing library invocations of this module in say
78! external methods it is vital to close units before exiting the
79! library to ensure that no dangling files are kept open.
80! To ensure the file is closed before return one can use
81! the shutdown=.true. argument to forcefully close (and free)
82! the unit used to report. If not present or .false. the unit
83! will be kept open.
84!   In parallel execution, the report sections that involve every
85! reallocation (levels 1, 3, and 4) are written only by node 0.
86! The section that is written upon request (level 2) is written
87! only by the node with the highest peak of memory up to that time,
88! but it contains a summary of the memory used by all other nodes.
89!   In parallel execution, the nodes that share the same file
90! system (e.g. different chip cores or NFS-connected nodes) write
91! on the same file. Otherwise they write on files with the same name
92! in their local disks.
93! ==================================================================
94! SUBROUTINE re_alloc( array, [i1min,] i1max,
95!                      [[i2min,] i2max, [[i3min,] i3max]],
96!                      name, routine, copy, shrink )
97! INPUT:
98!   integer       :: i1min        : Lower bound of first dimension
99!                                   If not present, it is fixed by
100!                                   the last call to alloc_default.
101!                                   If present and the rank is 2(3),
102!                                   then i2min(&i3min) must also be
103!                                   present
104!   integer       :: i1max        : Upper bound of first dimension
105!   integer       :: i2min,i2max  : Bounds of second dimension, if
106!                                   applicable
107!   integer       :: i3min,i3max  : Bounds of third dimension, if appl.
108!
109! INPUT (optional):
110!   character*(*) :: name         : Actual array name or a label for it
111!   character*(*) :: routine      : Name of the calling routine
112!                                     or routine section
113!   logical       :: copy         : Save (copy) contents of old array
114!                                   to new array?
115!   logical       :: shrink       : Reallocate if the new array bounds
116!                                   are contained within the old ones?
117!                                   If not present, copy and/or shrink
118!                                      are fixed by the last call to
119!                                      alloc_default.
120! INPUT/OUTPUT:
121!   TYPE, pointer :: array : Array to be allocated or reallocated.
122!                            Implemented types and ranks are:
123!                              integer,    rank 1, 2, 3
124!                              integer*8,  rank 1
125!                              real*4,     rank 1, 2, 3, 4
126!                              real*8,     rank 1, 2, 3, 4
127!                              complex*16, rank 1, 2
128!                              logical,    rank 1, 2, 3
129!                              character(len=*), rank 1
130! BEHAVIOR:
131!   Pointers MUST NOT enter in an undefined state. Before using them
132! for the first time, they must be nullified explicitly. Alternatively,
133! in f95, they can be initialized as null() upon declaration.
134!   If argument array is not associated on input, it is just allocated.
135!   If array is associated and has the same bounds (or smaller bonds
136! and shrink is false) nothing is done. Thus, it is perfectly safe and
137! efficient to call re_alloc repeatedly without deallocating the array.
138! However, subroutine dealloc is provided to eliminate large arrays
139! when they are not needed.
140!   In order to save (copy) the contents of the old array, the new array
141! needs to be allocated before deallocating the old one. Thus, if the
142! contents are not needed, or if reducing memory is a must, calling
143! re_alloc with copy=.false. makes it to deallocate before allocating.
144!   The elements that are not copied (because copy=.false. or because
145! they are outside the bounds of the input array) return with value
146! zero (integer and real), .false. (logical), or blank (character).
147!   If imin>imax for any dimension, the array pointer returns
148! associated to a zero-size array.
149!   Besides allocating or reallocating the array, re_alloc keeps a count
150! of memory usage (and prints a report in a file previously fixed by a
151! call to alloc_report). Thus, it is not recommended to call re_alloc
152! to reallocate an array not allocated by it the first time.
153!   If name is not present, the memory count associated to the
154! allocation/deallocation is still made, but the allocation report
155! will account the array under the generic name 'unknown'.
156!   The routine argument is NOT used to classify the counting of
157! memory usage. The classification uses only the name argument.
158! This is because a pointer may be allocated in one routine and
159! deallocated in a different one (i.e. when it is used to return an
160! array whose size is not known in advance). However, the routine
161! argument is reported when alloc_report=4, and it is also used to
162! report in which routine the memory peack occurs. If you want the
163! routine name to be used for classification, you should include it
164! as part of the name argument, like in name='matInvert '//'aux'.
165! ==================================================================---
166! SUBROUTINE de_alloc( array, name, routine )
167! INPUT (optional):
168!   character*(*) :: name    : Actual array name or a label for it
169!   character*(*) :: routine : Name of the calling routine
170!                                or routine section
171! INPUT/OUTPUT:
172!   TYPE, pointer :: array : Array be deallocated (same types and
173!                            kinds as in re_alloc).
174! BEHAVIOR:
175!   Besides deallocating the array, re_alloc decreases the count of
176! memory usage previously counted by re_alloc. Thus, dealloc should
177! not be called to deallocate an array not allocated by re_alloc.
178! Equally, arrays allocated or reallocated by re_alloc should be
179! deallocated by dealloc.
180! ==================================================================---
181! SUBROUTINE alloc_count( delta_size, type, name, routine )
182! INPUT:
183!   integer         :: delta_size : +/-size(array)
184!                                   + => allocation
185!                                   - => deallocation
186!   character       :: type       : 'I' => integer
187!                                   'E' => integer*8
188!                                   'R' => real*4
189!                                   'D' => real*8
190!                                   'L' => logical
191!                                   'S' => character (string)
192! INPUT (optional):
193!   character(len=*) :: name      : Actual array name or a label for it
194!   character(len=*) :: routine   : Name of the calling routine
195!                                   or routine section
196! USAGE:
197!   integer,         allocatable:: intArray(:)
198!   double precision,allocatable:: doubleArray(:)
199!   complex,         allocatable:: complexArray(:)
200!   allocate( intArray(n), doubleArray(n), complexArray(n) )
201!   call alloc_count(   +n, 'I', 'intArray',  programName )
202!   call alloc_count(   +n, 'D', 'doubleArray', programName )
203!   call alloc_count( +2*n, 'R', 'complexArray', programName )
204!   deallocate( intArray, doubleArray, complexArray )
205!   call alloc_count(   -n, 'I', 'intArray',  programName )
206!   call alloc_count(   -n, 'D', 'doubleArray', programName )
207!   call alloc_count( -2*n, 'R', 'complexArray', programName )
208!
209! ==================================================================---
210
211MODULE alloc
212
213  use precision, only: sp        ! Single precision real type
214  use precision, only: dp        ! Double precision real type
215  use precision, only: i8b       ! 8-byte integer
216  use parallel,  only: Node      ! My processor node index
217  use parallel,  only: Nodes     ! Number of parallel processors
218  use parallel,  only: ionode    ! Am I the I/O processor?
219  use parallel,  only: parallel_init  ! Initialize parallel variables
220  use sys,       only: die       ! Termination routine
221  use m_io,      only: io_assign ! Get and reserve an available IO unit
222#ifdef MPI
223  use mpi_siesta
224#endif
225
226  implicit none
227
228PUBLIC ::             &
229  alloc_default,      &! Sets allocation defaults
230  alloc_report,       &! Sets log report defaults
231  re_alloc,           &! Allocation/reallocation
232  de_alloc,           &! Deallocation
233  alloc_count,        &! Memory counting for external allocs
234  allocDefaults        ! Derived type to hold allocation defaults
235
236PRIVATE      ! Nothing is declared public beyond this point
237
238  interface de_alloc
239    module procedure &
240      dealloc_i1, dealloc_i2, dealloc_i3,             &
241      dealloc_E1,                                     &
242      dealloc_r1, dealloc_r2, dealloc_r3, dealloc_r4, &
243      dealloc_d1, dealloc_d2, dealloc_d3, dealloc_d4, &
244      dealloc_c1, dealloc_c2,                         &
245      dealloc_z1, dealloc_z2, dealloc_z3, dealloc_z4, &
246      dealloc_l1, dealloc_l2, dealloc_l3,             &
247      dealloc_s1
248  end interface
249
250  interface re_alloc
251    module procedure &
252      realloc_i1, realloc_i2, realloc_i3,             &
253      realloc_E1,                                     &
254      realloc_r1, realloc_r2, realloc_r3, realloc_r4, &
255      realloc_d1, realloc_d2, realloc_d3, realloc_d4, &
256      realloc_c1, realloc_c2,                         &
257      realloc_z1, realloc_z2, realloc_z3, realloc_z4, &
258      realloc_l1, realloc_l2, realloc_l3,             &
259      realloc_s1
260!    module procedure & ! AG: Dangerous!!!
261!      realloc_i1s, realloc_i2s, realloc_i3s,              &
262!      realloc_r1s, realloc_r2s, realloc_r3s, realloc_r4s, &
263!      realloc_d1s, realloc_d2s, realloc_d3s, realloc_d4s, &
264!      realloc_l1s, realloc_l2s, realloc_l3s
265  end interface
266
267  ! Initial default values
268  character(len=*), parameter :: &
269    DEFAULT_NAME = 'unknown_name'         ! Array name default
270  character(len=*), parameter :: &
271    DEFAULT_ROUTINE = 'unknown_routine'   ! Routine name default
272  integer, save ::               &
273    REPORT_LEVEL = 0,            &! Level (detail) of allocation report
274    REPORT_UNIT  = 0              ! Output file unit for report
275  character(len=50), save ::     &
276    REPORT_FILE = 'alloc_report'  ! Output file name for report
277  real(dp), save ::              &
278    REPORT_THRESHOLD = 0          ! Memory threshold (in bytes) to print
279                                  ! the memory use of any given array
280
281  ! Derived type to hold allocation default options
282  type allocDefaults
283    private
284    logical ::          copy = .true.              ! Copy default
285    logical ::          shrink = .true.            ! Shrink default
286    integer ::          imin = 1                   ! Imin default
287    character(len=32):: routine = DEFAULT_ROUTINE  ! Routine name default
288  end type allocDefaults
289
290  ! Object to hold present allocation default options
291  type(allocDefaults), save :: DEFAULT
292
293  ! Internal auxiliary type for a binary tree
294  type TREE
295    character(len=80)   :: name  ! Name of an allocated array
296    real(DP)            :: mem   ! Present memory use of the array
297    real(DP)            :: max   ! Maximum memory use of the array
298    real(DP)            :: peak  ! Memory use of the array during
299                                 !   peak of total memory
300    type(TREE), pointer :: left  ! Pointer to data of allocated arrays
301                                 !   preceeding in alphabetical order
302    type(TREE), pointer :: right ! Pointer to data of allocated arrays
303                                 !   trailing in alphabetical order
304  end type TREE
305
306  ! Global variables used to store allocation data
307  real(DP),   parameter     :: MBYTE = 1.e6_dp
308  type(TREE), pointer, save :: REPORT_TREE
309  real(DP),            save :: TOT_MEM  = 0._dp
310  real(DP),            save :: PEAK_MEM = 0._dp
311  character(len=80),   save :: PEAK_ARRAY = ' '
312  character(len=32),   save :: PEAK_ROUTINE = ' '
313  integer,             save :: MAX_LEN  = 0
314
315  ! Other common variables
316  integer :: IERR
317  logical :: ASSOCIATED_ARRAY, NEEDS_ALLOC, NEEDS_COPY, NEEDS_DEALLOC
318
319CONTAINS
320
321! ==================================================================
322
323SUBROUTINE alloc_default( old, new, restore,          &
324                          routine, copy, shrink, imin )
325implicit none
326type(allocDefaults), optional, intent(out) :: old, new
327type(allocDefaults), optional, intent(in)  :: restore
328character(len=*),    optional, intent(in)  :: routine
329logical,             optional, intent(in)  :: copy, shrink
330integer,             optional, intent(in)  :: imin
331
332if (present(old))     old = DEFAULT
333if (present(restore)) DEFAULT = restore
334if (present(copy))    DEFAULT%copy   = copy
335if (present(shrink))  DEFAULT%shrink = shrink
336if (present(imin))    DEFAULT%imin   = imin
337if (present(routine)) DEFAULT%routine = routine
338if (present(new))     new = DEFAULT
339
340END SUBROUTINE alloc_default
341
342! ==================================================================
343
344SUBROUTINE alloc_report( level, unit, file, printNow, threshold, &
345  shutdown)
346
347implicit none
348
349integer,          optional, intent(in) :: level, unit
350character(len=*), optional, intent(in) :: file
351logical,          optional, intent(in) :: printNow
352real(dp),         optional, intent(in) :: threshold
353logical,          optional, intent(in) :: shutdown
354
355logical :: is_open
356
357#ifdef MPI
358integer :: MPIerror
359#endif
360
361if (present(level)) then
362   REPORT_LEVEL = level
363end if
364
365if (REPORT_LEVEL <= 0) RETURN
366
367if (node == 0) then
368  if (present(unit)) then  ! Assume that unit has been open outside
369    if (unit > 0) then
370      REPORT_UNIT = unit
371      if (present(file)) then
372        REPORT_FILE = file
373      else
374        REPORT_FILE = 'unknown'
375      end if
376    end if
377  else if (present(file)) then    ! If file is the same, do nothing
378    if (file /= REPORT_FILE) then ! Check if file was open outside
379      REPORT_FILE = file
380      inquire( file=REPORT_FILE, opened=is_open, number=REPORT_UNIT )
381      if (.not.is_open) then         ! Open new file
382        call io_assign(REPORT_UNIT)
383        open( REPORT_UNIT, file=REPORT_FILE, status='unknown')
384        write(REPORT_UNIT,*) ' '  ! Overwrite previous reports
385      end if
386    end if
387  else if (REPORT_UNIT==0) then   ! No unit has been open yet
388    REPORT_FILE = 'alloc_report'
389    call io_assign(REPORT_UNIT)
390    open( REPORT_UNIT, file=REPORT_FILE, status='unknown')
391    write(REPORT_UNIT,*) ' '      ! Overwrite previous reports
392  end if
393end if
394
395#ifdef MPI
396! Distribute information to other nodes and open REPORT_UNIT
397! NP:
398!   I am not too happy about this
399!   This forces a certain unit to be Bcasted to the others.
400!   If that unit is already open on the other nodes then you will have
401!   this module and some other module write to the same file.
402!   Perhaps we should just open the file from this node with
403!   the node id appended?
404call MPI_Bcast(REPORT_FILE,50,MPI_character,0,MPI_Comm_World,MPIerror)
405! JMS: open file only in node 0
406!if (node > 0) then
407!  open( REPORT_UNIT, file=REPORT_FILE )
408!end if
409#endif
410
411if (present(threshold)) REPORT_THRESHOLD = threshold
412
413if (present(printNow)) then
414  if (printNow) call print_report( )
415end if
416
417if (present(shutdown)) then
418   if ( shutdown .and. REPORT_UNIT /= 0 ) then
419      inquire( REPORT_UNIT, opened=is_open )
420      if ( is_open ) call io_close(REPORT_UNIT)
421   end if
422end if
423
424END SUBROUTINE alloc_report
425
426! ==================================================================
427! Integer array reallocs
428! ==================================================================
429
430SUBROUTINE realloc_i1( array, i1min, i1max, &
431                       name, routine, copy, shrink )
432! Arguments
433implicit none
434integer, dimension(:),      pointer    :: array
435integer,                    intent(in) :: i1min
436integer,                    intent(in) :: i1max
437character(len=*), optional, intent(in) :: name
438character(len=*), optional, intent(in) :: routine
439logical,          optional, intent(in) :: copy
440logical,          optional, intent(in) :: shrink
441
442! Internal variables and arrays
443character, parameter           :: type='I'
444integer, parameter             :: rank=1
445integer, dimension(:), pointer :: old_array
446integer, dimension(2,rank)     :: b, c, new_bounds, old_bounds
447
448! Get old array bounds
449ASSOCIATED_ARRAY = associated(array)
450if (ASSOCIATED_ARRAY) then
451  old_array => array          ! Keep pointer to old array
452  old_bounds(1,:) = lbound(old_array)
453  old_bounds(2,:) = ubound(old_array)
454end if
455
456! Copy new requested array bounds
457new_bounds(1,:) = (/ i1min /)
458new_bounds(2,:) = (/ i1max /)
459
460! Find if it is a new allocation or a true reallocation,
461! and if the contents need to be copied (saved)
462! Argument b returns common bounds
463! Options routine also reads common variable ASSOCIATED_ARRAY,
464! and it sets NEEDS_ALLOC, NEEDS_DEALLOC, and NEEDS_COPY
465call options( b, c, old_bounds, new_bounds, copy, shrink )
466! Deallocate old space
467if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
468  call alloc_count( -size(old_array), type, name, routine )
469  deallocate(old_array,stat=IERR)
470  call alloc_err( IERR, name, routine, old_bounds )
471end if
472
473! Allocate new space
474if (NEEDS_ALLOC) then
475  allocate( array(b(1,1):b(2,1)), stat=IERR )
476  call alloc_err( IERR, name, routine, new_bounds )
477  call alloc_count( size(array), type, name, routine )
478  array = 0
479end if
480
481! Copy contents and deallocate old space
482if (NEEDS_COPY) then
483  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
484  call alloc_count( -size(old_array), type, name, routine )
485  deallocate(old_array,stat=IERR)
486  call alloc_err( IERR, name, routine, old_bounds )
487end if
488END SUBROUTINE realloc_i1
489
490! ==================================================================
491SUBROUTINE realloc_i2( array, i1min,i1max, i2min,i2max,       &
492                       name, routine, copy, shrink )
493implicit none
494character, parameter                   :: type='I'
495integer, parameter                     :: rank=2
496integer, dimension(:,:),    pointer    :: array, old_array
497integer,                    intent(in) :: i1min, i1max, i2min, i2max
498character(len=*), optional, intent(in) :: name, routine
499logical,          optional, intent(in) :: copy, shrink
500integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
501ASSOCIATED_ARRAY = associated(array)
502if (ASSOCIATED_ARRAY) then
503  old_array => array
504  old_bounds(1,:) = lbound(old_array)
505  old_bounds(2,:) = ubound(old_array)
506end if
507new_bounds(1,:) = (/ i1min, i2min /)
508new_bounds(2,:) = (/ i1max, i2max /)
509call options( b, c, old_bounds, new_bounds, copy, shrink )
510if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
511  call alloc_count( -size(old_array), type, name, routine )
512  deallocate(old_array,stat=IERR)
513  call alloc_err( IERR, name, routine, old_bounds )
514end if
515if (NEEDS_ALLOC) then
516  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
517  call alloc_err( IERR, name, routine, new_bounds )
518  call alloc_count( size(array), type, name, routine )
519  array = 0
520end if
521if (NEEDS_COPY) then
522      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
523  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
524  call alloc_count( -size(old_array), type, name, routine )
525  deallocate(old_array,stat=IERR)
526  call alloc_err( IERR, name, routine, old_bounds )
527end if
528END SUBROUTINE realloc_i2
529! ==================================================================
530
531SUBROUTINE realloc_i3( array, i1min,i1max, i2min,i2max, i3min,i3max, &
532                       name, routine, copy, shrink )
533implicit none
534character, parameter                   :: type='I'
535integer, parameter                     :: rank=3
536integer, dimension(:,:,:),  pointer    :: array, old_array
537integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
538                                          i3min,i3max
539character(len=*), optional, intent(in) :: name, routine
540logical,          optional, intent(in) :: copy, shrink
541integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
542ASSOCIATED_ARRAY = associated(array)
543if (ASSOCIATED_ARRAY) then
544  old_array => array
545  old_bounds(1,:) = lbound(old_array)
546  old_bounds(2,:) = ubound(old_array)
547end if
548new_bounds(1,:) = (/ i1min, i2min, i3min /)
549new_bounds(2,:) = (/ i1max, i2max, i3max /)
550call options( b, c, old_bounds, new_bounds, copy, shrink )
551if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
552  call alloc_count( -size(old_array), type, name, routine )
553  deallocate(old_array,stat=IERR)
554  call alloc_err( IERR, name, routine, old_bounds )
555end if
556if (NEEDS_ALLOC) then
557  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3)),stat=IERR)
558  call alloc_err( IERR, name, routine, new_bounds )
559  call alloc_count( size(array), type, name, routine )
560  array = 0
561end if
562if (NEEDS_COPY) then
563      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3)) =  &
564  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3))
565  call alloc_count( -size(old_array), type, name, routine )
566  deallocate(old_array,stat=IERR)
567  call alloc_err( IERR, name, routine, old_bounds )
568end if
569END SUBROUTINE realloc_i3
570! ==================================================================
571SUBROUTINE realloc_E1( array, i1min, i1max, &
572                       name, routine, copy, shrink )
573! Arguments
574implicit none
575integer(i8b), dimension(:),    pointer    :: array
576integer,                    intent(in) :: i1min
577integer,                    intent(in) :: i1max
578character(len=*), optional, intent(in) :: name
579character(len=*), optional, intent(in) :: routine
580logical,          optional, intent(in) :: copy
581logical,          optional, intent(in) :: shrink
582
583! Internal variables and arrays
584character, parameter                   :: type='I'
585integer, parameter                     :: rank=1
586integer(i8b), dimension(:), pointer       :: old_array
587integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
588
589! Get old array bounds
590ASSOCIATED_ARRAY = associated(array)
591if (ASSOCIATED_ARRAY) then
592  old_array => array          ! Keep pointer to old array
593  old_bounds(1,:) = lbound(old_array)
594  old_bounds(2,:) = ubound(old_array)
595end if
596
597! Copy new requested array bounds
598new_bounds(1,:) = (/ i1min /)
599new_bounds(2,:) = (/ i1max /)
600
601! Find if it is a new allocation or a true reallocation,
602! and if the contents need to be copied (saved)
603! Argument b returns common bounds
604! Options routine also reads common variable ASSOCIATED_ARRAY,
605! and it sets NEEDS_ALLOC, NEEDS_DEALLOC, and NEEDS_COPY
606call options( b, c, old_bounds, new_bounds, copy, shrink )
607
608! Deallocate old space
609if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
610  call alloc_count( -size(old_array), type, name, routine )
611  deallocate(old_array,stat=IERR)
612  call alloc_err( IERR, name, routine, old_bounds )
613end if
614
615! Allocate new space
616if (NEEDS_ALLOC) then
617  allocate( array(b(1,1):b(2,1)), stat=IERR )
618  call alloc_err( IERR, name, routine, new_bounds )
619  call alloc_count( size(array), type, name, routine )
620  array = 0
621end if
622
623! Copy contents and deallocate old space
624if (NEEDS_COPY) then
625  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
626  call alloc_count( -size(old_array), type, name, routine )
627  deallocate(old_array,stat=IERR)
628  call alloc_err( IERR, name, routine, old_bounds )
629end if
630
631END SUBROUTINE realloc_E1
632
633! ==================================================================
634! Single precision real array reallocs
635! ==================================================================
636SUBROUTINE realloc_r1( array, i1min, i1max,        &
637                       name, routine, copy, shrink )
638implicit none
639character, parameter                   :: type='R'
640integer, parameter                     :: rank=1
641real(SP), dimension(:),     pointer    :: array, old_array
642integer,                    intent(in) :: i1min, i1max
643character(len=*), optional, intent(in) :: name, routine
644logical,          optional, intent(in) :: copy, shrink
645integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
646ASSOCIATED_ARRAY = associated(array)
647if (ASSOCIATED_ARRAY) then
648  old_array => array
649  old_bounds(1,:) = lbound(old_array)
650  old_bounds(2,:) = ubound(old_array)
651end if
652new_bounds(1,:) = (/ i1min /)
653new_bounds(2,:) = (/ i1max /)
654call options( b, c, old_bounds, new_bounds, copy, shrink )
655if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
656  call alloc_count( -size(old_array), type, name, routine )
657  deallocate(old_array,stat=IERR)
658  call alloc_err( IERR, name, routine, old_bounds )
659end if
660if (NEEDS_ALLOC) then
661  allocate( array(b(1,1):b(2,1)), stat=IERR )
662  call alloc_err( IERR, name, routine, new_bounds )
663  call alloc_count( size(array), type, name, routine )
664  array = 0._sp
665end if
666if (NEEDS_COPY) then
667  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
668  call alloc_count( -size(old_array), type, name, routine )
669  deallocate(old_array,stat=IERR)
670  call alloc_err( IERR, name, routine, old_bounds )
671end if
672END SUBROUTINE realloc_r1
673! ==================================================================
674SUBROUTINE realloc_r2( array, i1min,i1max, i2min,i2max, &
675                       name, routine, copy, shrink )
676implicit none
677character, parameter             :: type='R'
678integer, parameter               :: rank=2
679real(SP), dimension(:,:),   pointer    :: array, old_array
680integer,                    intent(in) :: i1min, i1max, i2min, i2max
681character(len=*), optional, intent(in) :: name, routine
682logical,          optional, intent(in) :: copy, shrink
683integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
684ASSOCIATED_ARRAY = associated(array)
685if (ASSOCIATED_ARRAY) then
686  old_array => array
687  old_bounds(1,:) = lbound(old_array)
688  old_bounds(2,:) = ubound(old_array)
689end if
690new_bounds(1,:) = (/ i1min, i2min /)
691new_bounds(2,:) = (/ i1max, i2max /)
692call options( b, c, old_bounds, new_bounds, copy, shrink )
693if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
694  call alloc_count( -size(old_array), type, name, routine )
695  deallocate(old_array,stat=IERR)
696  call alloc_err( IERR, name, routine, old_bounds )
697end if
698if (NEEDS_ALLOC) then
699  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
700  call alloc_err( IERR, name, routine, new_bounds )
701  call alloc_count( size(array), type, name, routine )
702  array = 0._sp
703end if
704if (NEEDS_COPY) then
705      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
706  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
707  call alloc_count( -size(old_array), type, name, routine )
708  deallocate(old_array,stat=IERR)
709  call alloc_err( IERR, name, routine, old_bounds )
710end if
711END SUBROUTINE realloc_r2
712! ==================================================================
713SUBROUTINE realloc_r3( array, i1min,i1max, i2min,i2max, i3min,i3max, &
714                       name, routine, copy, shrink )
715implicit none
716character, parameter                   :: type='R'
717integer, parameter                     :: rank=3
718real(SP), dimension(:,:,:), pointer    :: array, old_array
719integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
720                                          i3min,i3max
721character(len=*), optional, intent(in) :: name, routine
722logical,          optional, intent(in) :: copy, shrink
723integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
724ASSOCIATED_ARRAY = associated(array)
725if (ASSOCIATED_ARRAY) then
726  old_array => array          ! Keep pointer to old array
727  old_bounds(1,:) = lbound(old_array)
728  old_bounds(2,:) = ubound(old_array)
729end if
730new_bounds(1,:) = (/ i1min, i2min, i3min /)
731new_bounds(2,:) = (/ i1max, i2max, i3max /)
732call options( b, c, old_bounds, new_bounds, copy, shrink )
733if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
734  call alloc_count( -size(old_array), type, name, routine )
735  deallocate(old_array,stat=IERR)
736  call alloc_err( IERR, name, routine, old_bounds )
737end if
738if (NEEDS_ALLOC) then
739  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3)),stat=IERR)
740  call alloc_err( IERR, name, routine, new_bounds )
741  call alloc_count( size(array), type, name, routine )
742  array = 0._sp
743end if
744if (NEEDS_COPY) then
745      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3)) =  &
746  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3))
747  call alloc_count( -size(old_array), type, name, routine )
748  deallocate(old_array,stat=IERR)
749  call alloc_err( IERR, name, routine, old_bounds )
750end if
751END SUBROUTINE realloc_r3
752! ==================================================================
753SUBROUTINE realloc_r4( array, i1min,i1max, i2min,i2max, &
754                              i3min,i3max, i4min,i4max, &
755                       name, routine, copy, shrink )
756implicit none
757character, parameter                   :: type='R'
758integer, parameter                     :: rank=4
759real(SP), dimension(:,:,:,:), pointer  :: array, old_array
760integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
761                                          i3min,i3max, i4min,i4max
762character(len=*), optional, intent(in) :: name, routine
763logical,          optional, intent(in) :: copy, shrink
764integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
765ASSOCIATED_ARRAY = associated(array)
766if (ASSOCIATED_ARRAY) then
767  old_array => array
768  old_bounds(1,:) = lbound(old_array)
769  old_bounds(2,:) = ubound(old_array)
770end if
771new_bounds(1,:) = (/ i1min, i2min, i3min, i4min /)
772new_bounds(2,:) = (/ i1max, i2max, i3max, i4max /)
773call options( b, c, old_bounds, new_bounds, copy, shrink )
774if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
775  call alloc_count( -size(old_array), type, name, routine )
776  deallocate(old_array,stat=IERR)
777  call alloc_err( IERR, name, routine, old_bounds )
778end if
779if (NEEDS_ALLOC) then
780  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2), &
781                 b(1,3):b(2,3),b(1,4):b(2,4)),stat=IERR)
782  call alloc_err( IERR, name, routine, new_bounds )
783  call alloc_count( size(array), type, name, routine )
784  array = 0._sp
785end if
786if (NEEDS_COPY) then
787      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4))= &
788  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4))
789  call alloc_count( -size(old_array), type, name, routine )
790  deallocate(old_array,stat=IERR)
791  call alloc_err( IERR, name, routine, old_bounds )
792end if
793END SUBROUTINE realloc_r4
794! ==================================================================
795! Double precision real array reallocs
796! ==================================================================
797SUBROUTINE realloc_d1( array, i1min, i1max,        &
798                       name, routine, copy, shrink )
799implicit none
800character, parameter                   :: type='D'
801integer, parameter                     :: rank=1
802real(DP), dimension(:),     pointer    :: array, old_array
803integer,                    intent(in) :: i1min, i1max
804character(len=*), optional, intent(in) :: name, routine
805logical,          optional, intent(in) :: copy, shrink
806integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
807ASSOCIATED_ARRAY = associated(array)
808if (ASSOCIATED_ARRAY) then
809  old_array => array
810  old_bounds(1,:) = lbound(old_array)
811  old_bounds(2,:) = ubound(old_array)
812end if
813new_bounds(1,:) = (/ i1min /)
814new_bounds(2,:) = (/ i1max /)
815call options( b, c, old_bounds, new_bounds, copy, shrink )
816if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
817  call alloc_count( -size(old_array), type, name, routine )
818  deallocate(old_array,stat=IERR)
819  call alloc_err( IERR, name, routine, old_bounds )
820end if
821if (NEEDS_ALLOC) then
822  allocate( array(b(1,1):b(2,1)), stat=IERR )
823  call alloc_err( IERR, name, routine, new_bounds )
824  call alloc_count( size(array), type, name, routine )
825  array = 0._dp
826end if
827if (NEEDS_COPY) then
828  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
829  call alloc_count( -size(old_array), type, name, routine )
830  deallocate(old_array,stat=IERR)
831  call alloc_err( IERR, name, routine, old_bounds )
832end if
833END SUBROUTINE realloc_d1
834! ==================================================================
835SUBROUTINE realloc_d2( array, i1min,i1max, i2min,i2max, &
836                       name, routine, copy, shrink )
837implicit none
838character, parameter                   :: type='D'
839integer, parameter                     :: rank=2
840real(DP), dimension(:,:),   pointer    :: array, old_array
841integer,                    intent(in) :: i1min, i1max, i2min, i2max
842character(len=*), optional, intent(in) :: name, routine
843logical,          optional, intent(in) :: copy, shrink
844integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
845integer                                :: i1, i2
846ASSOCIATED_ARRAY = associated(array)
847if (ASSOCIATED_ARRAY) then
848  old_array => array
849  old_bounds(1,:) = lbound(old_array)
850  old_bounds(2,:) = ubound(old_array)
851end if
852new_bounds(1,:) = (/ i1min, i2min /)
853new_bounds(2,:) = (/ i1max, i2max /)
854call options( b, c, old_bounds, new_bounds, copy, shrink )
855if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
856  call alloc_count( -size(old_array), type, name, routine )
857  deallocate(old_array,stat=IERR)
858  call alloc_err( IERR, name, routine, old_bounds )
859end if
860if (NEEDS_ALLOC) then
861  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
862  call alloc_err( IERR, name, routine, new_bounds )
863  call alloc_count( size(array), type, name, routine )
864  array = 0._dp
865end if
866if (NEEDS_COPY) then
867!      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
868!  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
869  do i2 = c(1,2),c(2,2)
870  do i1 = c(1,1),c(2,1)
871    array(i1,i2) = old_array(i1,i2)
872  end do
873  end do
874  call alloc_count( -size(old_array), type, name, routine )
875  deallocate(old_array,stat=IERR)
876  call alloc_err( IERR, name, routine, old_bounds )
877end if
878END SUBROUTINE realloc_d2
879! ==================================================================
880SUBROUTINE realloc_d3( array, i1min,i1max, i2min,i2max, i3min,i3max, &
881                       name, routine, copy, shrink )
882implicit none
883character, parameter                   :: type='D'
884integer, parameter                     :: rank=3
885real(DP), dimension(:,:,:), pointer    :: array, old_array
886integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
887                                          i3min,i3max
888character(len=*), optional, intent(in) :: name, routine
889logical,          optional, intent(in) :: copy, shrink
890integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
891integer                                :: i1, i2, i3
892ASSOCIATED_ARRAY = associated(array)
893if (ASSOCIATED_ARRAY) then
894  old_array => array
895  old_bounds(1,:) = lbound(old_array)
896  old_bounds(2,:) = ubound(old_array)
897end if
898new_bounds(1,:) = (/ i1min, i2min, i3min /)
899new_bounds(2,:) = (/ i1max, i2max, i3max /)
900call options( b, c, old_bounds, new_bounds, copy, shrink )
901if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
902  call alloc_count( -size(old_array), type, name, routine )
903  deallocate(old_array,stat=IERR)
904  call alloc_err( IERR, name, routine, old_bounds )
905end if
906if (NEEDS_ALLOC) then
907  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3)),stat=IERR)
908  call alloc_err( IERR, name, routine, new_bounds )
909  call alloc_count( size(array), type, name, routine )
910  array = 0._dp
911end if
912if (NEEDS_COPY) then
913!      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3)) =  &
914!  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3))
915  do i3 = c(1,3),c(2,3)
916  do i2 = c(1,2),c(2,2)
917  do i1 = c(1,1),c(2,1)
918    array(i1,i2,i3) = old_array(i1,i2,i3)
919  end do
920  end do
921  end do
922  call alloc_count( -size(old_array), type, name, routine )
923  deallocate(old_array,stat=IERR)
924  call alloc_err( IERR, name, routine, old_bounds )
925end if
926END SUBROUTINE realloc_d3
927! ==================================================================
928SUBROUTINE realloc_d4( array, i1min,i1max, i2min,i2max, &
929                              i3min,i3max, i4min,i4max, &
930                       name, routine, copy, shrink )
931implicit none
932character, parameter                   :: type='D'
933integer, parameter                     :: rank=4
934real(DP), dimension(:,:,:,:), pointer  :: array, old_array
935integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
936                                          i3min,i3max, i4min,i4max
937character(len=*), optional, intent(in) :: name, routine
938logical,          optional, intent(in) :: copy, shrink
939integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
940ASSOCIATED_ARRAY = associated(array)
941if (ASSOCIATED_ARRAY) then
942  old_array => array
943  old_bounds(1,:) = lbound(old_array)
944  old_bounds(2,:) = ubound(old_array)
945end if
946new_bounds(1,:) = (/ i1min, i2min, i3min, i4min /)
947new_bounds(2,:) = (/ i1max, i2max, i3max, i4max /)
948call options( b, c, old_bounds, new_bounds, copy, shrink )
949if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
950  call alloc_count( -size(old_array), type, name, routine )
951  deallocate(old_array,stat=IERR)
952  call alloc_err( IERR, name, routine, old_bounds )
953end if
954if (NEEDS_ALLOC) then
955  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2), &
956                 b(1,3):b(2,3),b(1,4):b(2,4)),stat=IERR)
957  call alloc_err( IERR, name, routine, new_bounds )
958  call alloc_count( size(array), type, name, routine )
959  array = 0._dp
960end if
961if (NEEDS_COPY) then
962      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4))= &
963  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4))
964  call alloc_count( -size(old_array), type, name, routine )
965  deallocate(old_array,stat=IERR)
966  call alloc_err( IERR, name, routine, old_bounds )
967end if
968END SUBROUTINE realloc_d4
969
970! ==================================================================
971! Single precision complex array reallocs
972! ==================================================================
973SUBROUTINE realloc_c1( array, i1min, i1max,        &
974                       name, routine, copy, shrink )
975implicit none
976character, parameter                   :: type='S'
977integer, parameter                     :: rank=1
978complex(SP), dimension(:),  pointer    :: array, old_array
979integer,                    intent(in) :: i1min, i1max
980character(len=*), optional, intent(in) :: name, routine
981logical,          optional, intent(in) :: copy, shrink
982integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
983ASSOCIATED_ARRAY = associated(array)
984if (ASSOCIATED_ARRAY) then
985  old_array => array
986  old_bounds(1,:) = lbound(old_array)
987  old_bounds(2,:) = ubound(old_array)
988end if
989new_bounds(1,:) = (/ i1min /)
990new_bounds(2,:) = (/ i1max /)
991call options( b, c, old_bounds, new_bounds, copy, shrink )
992if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
993  call alloc_count( -2*size(old_array), type, name, routine )
994  deallocate(old_array,stat=IERR)
995  call alloc_err( IERR, name, routine, old_bounds )
996end if
997if (NEEDS_ALLOC) then
998  allocate( array(b(1,1):b(2,1)), stat=IERR )
999  call alloc_err( IERR, name, routine, new_bounds )
1000  call alloc_count( 2*size(array), type, name, routine )
1001  array = 0._dp
1002end if
1003if (NEEDS_COPY) then
1004  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
1005  call alloc_count( -2*size(old_array), type, name, routine )
1006  deallocate(old_array,stat=IERR)
1007  call alloc_err( IERR, name, routine, old_bounds )
1008end if
1009END SUBROUTINE realloc_c1
1010! ==================================================================
1011SUBROUTINE realloc_c2( array, i1min,i1max, i2min,i2max, &
1012                       name, routine, copy, shrink )
1013implicit none
1014character, parameter                   :: type='S'
1015integer, parameter                     :: rank=2
1016complex(SP), dimension(:,:),  pointer  :: array, old_array
1017integer,                    intent(in) :: i1min, i1max, i2min, i2max
1018character(len=*), optional, intent(in) :: name, routine
1019logical,          optional, intent(in) :: copy, shrink
1020integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1021integer                                :: i1, i2
1022ASSOCIATED_ARRAY = associated(array)
1023if (ASSOCIATED_ARRAY) then
1024  old_array => array
1025  old_bounds(1,:) = lbound(old_array)
1026  old_bounds(2,:) = ubound(old_array)
1027end if
1028new_bounds(1,:) = (/ i1min, i2min /)
1029new_bounds(2,:) = (/ i1max, i2max /)
1030call options( b, c, old_bounds, new_bounds, copy, shrink )
1031if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1032  call alloc_count( -2*size(old_array), type, name, routine )
1033  deallocate(old_array,stat=IERR)
1034  call alloc_err( IERR, name, routine, old_bounds )
1035end if
1036if (NEEDS_ALLOC) then
1037  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
1038  call alloc_err( IERR, name, routine, new_bounds )
1039  call alloc_count( 2*size(array), type, name, routine )
1040  array = 0._dp
1041end if
1042if (NEEDS_COPY) then
1043!      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
1044!  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
1045  do i2 = c(1,2),c(2,2)
1046  do i1 = c(1,1),c(2,1)
1047    array(i1,i2) = old_array(i1,i2)
1048  end do
1049  end do
1050  call alloc_count( -2*size(old_array), type, name, routine )
1051  deallocate(old_array,stat=IERR)
1052  call alloc_err( IERR, name, routine, old_bounds )
1053end if
1054END SUBROUTINE realloc_c2
1055
1056! ==================================================================
1057! Double precision complex array reallocs
1058! ==================================================================
1059SUBROUTINE realloc_z1( array, i1min, i1max,        &
1060                       name, routine, copy, shrink )
1061implicit none
1062character, parameter                   :: type='D'
1063integer, parameter                     :: rank=1
1064complex(DP), dimension(:),  pointer    :: array, old_array
1065integer,                    intent(in) :: i1min, i1max
1066character(len=*), optional, intent(in) :: name, routine
1067logical,          optional, intent(in) :: copy, shrink
1068integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1069ASSOCIATED_ARRAY = associated(array)
1070if (ASSOCIATED_ARRAY) then
1071  old_array => array
1072  old_bounds(1,:) = lbound(old_array)
1073  old_bounds(2,:) = ubound(old_array)
1074end if
1075new_bounds(1,:) = (/ i1min /)
1076new_bounds(2,:) = (/ i1max /)
1077call options( b, c, old_bounds, new_bounds, copy, shrink )
1078if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1079  call alloc_count( -2*size(old_array), type, name, routine )
1080  deallocate(old_array,stat=IERR)
1081  call alloc_err( IERR, name, routine, old_bounds )
1082end if
1083if (NEEDS_ALLOC) then
1084  allocate( array(b(1,1):b(2,1)), stat=IERR )
1085  call alloc_err( IERR, name, routine, new_bounds )
1086  call alloc_count( 2*size(array), type, name, routine )
1087  array = 0._dp
1088end if
1089if (NEEDS_COPY) then
1090  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
1091  call alloc_count( -2*size(old_array), type, name, routine )
1092  deallocate(old_array,stat=IERR)
1093  call alloc_err( IERR, name, routine, old_bounds )
1094end if
1095END SUBROUTINE realloc_z1
1096! ==================================================================
1097SUBROUTINE realloc_z2( array, i1min,i1max, i2min,i2max, &
1098                       name, routine, copy, shrink )
1099implicit none
1100character, parameter                   :: type='D'
1101integer, parameter                     :: rank=2
1102complex(DP), dimension(:,:),  pointer  :: array, old_array
1103integer,                    intent(in) :: i1min, i1max, i2min, i2max
1104character(len=*), optional, intent(in) :: name, routine
1105logical,          optional, intent(in) :: copy, shrink
1106integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1107integer                                :: i1, i2
1108ASSOCIATED_ARRAY = associated(array)
1109if (ASSOCIATED_ARRAY) then
1110  old_array => array
1111  old_bounds(1,:) = lbound(old_array)
1112  old_bounds(2,:) = ubound(old_array)
1113end if
1114new_bounds(1,:) = (/ i1min, i2min /)
1115new_bounds(2,:) = (/ i1max, i2max /)
1116call options( b, c, old_bounds, new_bounds, copy, shrink )
1117if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1118  call alloc_count( -2*size(old_array), type, name, routine )
1119  deallocate(old_array,stat=IERR)
1120  call alloc_err( IERR, name, routine, old_bounds )
1121end if
1122if (NEEDS_ALLOC) then
1123  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
1124  call alloc_err( IERR, name, routine, new_bounds )
1125  call alloc_count( 2*size(array), type, name, routine )
1126  array = 0._dp
1127end if
1128if (NEEDS_COPY) then
1129!      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
1130!  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
1131  do i2 = c(1,2),c(2,2)
1132  do i1 = c(1,1),c(2,1)
1133    array(i1,i2) = old_array(i1,i2)
1134  end do
1135  end do
1136  call alloc_count( -2*size(old_array), type, name, routine )
1137  deallocate(old_array,stat=IERR)
1138  call alloc_err( IERR, name, routine, old_bounds )
1139end if
1140END SUBROUTINE realloc_z2
1141SUBROUTINE realloc_z3( array, i1min,i1max, i2min,i2max, i3min,i3max, &
1142                       name, routine, copy, shrink )
1143implicit none
1144character, parameter                   :: type='D'
1145integer, parameter                     :: rank=3
1146complex(DP), dimension(:,:,:), pointer :: array, old_array
1147integer,                    intent(in) :: i1min, i1max, i2min, i2max
1148integer,                    intent(in) :: i3min, i3max
1149character(len=*), optional, intent(in) :: name, routine
1150logical,          optional, intent(in) :: copy, shrink
1151integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1152integer                                :: i1, i2, i3
1153ASSOCIATED_ARRAY = associated(array)
1154if (ASSOCIATED_ARRAY) then
1155  old_array => array
1156  old_bounds(1,:) = lbound(old_array)
1157  old_bounds(2,:) = ubound(old_array)
1158end if
1159new_bounds(1,:) = (/ i1min, i2min, i3min /)
1160new_bounds(2,:) = (/ i1max, i2max, i3max /)
1161call options( b, c, old_bounds, new_bounds, copy, shrink )
1162if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1163  call alloc_count( -2*size(old_array), type, name, routine )
1164  deallocate(old_array,stat=IERR)
1165  call alloc_err( IERR, name, routine, old_bounds )
1166end if
1167if (NEEDS_ALLOC) then
1168  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3)), stat=IERR )
1169  call alloc_err( IERR, name, routine, new_bounds )
1170  call alloc_count( 2*size(array), type, name, routine )
1171  array = 0._dp
1172end if
1173if (NEEDS_COPY) then
1174!      array(c(1,1):c(2,1),c(1,2):c(2,2)) =  &
1175!  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
1176  do i3 = c(1,3),c(2,3)
1177  do i2 = c(1,2),c(2,2)
1178  do i1 = c(1,1),c(2,1)
1179    array(i1,i2,i3) = old_array(i1,i2,i3)
1180  end do
1181  end do
1182  end do
1183  call alloc_count( -2*size(old_array), type, name, routine )
1184  deallocate(old_array,stat=IERR)
1185  call alloc_err( IERR, name, routine, old_bounds )
1186end if
1187END SUBROUTINE realloc_z3
1188SUBROUTINE realloc_z4( array, i1min,i1max, i2min,i2max, &
1189                              i3min,i3max, i4min,i4max, &
1190                       name, routine, copy, shrink )
1191implicit none
1192character, parameter                       :: type='D'
1193integer, parameter                         :: rank=4
1194complex(DP), dimension(:,:,:,:),  pointer  :: array, old_array
1195integer,                    intent(in)     :: i1min, i1max, i2min, i2max, &
1196                                              i3min, i3max, i4min, i4max
1197character(len=*), optional, intent(in)     :: name, routine
1198logical,          optional, intent(in)     :: copy, shrink
1199integer, dimension(2,rank)                 :: b, c, new_bounds, old_bounds
1200integer                                    :: i1, i2, i3, i4
1201ASSOCIATED_ARRAY = associated(array)
1202if (ASSOCIATED_ARRAY) then
1203  old_array => array
1204  old_bounds(1,:) = lbound(old_array)
1205  old_bounds(2,:) = ubound(old_array)
1206end if
1207new_bounds(1,:) = (/ i1min, i2min, i3min, i4min /)
1208new_bounds(2,:) = (/ i1max, i2max, i3max, i4max /)
1209call options( b, c, old_bounds, new_bounds, copy, shrink )
1210if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1211  call alloc_count( -2*size(old_array), type, name, routine )
1212  deallocate(old_array,stat=IERR)
1213  call alloc_err( IERR, name, routine, old_bounds )
1214end if
1215if (NEEDS_ALLOC) then
1216  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3),b(1,4):b(2,4)), &
1217            stat=IERR )
1218  call alloc_err( IERR, name, routine, new_bounds )
1219  call alloc_count( 2*size(array), type, name, routine )
1220  array = 0._dp
1221end if
1222if (NEEDS_COPY) then
1223!      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4)) =  &
1224!  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3),c(1,4):c(2,4))
1225  do i4 = c(1,4),c(2,4)
1226  do i3 = c(1,3),c(2,3)
1227  do i2 = c(1,2),c(2,2)
1228  do i1 = c(1,1),c(2,1)
1229    array(i1,i2,i3,i4) = old_array(i1,i2,i3,i4)
1230  end do
1231  end do
1232  end do
1233  end do
1234  call alloc_count( -2*size(old_array), type, name, routine )
1235  deallocate(old_array,stat=IERR)
1236  call alloc_err( IERR, name, routine, old_bounds )
1237end if
1238END SUBROUTINE realloc_z4
1239! ==================================================================
1240! Logical array reallocs
1241! ==================================================================
1242SUBROUTINE realloc_l1( array, i1min,i1max,  &
1243                       name, routine, copy, shrink )
1244implicit none
1245character, parameter                   :: type='L'
1246integer, parameter                     :: rank=1
1247logical, dimension(:),      pointer    :: array, old_array
1248integer,                    intent(in) :: i1min,i1max
1249character(len=*), optional, intent(in) :: name, routine
1250logical,          optional, intent(in) :: copy, shrink
1251integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1252ASSOCIATED_ARRAY = associated(array)
1253if (ASSOCIATED_ARRAY) then
1254  old_array => array
1255  old_bounds(1,:) = lbound(old_array)
1256  old_bounds(2,:) = ubound(old_array)
1257end if
1258new_bounds(1,:) = (/ i1min /)
1259new_bounds(2,:) = (/ i1max /)
1260call options( b, c, old_bounds, new_bounds, copy, shrink )
1261if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1262  call alloc_count( -size(old_array), type, name, routine )
1263  deallocate(old_array,stat=IERR)
1264  call alloc_err( IERR, name, routine, old_bounds )
1265end if
1266if (NEEDS_ALLOC) then
1267  allocate( array(b(1,1):b(2,1)), stat=IERR )
1268  call alloc_err( IERR, name, routine, new_bounds )
1269  call alloc_count( size(array), type, name, routine )
1270  array = .false.
1271end if
1272if (NEEDS_COPY) then
1273  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
1274  call alloc_count( -size(old_array), type, name, routine )
1275  deallocate(old_array,stat=IERR)
1276  call alloc_err( IERR, name, routine, old_bounds )
1277end if
1278END SUBROUTINE realloc_l1
1279! ==================================================================
1280SUBROUTINE realloc_l2( array, i1min,i1max, i2min,i2max, &
1281                       name, routine, copy, shrink )
1282implicit none
1283character, parameter                   :: type='L'
1284integer, parameter                     :: rank=2
1285logical, dimension(:,:),    pointer    :: array, old_array
1286integer,                    intent(in) :: i1min,i1max, i2min,i2max
1287character(len=*), optional, intent(in) :: name, routine
1288logical,          optional, intent(in) :: copy, shrink
1289integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1290ASSOCIATED_ARRAY = associated(array)
1291if (ASSOCIATED_ARRAY) then
1292  old_array => array
1293  old_bounds(1,:) = lbound(old_array)
1294  old_bounds(2,:) = ubound(old_array)
1295end if
1296new_bounds(1,:) = (/ i1min, i2min /)
1297new_bounds(2,:) = (/ i1max, i2max /)
1298call options( b, c, old_bounds, new_bounds, copy, shrink )
1299if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1300  call alloc_count( -size(old_array), type, name, routine )
1301  deallocate(old_array,stat=IERR)
1302  call alloc_err( IERR, name, routine, old_bounds )
1303end if
1304if (NEEDS_ALLOC) then
1305  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2)), stat=IERR )
1306  call alloc_err( IERR, name, routine, new_bounds )
1307  call alloc_count( size(array), type, name, routine )
1308  array = .false.
1309end if
1310if (NEEDS_COPY) then
1311      array(c(1,1):c(2,1),c(1,2):c(2,2)) = &
1312  old_array(c(1,1):c(2,1),c(1,2):c(2,2))
1313  call alloc_count( -size(old_array), type, name, routine )
1314  deallocate(old_array,stat=IERR)
1315  call alloc_err( IERR, name, routine, old_bounds )
1316end if
1317END SUBROUTINE realloc_l2
1318! ==================================================================
1319SUBROUTINE realloc_l3( array, i1min,i1max, i2min,i2max, i3min,i3max, &
1320                       name, routine, copy, shrink )
1321implicit none
1322character, parameter                   :: type='L'
1323integer, parameter                     :: rank=3
1324logical, dimension(:,:,:),  pointer    :: array, old_array
1325integer,                    intent(in) :: i1min,i1max, i2min,i2max, &
1326                                          i3min,i3max
1327character(len=*), optional, intent(in) :: name, routine
1328logical,          optional, intent(in) :: copy, shrink
1329integer, dimension(2,rank)             :: b, c, new_bounds, old_bounds
1330ASSOCIATED_ARRAY = associated(array)
1331if (ASSOCIATED_ARRAY) then
1332  old_array => array
1333  old_bounds(1,:) = lbound(old_array)
1334  old_bounds(2,:) = ubound(old_array)
1335end if
1336new_bounds(1,:) = (/ i1min, i2min, i3min /)
1337new_bounds(2,:) = (/ i1max, i2max, i3max /)
1338call options( b, c, old_bounds, new_bounds, copy, shrink )
1339if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1340  call alloc_count( -size(old_array), type, name, routine )
1341  deallocate(old_array,stat=IERR)
1342  call alloc_err( IERR, name, routine, old_bounds )
1343end if
1344if (NEEDS_ALLOC) then
1345  allocate( array(b(1,1):b(2,1),b(1,2):b(2,2),b(1,3):b(2,3)),stat=IERR)
1346  call alloc_err( IERR, name, routine, new_bounds )
1347  call alloc_count( size(array), type, name, routine )
1348  array = .false.
1349end if
1350if (NEEDS_COPY) then
1351      array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3)) = &
1352  old_array(c(1,1):c(2,1),c(1,2):c(2,2),c(1,3):c(2,3))
1353  call alloc_count( -size(old_array), type, name, routine )
1354  deallocate(old_array,stat=IERR)
1355  call alloc_err( IERR, name, routine, old_bounds )
1356end if
1357END SUBROUTINE realloc_l3
1358! ==================================================================
1359! Realloc routines with assumed lower bound = 1
1360!AG: Extremely dangerous -- do not use.
1361! ==================================================================
1362SUBROUTINE realloc_i1s( array, i1max, &
1363                        name, routine, copy, shrink )
1364! Arguments
1365implicit none
1366integer, dimension(:),      pointer    :: array
1367integer,                    intent(in) :: i1max
1368character(len=*), optional, intent(in) :: name
1369character(len=*), optional, intent(in) :: routine
1370logical,          optional, intent(in) :: copy
1371logical,          optional, intent(in) :: shrink
1372
1373call realloc_i1( array, DEFAULT%imin, i1max, &
1374                 name, routine, copy, shrink )
1375
1376END SUBROUTINE realloc_i1s
1377! ==================================================================
1378SUBROUTINE realloc_i2s( array, i1max, i2max,  &
1379                        name, routine, copy, shrink )
1380implicit none
1381integer, dimension(:,:),    pointer    :: array
1382integer,                    intent(in) :: i1max, i2max
1383character(len=*), optional, intent(in) :: name, routine
1384logical,          optional, intent(in) :: copy, shrink
1385call realloc_i2( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1386                 name, routine, copy, shrink )
1387END SUBROUTINE realloc_i2s
1388! ==================================================================
1389SUBROUTINE realloc_i3s( array, i1max, i2max, i3max,  &
1390                        name, routine, copy, shrink )
1391implicit none
1392integer, dimension(:,:,:),  pointer    :: array
1393integer,                    intent(in) :: i1max, i2max, i3max
1394character(len=*), optional, intent(in) :: name, routine
1395logical,          optional, intent(in) :: copy, shrink
1396call realloc_i3( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1397                 DEFAULT%imin, i3max,                             &
1398                 name, routine, copy, shrink )
1399END SUBROUTINE realloc_i3s
1400! ==================================================================
1401SUBROUTINE realloc_r1s( array, i1max, &
1402                        name, routine, copy, shrink )
1403implicit none
1404real(SP), dimension(:),     pointer    :: array
1405integer,                    intent(in) :: i1max
1406character(len=*), optional, intent(in) :: name, routine
1407logical,          optional, intent(in) :: copy, shrink
1408call realloc_r1( array, DEFAULT%imin, i1max, &
1409                 name, routine, copy, shrink )
1410END SUBROUTINE realloc_r1s
1411! ==================================================================
1412SUBROUTINE realloc_r2s( array, i1max, i2max, &
1413                        name, routine, copy, shrink )
1414implicit none
1415real(SP), dimension(:,:),   pointer    :: array
1416integer,                    intent(in) :: i1max, i2max
1417character(len=*), optional, intent(in) :: name, routine
1418logical,          optional, intent(in) :: copy, shrink
1419call realloc_r2( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1420                 name, routine, copy, shrink )
1421END SUBROUTINE realloc_r2s
1422! ==================================================================
1423SUBROUTINE realloc_r3s( array, i1max, i2max, i3max, &
1424                        name, routine, copy, shrink )
1425implicit none
1426real(SP), dimension(:,:,:), pointer    :: array
1427integer,                    intent(in) :: i1max, i2max, i3max
1428character(len=*), optional, intent(in) :: name, routine
1429logical,          optional, intent(in) :: copy, shrink
1430call realloc_r3( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1431                 DEFAULT%imin, i3max,                             &
1432                 name, routine, copy, shrink )
1433END SUBROUTINE realloc_r3s
1434! ==================================================================
1435SUBROUTINE realloc_r4s( array, i1max, i2max, i3max, i4max, &
1436                        name, routine, copy, shrink )
1437implicit none
1438real(SP), dimension(:,:,:,:), pointer  :: array
1439integer,                    intent(in) :: i1max, i2max, i3max, i4max
1440character(len=*), optional, intent(in) :: name, routine
1441logical,          optional, intent(in) :: copy, shrink
1442call realloc_r4( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1443                        DEFAULT%imin, i3max, DEFAULT%imin, i4max, &
1444                 name, routine, copy, shrink )
1445END SUBROUTINE realloc_r4s
1446! ==================================================================
1447SUBROUTINE realloc_d1s( array, i1max, &
1448                        name, routine, copy, shrink )
1449implicit none
1450real(DP), dimension(:),     pointer    :: array
1451integer,                    intent(in) :: i1max
1452character(len=*), optional, intent(in) :: name, routine
1453logical,          optional, intent(in) :: copy, shrink
1454call realloc_d1( array, DEFAULT%imin, i1max, &
1455                 name, routine, copy, shrink )
1456END SUBROUTINE realloc_d1s
1457! ==================================================================
1458SUBROUTINE realloc_d2s( array, i1max, i2max, &
1459                        name, routine, copy, shrink )
1460implicit none
1461real(DP), dimension(:,:),   pointer    :: array
1462integer,                    intent(in) :: i1max, i2max
1463character(len=*), optional, intent(in) :: name, routine
1464logical,          optional, intent(in) :: copy, shrink
1465call realloc_d2( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1466                 name, routine, copy, shrink )
1467END SUBROUTINE realloc_d2s
1468! ==================================================================
1469SUBROUTINE realloc_d3s( array, i1max, i2max, i3max, &
1470                        name, routine, copy, shrink )
1471implicit none
1472real(DP), dimension(:,:,:), pointer    :: array
1473integer,                    intent(in) :: i1max, i2max, i3max
1474character(len=*), optional, intent(in) :: name, routine
1475logical,          optional, intent(in) :: copy, shrink
1476call realloc_d3( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1477                 DEFAULT%imin, i3max,                             &
1478                 name, routine, copy, shrink )
1479END SUBROUTINE realloc_d3s
1480! ==================================================================
1481SUBROUTINE realloc_d4s( array, i1max, i2max, i3max, i4max, &
1482                        name, routine, copy, shrink )
1483implicit none
1484real(DP), dimension(:,:,:,:), pointer  :: array
1485integer,                    intent(in) :: i1max, i2max, i3max, i4max
1486character(len=*), optional, intent(in) :: name, routine
1487logical,          optional, intent(in) :: copy, shrink
1488call realloc_d4( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1489                        DEFAULT%imin, i3max, DEFAULT%imin, i4max, &
1490                 name, routine, copy, shrink )
1491END SUBROUTINE realloc_d4s
1492! ==================================================================
1493SUBROUTINE realloc_l1s( array, i1max, &
1494                        name, routine, copy, shrink )
1495implicit none
1496logical, dimension(:),      pointer    :: array
1497integer,                    intent(in) :: i1max
1498character(len=*), optional, intent(in) :: name
1499character(len=*), optional, intent(in) :: routine
1500logical,          optional, intent(in) :: copy
1501logical,          optional, intent(in) :: shrink
1502call realloc_l1( array, DEFAULT%imin, i1max, &
1503                 name, routine, copy, shrink )
1504END SUBROUTINE realloc_l1s
1505! ==================================================================
1506SUBROUTINE realloc_l2s( array, i1max, i2max, &
1507                        name, routine, copy, shrink )
1508implicit none
1509logical, dimension(:,:),    pointer    :: array
1510integer,                    intent(in) :: i1max, i2max
1511character(len=*), optional, intent(in) :: name
1512character(len=*), optional, intent(in) :: routine
1513logical,          optional, intent(in) :: copy
1514logical,          optional, intent(in) :: shrink
1515call realloc_l2( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1516                 name, routine, copy, shrink )
1517END SUBROUTINE realloc_l2s
1518! ==================================================================
1519SUBROUTINE realloc_l3s( array, i1max, i2max, i3max, &
1520                        name, routine, copy, shrink )
1521implicit none
1522logical, dimension(:,:,:),  pointer    :: array
1523integer,                    intent(in) :: i1max, i2max, i3max
1524character(len=*), optional, intent(in) :: name
1525character(len=*), optional, intent(in) :: routine
1526logical,          optional, intent(in) :: copy
1527logical,          optional, intent(in) :: shrink
1528call realloc_l3( array, DEFAULT%imin, i1max, DEFAULT%imin, i2max, &
1529                 DEFAULT%imin, i3max, name, routine, copy, shrink )
1530END SUBROUTINE realloc_l3s
1531! ==================================================================
1532! Character vector realloc
1533! ==================================================================
1534SUBROUTINE realloc_s1( array, i1min, i1max, &
1535                       name, routine, copy, shrink )
1536! Arguments
1537implicit none
1538character(len=*), dimension(:),      pointer    :: array
1539integer,                    intent(in) :: i1min
1540integer,                    intent(in) :: i1max
1541character(len=*), optional, intent(in) :: name
1542character(len=*), optional, intent(in) :: routine
1543logical,          optional, intent(in) :: copy
1544logical,          optional, intent(in) :: shrink
1545
1546! Internal variables and arrays
1547character, parameter           :: type='S'
1548integer, parameter             :: rank=1
1549character(len=len(array)), dimension(:), pointer :: old_array
1550integer, dimension(2,rank)     :: b, c, new_bounds, old_bounds
1551
1552! Get old array bounds
1553ASSOCIATED_ARRAY = associated(array)
1554if (ASSOCIATED_ARRAY) then
1555  old_array => array          ! Keep pointer to old array
1556  old_bounds(1,:) = lbound(old_array)
1557  old_bounds(2,:) = ubound(old_array)
1558end if
1559
1560! Copy new requested array bounds
1561new_bounds(1,:) = (/ i1min /)
1562new_bounds(2,:) = (/ i1max /)
1563
1564! Find if it is a new allocation or a true reallocation,
1565! and if the contents need to be copied (saved)
1566! Argument b returns common bounds
1567! Options routine also reads common variable ASSOCIATED_ARRAY,
1568! and it sets NEEDS_ALLOC, NEEDS_DEALLOC, and NEEDS_COPY
1569call options( b, c, old_bounds, new_bounds, copy, shrink )
1570
1571! Deallocate old space
1572if (NEEDS_DEALLOC .and. .not.NEEDS_COPY) then
1573  call alloc_count( -size(old_array)*len(old_array), type, name, routine )
1574  deallocate(old_array,stat=IERR)
1575  call alloc_err( IERR, name, routine, old_bounds )
1576end if
1577
1578! Allocate new space
1579if (NEEDS_ALLOC) then
1580  allocate( array(b(1,1):b(2,1)), stat=IERR )
1581  call alloc_err( IERR, name, routine, new_bounds )
1582  call alloc_count( size(array)*len(array), type, name, routine )
1583  array = ''
1584end if
1585
1586! Copy contents and deallocate old space
1587if (NEEDS_COPY) then
1588  array(c(1,1):c(2,1)) = old_array(c(1,1):c(2,1))
1589  call alloc_count( -size(old_array)*len(old_array), type, name, routine )
1590  deallocate(old_array,stat=IERR)
1591  call alloc_err( IERR, name, routine, old_bounds )
1592end if
1593
1594END SUBROUTINE realloc_s1
1595! ==================================================================
1596! Dealloc routines
1597! ==================================================================
1598SUBROUTINE dealloc_i1( array, name, routine )
1599
1600! Arguments
1601implicit none
1602integer, dimension(:),      pointer    :: array
1603character(len=*), optional, intent(in) :: name
1604character(len=*), optional, intent(in) :: routine
1605
1606if (associated(array)) then
1607  call alloc_count( -size(array), 'I', name, routine )
1608  deallocate(array,stat=IERR)
1609  call alloc_err( IERR, name, routine )
1610end if
1611
1612END SUBROUTINE dealloc_i1
1613
1614! ==================================================================
1615SUBROUTINE dealloc_i2( array, name, routine )
1616implicit none
1617integer, dimension(:,:),    pointer    :: array
1618character(len=*), optional, intent(in) :: name, routine
1619if (associated(array)) then
1620  call alloc_count( -size(array), 'I', name, routine )
1621  deallocate(array,stat=IERR)
1622  call alloc_err( IERR, name, routine )
1623
1624end if
1625END SUBROUTINE dealloc_i2
1626! ==================================================================
1627SUBROUTINE dealloc_i3( array, name, routine )
1628implicit none
1629integer, dimension(:,:,:),  pointer    :: array
1630character(len=*), optional, intent(in) :: name, routine
1631if (associated(array)) then
1632  call alloc_count( -size(array), 'I', name, routine )
1633  deallocate(array,stat=IERR)
1634  call alloc_err( IERR, name, routine )
1635
1636end if
1637END SUBROUTINE dealloc_i3
1638! ==================================================================
1639SUBROUTINE dealloc_E1( array, name, routine )
1640
1641! Arguments
1642implicit none
1643integer(i8b), dimension(:),      pointer    :: array
1644character(len=*), optional, intent(in) :: name
1645character(len=*), optional, intent(in) :: routine
1646
1647if (associated(array)) then
1648  call alloc_count( -size(array), 'I', name, routine )
1649  deallocate(array,stat=IERR)
1650  call alloc_err( IERR, name, routine )
1651
1652end if
1653
1654END SUBROUTINE dealloc_E1
1655
1656! ==================================================================
1657SUBROUTINE dealloc_r1( array, name, routine )
1658implicit none
1659real(SP), dimension(:),     pointer    :: array
1660character(len=*), optional, intent(in) :: name, routine
1661if (associated(array)) then
1662  call alloc_count( -size(array), 'R', name, routine )
1663  deallocate(array,stat=IERR)
1664  call alloc_err( IERR, name, routine )
1665end if
1666END SUBROUTINE dealloc_r1
1667! ==================================================================
1668SUBROUTINE dealloc_r2( array, name, routine )
1669implicit none
1670real(SP), dimension(:,:),   pointer    :: array
1671character(len=*), optional, intent(in) :: name, routine
1672if (associated(array)) then
1673  call alloc_count( -size(array), 'R', name, routine )
1674  deallocate(array,stat=IERR)
1675  call alloc_err( IERR, name, routine )
1676end if
1677END SUBROUTINE dealloc_r2
1678! ==================================================================
1679SUBROUTINE dealloc_r3( array, name, routine )
1680implicit none
1681real(SP), dimension(:,:,:), pointer    :: array
1682character(len=*), optional, intent(in) :: name, routine
1683if (associated(array)) then
1684  call alloc_count( -size(array), 'R', name, routine )
1685  deallocate(array,stat=IERR)
1686  call alloc_err( IERR, name, routine )
1687end if
1688END SUBROUTINE dealloc_r3
1689! ==================================================================
1690SUBROUTINE dealloc_r4( array, name, routine )
1691implicit none
1692real(SP), dimension(:,:,:,:), pointer  :: array
1693character(len=*), optional, intent(in) :: name, routine
1694if (associated(array)) then
1695  call alloc_count( -size(array), 'R', name, routine )
1696  deallocate(array,stat=IERR)
1697  call alloc_err( IERR, name, routine )
1698end if
1699END SUBROUTINE dealloc_r4
1700! ==================================================================
1701SUBROUTINE dealloc_d1( array, name, routine )
1702implicit none
1703real(DP), dimension(:),     pointer    :: array
1704character(len=*), optional, intent(in) :: name, routine
1705if (associated(array)) then
1706  call alloc_count( -size(array), 'D', name, routine )
1707  deallocate(array,stat=IERR)
1708  call alloc_err( IERR, name, routine )
1709end if
1710END SUBROUTINE dealloc_d1
1711! ==================================================================
1712SUBROUTINE dealloc_d2( array, name, routine )
1713implicit none
1714real(DP), dimension(:,:),   pointer    :: array
1715character(len=*), optional, intent(in) :: name, routine
1716if (associated(array)) then
1717  call alloc_count( -size(array), 'D', name, routine )
1718  deallocate(array,stat=IERR)
1719  call alloc_err( IERR, name, routine )
1720end if
1721END SUBROUTINE dealloc_d2
1722! ==================================================================
1723SUBROUTINE dealloc_d3( array, name, routine )
1724implicit none
1725real(DP), dimension(:,:,:), pointer    :: array
1726character(len=*), optional, intent(in) :: name, routine
1727if (associated(array)) then
1728  call alloc_count( -size(array), 'D', name, routine )
1729  deallocate(array,stat=IERR)
1730  call alloc_err( IERR, name, routine )
1731end if
1732END SUBROUTINE dealloc_d3
1733! ==================================================================
1734SUBROUTINE dealloc_d4( array, name, routine )
1735implicit none
1736real(DP), dimension(:,:,:,:), pointer  :: array
1737character(len=*), optional, intent(in) :: name, routine
1738if (associated(array)) then
1739  call alloc_count( -size(array), 'D', name, routine )
1740  deallocate(array,stat=IERR)
1741  call alloc_err( IERR, name, routine )
1742end if
1743END SUBROUTINE dealloc_d4
1744! ==================================================================
1745! COMPLEX versions
1746!
1747SUBROUTINE dealloc_c1( array, name, routine )
1748implicit none
1749complex(SP), dimension(:),   pointer   :: array
1750character(len=*), optional, intent(in) :: name, routine
1751if (associated(array)) then
1752  call alloc_count( -2*size(array), 'S', name, routine )
1753  deallocate(array,stat=IERR)
1754  call alloc_err( IERR, name, routine )
1755end if
1756END SUBROUTINE dealloc_c1
1757! ==================================================================
1758SUBROUTINE dealloc_c2( array, name, routine )
1759implicit none
1760complex(SP), dimension(:,:),  pointer  :: array
1761character(len=*), optional, intent(in) :: name, routine
1762if (associated(array)) then
1763  call alloc_count( -2*size(array), 'S', name, routine )
1764  deallocate(array,stat=IERR)
1765  call alloc_err( IERR, name, routine )
1766end if
1767END SUBROUTINE dealloc_c2
1768! ==================================================================
1769SUBROUTINE dealloc_z1( array, name, routine )
1770implicit none
1771complex(DP), dimension(:),   pointer   :: array
1772character(len=*), optional, intent(in) :: name, routine
1773if (associated(array)) then
1774  call alloc_count( -2*size(array), 'D', name, routine )
1775  deallocate(array,stat=IERR)
1776  call alloc_err( IERR, name, routine )
1777end if
1778END SUBROUTINE dealloc_z1
1779! ==================================================================
1780SUBROUTINE dealloc_z2( array, name, routine )
1781implicit none
1782complex(DP), dimension(:,:),  pointer  :: array
1783character(len=*), optional, intent(in) :: name, routine
1784if (associated(array)) then
1785  call alloc_count( -2*size(array), 'D', name, routine )
1786  deallocate(array,stat=IERR)
1787  call alloc_err( IERR, name, routine )
1788end if
1789END SUBROUTINE dealloc_z2
1790! ==================================================================
1791SUBROUTINE dealloc_z3( array, name, routine )
1792implicit none
1793complex(DP), dimension(:,:,:), pointer :: array
1794character(len=*), optional, intent(in) :: name, routine
1795if (associated(array)) then
1796  call alloc_count( -2*size(array), 'D', name, routine )
1797  deallocate(array,stat=IERR)
1798  call alloc_err( IERR, name, routine )
1799end if
1800END SUBROUTINE dealloc_z3
1801! ==================================================================
1802SUBROUTINE dealloc_z4( array, name, routine )
1803implicit none
1804complex(DP), dimension(:,:,:,:),  pointer  :: array
1805character(len=*), optional, intent(in) :: name, routine
1806if (associated(array)) then
1807  call alloc_count( -2*size(array), 'D', name, routine )
1808  deallocate(array,stat=IERR)
1809  call alloc_err( IERR, name, routine )
1810end if
1811END SUBROUTINE dealloc_z4
1812! ==================================================================
1813SUBROUTINE dealloc_l1( array, name, routine )
1814implicit none
1815logical, dimension(:),      pointer    :: array
1816character(len=*), optional, intent(in) :: name, routine
1817if (associated(array)) then
1818  call alloc_count( -size(array), 'L', name, routine )
1819  deallocate(array,stat=IERR)
1820  call alloc_err( IERR, name, routine )
1821end if
1822END SUBROUTINE dealloc_l1
1823! ==================================================================
1824SUBROUTINE dealloc_l2( array, name, routine )
1825implicit none
1826logical, dimension(:,:),    pointer    :: array
1827character(len=*), optional, intent(in) :: name, routine
1828if (associated(array)) then
1829  call alloc_count( -size(array), 'L', name, routine )
1830  deallocate(array,stat=IERR)
1831  call alloc_err( IERR, name, routine )
1832end if
1833END SUBROUTINE dealloc_l2
1834! ==================================================================
1835SUBROUTINE dealloc_l3( array, name, routine )
1836implicit none
1837logical, dimension(:,:,:),  pointer    :: array
1838character(len=*), optional, intent(in) :: name, routine
1839if (associated(array)) then
1840  call alloc_count( -size(array), 'L', name, routine )
1841  deallocate(array,stat=IERR)
1842  call alloc_err( IERR, name, routine )
1843end if
1844END SUBROUTINE dealloc_l3
1845! ==================================================================
1846SUBROUTINE dealloc_s1( array, name, routine )
1847implicit none
1848character(len=*), dimension(:), pointer :: array
1849character(len=*), optional, intent(in)  :: name, routine
1850if (associated(array)) then
1851  call alloc_count( -size(array)*len(array), 'S', name, routine )
1852  deallocate(array,stat=IERR)
1853  call alloc_err( IERR, name, routine )
1854end if
1855END SUBROUTINE dealloc_s1
1856
1857! ==================================================================
1858! Internal subroutines
1859! ==================================================================
1860
1861SUBROUTINE options( final_bounds, common_bounds, &
1862                    old_bounds, new_bounds, copy, shrink )
1863! Arguments
1864integer, dimension(:,:), intent(out) :: final_bounds
1865integer, dimension(:,:), intent(out) :: common_bounds
1866integer, dimension(:,:), intent(in)  :: old_bounds
1867integer, dimension(:,:), intent(in)  :: new_bounds
1868logical,       optional, intent(in)  :: copy
1869logical,       optional, intent(in)  :: shrink
1870
1871! Internal variables and arrays
1872logical want_shrink
1873
1874
1875!! AG*****
1876!  It might be worthwhile to check whether the user
1877!  atttemps to use bounds which do not make sense,
1878!  such as zero, or with upper<lower...
1879!!***
1880
1881! Find if it is a new allocation or a true reallocation,
1882! and if the contents need to be copied (saved)
1883if (ASSOCIATED_ARRAY) then
1884
1885  ! Check if array bounds have changed
1886  if ( all(new_bounds==old_bounds) ) then
1887    ! Old and new arrays are equal. Nothing needs to be done
1888    NEEDS_ALLOC   = .false.
1889    NEEDS_DEALLOC = .false.
1890    NEEDS_COPY    = .false.
1891  else
1892
1893    ! Want to shrink?
1894    if (present(shrink)) then
1895      want_shrink = shrink
1896    else
1897      want_shrink = DEFAULT%shrink
1898    end if
1899
1900    if (.not. want_shrink  &
1901        .and. all(new_bounds(1,:)>=old_bounds(1,:)) &
1902        .and. all(new_bounds(2,:)<=old_bounds(2,:)) ) then
1903      ! Old array is already fine. Nothing needs to be done
1904      NEEDS_ALLOC   = .false.
1905      NEEDS_DEALLOC = .false.
1906      NEEDS_COPY    = .false.
1907    else
1908      ! Old array needs to be substituted by a new array
1909      NEEDS_ALLOC   = .true.
1910      NEEDS_DEALLOC = .true.
1911      if (present(copy)) then
1912        NEEDS_COPY = copy
1913      else
1914        NEEDS_COPY = DEFAULT%copy
1915      end if
1916
1917      ! Ensure that bounds shrink only if desired
1918      if (want_shrink) then
1919        final_bounds(1,:) = new_bounds(1,:)
1920        final_bounds(2,:) = new_bounds(2,:)
1921      else
1922        final_bounds(1,:) = min( old_bounds(1,:), new_bounds(1,:) )
1923        final_bounds(2,:) = max( old_bounds(2,:), new_bounds(2,:) )
1924      end if
1925
1926      ! Find common section of old and new arrays
1927      common_bounds(1,:) = max( old_bounds(1,:), final_bounds(1,:) )
1928      common_bounds(2,:) = min( old_bounds(2,:), final_bounds(2,:) )
1929    end if
1930
1931  end if
1932
1933else
1934  ! Old array does not exist. Allocate new one
1935  NEEDS_ALLOC   = .true.
1936  NEEDS_DEALLOC = .false.
1937  NEEDS_COPY    = .false.
1938  final_bounds(1,:) = new_bounds(1,:)
1939  final_bounds(2,:) = new_bounds(2,:)
1940end if
1941
1942END SUBROUTINE options
1943
1944! ==================================================================
1945
1946SUBROUTINE alloc_count( delta_size, type, name, routine )
1947
1948implicit none
1949
1950integer, intent(in)          :: delta_size  ! +/-size(array)
1951character, intent(in)        :: type        ! 'I' => integer
1952                                            ! 'E' => integer*8
1953                                            ! 'R' => real*4
1954                                            ! 'D' => real*8
1955                                            ! 'L' => logical
1956                                            ! 'S' => character (string)
1957character(len=*), optional, intent(in) :: name
1958character(len=*), optional, intent(in) :: routine
1959
1960character(len=32)   :: aname, rname
1961character(len=1)    :: memType, task
1962real(DP)            :: delta_mem
1963logical             :: newPeak
1964logical,  save      :: header_written = .false.
1965logical,  save      :: tree_nullified = .false.
1966integer             :: memSize
1967
1968! Set routine name
1969if (present(routine)) then
1970  rname = routine
1971else
1972  rname = DEFAULT%routine
1973end if
1974
1975! Call siesta counting routine 'memory'
1976! Switched off in Aug.2009, and made 'memory' call alloc_count instead
1977! if (delta_size > 0) then
1978!   task = 'A'
1979! else
1980!   task = 'D'
1981! end if
1982! select case( type )
1983! case ('R')         ! Real --> single
1984!   memType = 'S'
1985!   memSize = abs(delta_size)
1986! case ('S')         ! Character (string) --> integer/4
1987!   memType = 'I'
1988!   memSize = abs(delta_size) / type_mem('I')
1989! case default
1990!   memType = type
1991!   memSize = abs(delta_size)
1992! end select
1993! call memory( task, memType, memSize, trim(rname) )
1994
1995
1996if (REPORT_LEVEL <= 0) return
1997
1998! Compound routine+array name
1999if (present(name) .and. present(routine)) then
2000  aname = trim(routine)//' '//name
2001else if (present(name) .and. DEFAULT%routine/=DEFAULT_ROUTINE) then
2002  aname = trim(DEFAULT%routine)//' '//name
2003else if (present(name)) then
2004  aname = name
2005else if (present(routine)) then
2006  aname = trim(routine)//' '//DEFAULT_NAME
2007else if (DEFAULT%routine/=DEFAULT_ROUTINE) then
2008  aname = trim(DEFAULT%routine)//' '//DEFAULT_NAME
2009else
2010  aname = DEFAULT_ROUTINE//' '//DEFAULT_NAME
2011end if
2012
2013MAX_LEN = max( MAX_LEN, len(trim(aname)) )
2014
2015! Find memory increment and total allocated memory
2016delta_mem = delta_size * type_mem(type)
2017TOT_MEM = TOT_MEM + delta_mem
2018if (TOT_MEM > PEAK_MEM+0.5_dp) then
2019  newPeak = .true.
2020  PEAK_MEM = TOT_MEM
2021  PEAK_ARRAY = aname
2022  PEAK_ROUTINE = rname
2023!  print'(/,a,f18.6),a,/)',
2024!    'alloc_count: Memory peak =', PEAK_MEM/MBYTE, ' Mbytes'
2025else
2026  newPeak = .false.
2027end if
2028
2029! Add/subtract/classify array memory
2030if (REPORT_LEVEL > 1) then
2031  if (.not.tree_nullified) then
2032    nullify(report_tree)
2033    tree_nullified = .true.
2034  end if
2035  call tree_add( report_tree, aname, delta_mem )
2036  if (newPeak) call tree_peak( report_tree )
2037end if
2038
2039! Print report, but only in node 0, as not all
2040! processors may follow the same route here
2041!   The detail/extent of the report increses with the value of level:
2042! level=0 : no report at all (the default)
2043! level=1 : only total memory peak and where it occurred
2044! level=2 : detailed report created but printed only upon request
2045! level=3 : detailed report printed at every new memory peak
2046! level=4 : print every individual reallocation or deallocation
2047
2048
2049if (newPeak .and. (REPORT_LEVEL==1 .or. REPORT_LEVEL==3) .and. &
2050    node == 0) then
2051  call print_report(.false.)
2052end if
2053
2054if (REPORT_LEVEL == 4 .and. node == 0) then
2055  if (.not.header_written) then
2056    write(REPORT_UNIT,'(/,a7,9x,1x,a4,28x,1x,2a15)') &
2057     'Routine', 'Name', 'Incr. (MB)', 'Total (MB)'
2058    header_written = .true.
2059  end if
2060  write(REPORT_UNIT,'(a16,1x,a32,1x,2f15.6)') &
2061     rname, aname, delta_mem/MBYTE, TOT_MEM/MBYTE
2062end if
2063END SUBROUTINE alloc_count
2064
2065! ==================================================================
2066
2067INTEGER FUNCTION type_mem( var_type )
2068!
2069! It is not clear that the sizes assumed are universal for
2070! non-Cray machines...
2071!
2072implicit none
2073character, intent(in) :: var_type
2074character(len=40)     :: message
2075
2076select case( var_type )
2077#ifdef OLD_CRAY
2078  case('I')
2079    type_mem = 8
2080  case('R')
2081    type_mem = 8
2082  case('L')
2083    type_mem = 8
2084#else
2085  case('I')
2086    type_mem = 4
2087  case('R')
2088    type_mem = 4
2089  case('L')
2090    type_mem = 4
2091#endif
2092case('E')
2093  type_mem = 8
2094case('D')
2095  type_mem = 8
2096case('S')
2097  type_mem = 1
2098case default
2099  write(message,"(2a)") &
2100    'alloc_count: ERROR: unknown type = ', var_type
2101  call die(trim(message))
2102end select
2103
2104END FUNCTION type_mem
2105
2106! ==================================================================
2107
2108RECURSIVE SUBROUTINE tree_add( t, name, delta_mem )
2109
2110implicit none
2111type(TREE),       pointer    :: t
2112character(len=*), intent(in) :: name
2113real(DP),         intent(in) :: delta_mem
2114
2115logical, save :: warn_negative = .true.
2116
2117if (.not.associated(t)) then
2118  allocate( t )
2119  t%name = name
2120  t%mem  = delta_mem
2121  t%max  = delta_mem
2122  t%peak = 0._dp
2123  nullify( t%left, t%right )
2124else if (name == t%name) then
2125  t%mem = t%mem + delta_mem
2126  ! The abs is to handle the case of apparent de_allocs without re_allocs,
2127  ! caused by routine/name argument mismatches
2128  if (abs(t%mem) > abs(t%max)) t%max = t%mem
2129else if ( llt(name,t%name) ) then
2130  call tree_add( t%left, name, delta_mem )
2131else
2132  call tree_add( t%right, name, delta_mem )
2133end if
2134
2135if (warn_negative .and. t%mem<0._dp) then
2136  call parallel_init()   ! Make sure that node and Nodes are initialized
2137  if (Node==0) then
2138   write(6,'(/,a,/,2a,/,a,f18.0,a)')  &
2139      'WARNING: alloc-realloc-dealloc name mismatch',  &
2140      '         Name: ', trim(name),                   &
2141      '         Size: ', t%mem, ' Bytes'
2142    if (Nodes>1) write(6,'(9x,a,i6)') 'Node:', Node
2143    write(6,'(9x,a)') 'Subsequent mismatches will not be reported'
2144    warn_negative = .false.  ! Print this warning only once
2145  end if
2146end if
2147
2148END SUBROUTINE tree_add
2149
2150! ==================================================================
2151
2152RECURSIVE SUBROUTINE tree_peak( t )
2153
2154implicit none
2155type(TREE), pointer :: t
2156
2157if (.not.associated(t)) return
2158
2159t%peak = t%mem
2160call tree_peak( t%left )
2161call tree_peak( t%right )
2162
2163END SUBROUTINE tree_peak
2164
2165! ==================================================================
2166
2167RECURSIVE SUBROUTINE tree_print( t )
2168
2169implicit none
2170type(TREE), pointer :: t
2171
2172if (.not.associated(t)) return
2173
2174call tree_print( t%left )
2175
2176if (abs(t%max) >= REPORT_THRESHOLD) then
2177  write(REPORT_UNIT,'(a,1x,3f15.6,f9.2)') &
2178    t%name(1:MAX_LEN), t%mem/MBYTE, t%max/MBYTE, t%peak/MBYTE, &
2179    100._dp * t%peak / (PEAK_MEM + tiny(PEAK_MEM) )
2180end if
2181
2182call tree_print( t%right )
2183
2184END SUBROUTINE tree_print
2185
2186! ==================================================================
2187
2188SUBROUTINE print_report(all)
2189
2190implicit none
2191
2192! Whether MPI reductions should be performed
2193! If, not then ensure that no MPI calls are performed
2194logical, intent(in), optional :: all
2195
2196character(len=80)   :: string = 'Name'
2197character           :: date*8, time*10, zone*5
2198integer             :: iNode, peakNode
2199real(dp)            :: maxPeak
2200real(dp),allocatable:: nodeMem(:), nodePeak(:)
2201logical             :: lall
2202
2203#ifdef MPI
2204integer           :: MPIerror
2205#endif
2206
2207! Enables parallel call of print-report
2208lall = .true.
2209if ( present(all) ) lall = all
2210
2211! Only if MPI should all be used
2212if ( lall ) then
2213! Make sure that variables node and Nodes are initialized
2214call parallel_init()
2215end if
2216
2217! Allocate and initialize two small arrays
2218allocate( nodeMem(0:Nodes-1), nodePeak(0:Nodes-1) )
2219! initialize (in case all == .false.)
2220nodeMem = 0._dp
2221nodePeak = 0._dp
2222
2223! Initializations for Nodes=1 (serial case)
2224nodeMem(node) = TOT_MEM
2225nodePeak(node) = PEAK_MEM
2226peakNode = node
2227
2228! In parallel, find the memory values of all nodes
2229#ifdef MPI
2230if (Nodes > 1 .and. lall ) then
2231  ! Gather the present and peak memories of all nodes
2232  call MPI_AllGather( TOT_MEM, 1, MPI_double_precision, &
2233                      nodeMem, 1, MPI_double_precision, &
2234                      MPI_COMM_WORLD, MPIerror )
2235  call MPI_AllGather( PEAK_MEM, 1, MPI_double_precision, &
2236                      nodePeak, 1, MPI_double_precision, &
2237                      MPI_COMM_WORLD, MPIerror )
2238  ! Find the node with the highest peak of memory
2239  maxPeak = 0
2240  do iNode = 0,Nodes-1
2241    if (nodePeak(iNode) > maxPeak) then
2242      peakNode = iNode
2243      maxPeak = nodePeak(iNode)
2244    end if
2245  end do ! iNode
2246  ! Change the writing node for the peak-node information
2247  if (node==0 .and. peakNode/=0) close( REPORT_UNIT )
2248  call MPI_Barrier( MPI_COMM_WORLD, MPIerror )
2249  if (node==peakNode .and. peakNode/=0) then
2250     call io_assign(REPORT_UNIT)
2251     open( unit=REPORT_UNIT, file=REPORT_FILE, &
2252          status='unknown', position='append' )
2253  end if
2254end if ! (Nodes>1)
2255#endif
2256
2257! The report is printed by the highest-peak node
2258if (node == peakNode) then
2259
2260  ! AG: Commented out to allow multiple batches of information
2261  ! if (REPORT_LEVEL < 4) rewind(REPORT_UNIT)
2262
2263  call date_and_time( date, time, zone )
2264
2265  write(REPORT_UNIT,'(/,a,16a)')                &
2266    'Allocation summary at ',                   &
2267    date(1:4),'/',date(5:6),'/',date(7:8),' ',  &
2268    time(1:2),':',time(3:4),':',time(5:10),' ', &
2269    zone(1:3),':',zone(4:5)
2270
2271  if (Nodes > 1) then
2272    write(REPORT_UNIT,'(/,(a,f18.6,a))')            &
2273      'Present memory all nodes : ', sum(nodeMem)/MBYTE,  ' MB', &
2274      'Added peak mem all nodes : ', sum(nodePeak)/MBYTE, ' MB', &
2275      'Min peak memory in a node: ', minval(nodePeak)/MBYTE, ' MB', &
2276      'Max peak memory in a node: ', maxval(nodePeak)/MBYTE, ' MB'
2277!   Impractical for many nodes:
2278!    write(REPORT_UNIT,'(/,a,/,(i6,f12.6))') &
2279!      'Memory peaks of individual nodes (Mb):', &
2280!      (iNode,nodePeak(iNode)/MBYTE,iNode=0,Nodes-1)
2281    write(REPORT_UNIT,'(/,a,i6)') &
2282      'Maximum peak of memory occurred in node:', peakNode
2283  end if
2284
2285  write(REPORT_UNIT,'(2(/,a,f18.6,a),/,2a,/,2a)')            &
2286    'Present memory allocation: ', TOT_MEM/MBYTE,  ' MB', &
2287    'Maximum memory allocation: ', PEAK_MEM/MBYTE, ' MB', &
2288    'Occurred after allocating: ', trim(PEAK_ARRAY),         &
2289    'In routine:                ', trim(PEAK_ROUTINE)
2290
2291  if (REPORT_LEVEL > 1) then
2292    if (REPORT_THRESHOLD > 0._dp) then
2293      write(REPORT_UNIT,'(/,a,f12.6,a,/,a,1x,3a15,a9)') &
2294        'Allocated sizes (in MByte) of arrays larger than ', &
2295        REPORT_THRESHOLD/MBYTE, ' MB:', &
2296        string(1:MAX_LEN), 'Present', 'Maximum', 'At peak', '%'
2297    else
2298      write(REPORT_UNIT,'(/,a,/,a,1x,3a15,a9)') &
2299        'Allocated array sizes (in MByte):', &
2300        string(1:MAX_LEN), 'Present', 'Maximum', 'At peak', '%'
2301    end if
2302    call tree_print( report_tree )
2303  end if
2304
2305  ! Close file if not common IO node
2306  if ( node /= 0 ) call io_close(REPORT_UNIT)
2307end if ! (node == peakNode)
2308
2309! Change again the writing node for the rest of the report
2310#ifdef MPI
2311if (lall) call MPI_Barrier( MPI_COMM_WORLD, MPIerror )
2312if (node==0 .and. peakNode/=0) &
2313     open( unit=REPORT_UNIT, file=REPORT_FILE, &
2314     status='unknown', position='append' )
2315#endif
2316
2317deallocate( nodeMem, nodePeak )
2318
2319END SUBROUTINE print_report
2320
2321! ==================================================================
2322
2323SUBROUTINE alloc_err( ierr, name, routine, bounds )
2324
2325implicit none
2326
2327integer,                    intent(in) :: ierr
2328character(len=*), optional, intent(in) :: name
2329character(len=*), optional, intent(in) :: routine
2330integer, dimension(:,:), optional, intent(in) :: bounds
2331
2332integer i
2333
2334if (ierr/=0) then
2335  if (ionode) print*, 'alloc_err: allocate status error', ierr
2336  if (present(name).and.present(routine)) then
2337    if (ionode) print*, 'alloc_err: array ', name, &
2338               ' requested by ', routine
2339  elseif (present(name)) then
2340    if (ionode) print*, 'alloc_err: array ', name, &
2341               ' requested by unknown'
2342  elseif (present(routine)) then
2343    if (ionode) print* , 'alloc_err: array unknown', &
2344               ' requested by ', routine
2345  endif
2346  if (ionode.and.present(bounds))                            &
2347    print '(a,i3,2i10)', ('alloc_err: dim, lbound, ubound:',  &
2348          i,bounds(1,i),bounds(2,i),                         &
2349          i=1,size(bounds,dim=2))
2350
2351  call die('alloc_err: allocate error')
2352end if
2353
2354END SUBROUTINE alloc_err
2355
2356! ==================================================================
2357
2358END MODULE alloc
2359