Discussion of coding approach/style for dynamic memory

As a basis for the discussion, here's how I've currently coded NEMO (v.3.2) to use dynamic memory. I've used ''...'' to indicate places where I've missed-out chunks of code for clarity/brevity.

Highest-level changes are in OPA_SRC/opa.F90 and the routine opa_init():

#if defined key_mpp_dyndist
      ! Decide on size of grid now that we have our communicator size
      ! If we're not using dynamic memory then mpp_partition does nothing.

#if   defined key_mpp_mpi   ||   defined key_mpp_shmem
      CALL opa_partition()
#else
      jpni = 1
      jpnj = 1
      jpnij = jpni*jpnj
#endif
      ! Calculate domain dimensions given calculated jpi and jpj
      jpi = ( jpiglo-2*jpreci + (jpni-1) ) / jpni + 2*jpreci !: first  dim.
      jpj = ( jpjglo-2*jprecj + (jpnj-1) ) / jpnj + 2*jprecj !: second dim.
      jpim1 = jpi-1                                          !: inner domain indices
      jpjm1 = jpj-1                                          !:   "           "
      jpkm1 = jpk-1                                          !:   "           "
      jpij  = jpi*jpj                                        !:  jpi x j

      ! Now we know the dimensions of the grid, allocate arrays
      CALL opa_alloc()
#endif

The inclusion of the new code is currently controlled by the cpp key 'key_mpp_dyndist' but it would greatly improve the cleanliness of the code if we make a complete break from the static-memory version and thus can drop the use of this cpp key. My understanding of the conclusion of the meeting in Southampton was that this is what we're going to do.

This addition to opa_init() calls two further new routines, opa_partition() and opa_alloc():

   SUBROUTINE opa_partition
     USE par_oce
     IMPLICIT none
     INTEGER, PARAMETER :: nfactmax = 20
     INTEGER :: nfact ! The no. of factors returned
     INTEGER :: ierr  ! Error flag
     INTEGER :: i
     INTEGER :: idiff, mindiff, imin ! For choosing pair of factors that are
                                     ! closest in value
     INTEGER, DIMENSION(nfactmax) :: ifact ! Array of factors
     ierr = 0

#if ! defined key_mpp_dyndist
     RETURN ! If we aren't using dynamic memory then jpnj and jpni are set
            ! at compile time
#else

     CALL factorise(ifact, nfactmax, nfact, mppsize, ierr)

     IF(nfact <= 1)THEN
        WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
        WRITE (numout, *) '       : using grid of ',mppsize,' x 1'
        jpnj = 1
        jpni = mppsize
     ELSE
        ! Search through factors for the pair that are closest in value
        mindiff = 1000000
        imin    = 1
        DO i=1,nfact-1,2
           idiff = ABS(ifact(i) - ifact(i+1))
           IF(idiff < mindiff)THEN
              mindiff = idiff
              imin = i
           END IF
        END DO
        jpnj = ifact(imin)
        jpni = ifact(imin + 1)
     ENDIF
     jpnij = jpni*jpnj

     WRITE(*,*) 'ARPDBG: jpni = ',jpni,'jpnj = ',jpnj,'jpnij = ',jpnij

#endif

   END SUBROUTINE opa_partition

where factorise() returns the prime factors of its argument.

opa_alloc() oversees the allocation of all of the module-level arrays:

   SUBROUTINE opa_alloc
#if defined key_mpp_dyndist

#if defined key_lim2
     USE limwri_2,   ONLY: lim_wri_2_malloc
     USE limdmp_2,   ONLY: lim_dmp_2_malloc
     USE limhdf_2,   ONLY: lim_hdf_2_malloc
     USE ice_2,      ONLY: ice_2_malloc
     USE limsbc_2,   ONLY: lim_sbc_2_malloc
     USE thd_ice_2,  ONLY: thd_ice_2_malloc
     USE limdia_2,   ONLY: lim_dia_2_malloc
     USE dom_ice_2,  ONLY: dom_ice_2_malloc
#endif
     USE ldfslp,     ONLY: ldf_slp_malloc
     USE ldftra_oce, ONLY: ldftra_oce_malloc
     USE ldfdyn_oce, ONLY: ldfdyn_oce_malloc
#if defined key_vvl
     USE domvvl,     ONLY: dom_vvl_malloc
#endif
     USE dom_oce,    ONLY: dom_oce_malloc
...
...
     IMPLICIT none
     INTEGER, PARAMETER           :: NUMCALLS = 58
     INTEGER, DIMENSION(NUMCALLS) :: ierr
     INTEGER                      :: i, icount

     ierr = 0
     icount = 1

#if defined key_lim2
     ierr = ierr + lim_wri_2_malloc()
     ierr = ierr + lim_dmp_2_malloc()
     ierr = ierr + lim_hdf_2_malloc()
     ierr = ierr + ice_2_malloc()
     ierr = ierr + lim_sbc_2_malloc()
     ierr = ierr + thd_ice_2_malloc()
     ierr = ierr + lim_dia_2_malloc()
     ierr = ierr + dom_ice_2_malloc()
#endif
     ierr = ierr + ldf_slp_malloc()
     ierr = ierr + ldftra_oce_malloc()
     ierr = ierr + ldfdyn_oce_malloc()
#if defined key_vvl
     ierr = ierr + dom_vvl_malloc()
#endif
     ierr = ierr + dom_oce_malloc()
...
...
! Should do an MPI_SUM on ierr and then everyone can MPI_FINALIZE() if the
! memory allocation has failed on any one PE
     IF(ierr > 0)THEN
        WRITE(numout,*) 
        WRITE(numout,*) 'ERROR: Allocation of memory failed in opa_alloc'
        STOP 'ERROR: Allocation of memory failed in opa_alloc'
     END IF
#else
     RETURN ! Not using dynamic memory therefore do nothing
#endif

   END SUBROUTINE opa_alloc

Each of the modules must then be edited and the declaration of arrays with dimensions that depend upon jpi or jpj be changed. e.g. for OPA_SRC/DOM/dom_oce.F90:

MODULE dom_oce
   !!======================================================================
   !!                       ***  MODULE dom_oce  ***
   !!       
   !! ** Purpose :   Define in memory all the ocean space domain variables
   !!======================================================================
   !! History :  1.0  ! 2005-10  (A. Beckmann, G. Madec)  reactivate s-coordinate 
   !!----------------------------------------------------------------------
   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 
   !! $Id: dom_oce.F90 1601 2009-08-11 10:09:19Z ctlod $ 
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------
   USE par_oce      ! ocean parameters

   IMPLICIT NONE
...
...
#if defined key_mpp_dyndist
   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: & !: 
#else
   REAL(wp), PUBLIC,              DIMENSION(jpi,jpj,jpk) :: & !:
#endif
      tmask, umask,    &  !: land/ocean mask at T-, U-, V- and F-points
      vmask, fmask        !:

   REAL(wp), PUBLIC, DIMENSION(jpiglo) ::   tpol, fpol          !: north fold mask (nperio= 3 or 4)

#if defined key_noslip_accurate
   INTEGER, PUBLIC,                    DIMENSION(4,jpk) ::   npcoa   !: ???
#if defined key_mpp_dyndist
   INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nicoa, njcoa        !: ???
#else
   INTEGER, PUBLIC,        DIMENSION(2*(jpi+jpj),4,jpk) ::   nicoa, njcoa        !: ???
#endif
                                         
#endif

   !!----------------------------------------------------------------------
   !! agrif domain
   !!----------------------------------------------------------------------
#if defined key_agrif
   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .TRUE.    !: agrif flag
#else
   LOGICAL, PUBLIC, PARAMETER ::   lk_agrif = .FALSE.   !: agrif flag
#endif


#if defined key_mpp_dyndist
   PUBLIC dom_oce_malloc
#endif

#if defined key_mpp_dyndist

CONTAINS

  FUNCTION dom_oce_malloc()
    USE par_oce, Only: jpi, jpj, jpk, jpnij
    IMPLICIT none
    INTEGER :: dom_oce_malloc
    INTEGER :: ierr, ierr_tmp
    
    ierr = 0
...
...
    ALLOCATE(tmask(jpi,jpj,jpk), umask(jpi,jpj,jpk),    & 
             vmask(jpi,jpj,jpk), fmask(jpi,jpj,jpk), Stat=ierr_tmp)
    ierr = ierr + ierr_tmp

#if defined key_noslip_accurate
    ALLOCATE(nicoa(2*(jpi+jpj),4,jpk), njcoa(2*(jpi+jpj),4,jpk), &
             Stat=ierr_tmp)
    ierr = ierr + ierr_tmp
#endif

    IF(ierr != 0)THEN
      WRITE() "Some error message detailing which module's ..._malloc() routine failed"
    END IF

    ! Error status passed back up
    dom_oce_malloc = ierr

  END FUNCTION dom_oce_malloc
#endif

Finally (I think), there's the issue of work-space or local arrays. These are arrays local to subroutines that are declared with their dimensions in terms of jpi and jpj. Since jpi and jpj are no longer parameters, these arrays are now 'automatic' -- that is, they are allocated on the stack when the program enters the subroutine and free'd again when it leaves. Running out of stack memory will cause the program to abort abruptly and therefore it is more controlled to use allocate for big workspace arrays. (Although the amount of available stack can be set to 'unlimited' by using the shell's ulimit command, what limit this actually results in will be OS/set-up dependent.) I suggest therefore that we edit the code such that any large (3d or higher?) workspace arrays are made explicitly allocatable and then allocated the first time the subroutine is called. If we declare them with the SAVE attribute then the allocated memory will remain available to us in future calls, e.g.:

   SUBROUTINE dia_ptr( kt )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE dia_ptr  ***
      !!----------------------------------------------------------------------
      IMPLICIT none
      INTEGER, INTENT(in) ::   kt   ! ocean time step index
      !!
      INTEGER  ::   jk, jj, ji   ! dummy loop
      REAL(wp) ::   zsverdrup    ! conversion from m3/s to Sverdrup
      REAL(wp) ::   zpwatt       ! conversion from W    to PW
      REAL(wp) ::   zggram       ! conversion from g    to Pg
#if defined key_mpp_dyndist
      REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vt, vs   ! 3D workspace
#else
      REAL(wp),              DIMENSION(jpi,jpj,jpk) :: vt, vs   ! 3D workspace
#endif
      !!----------------------------------------------------------------------

      IF( kt == nit000 .OR. MOD( kt, nf_ptr ) == 0 )   THEN

#if defined key_mpp_dyndist
         IF(.not. ALLOCATED(vt))THEN
            ALLOCATE(vt(jpi,jpj,jpk), vs(jpi,jpj,jpk), Stat=ji)
            IF(ji .ne. 0)THEN
               ! Need an interface to MPI_Abort() here really
               STOP 'ERROR: Allocation of vt/vs failed in dia_ptr'
            END IF
         END IF
#endif
...
...

If an allocate() is done in this way then handling any failure is difficult since there's no guarantee that the allocate() call will have failed for all MPI processes. Therefore, I think the best bet is to print out an error message and then do an MPI_Abort on MPI_COMM_WORLD.