Index: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dommsk.F90
===================================================================
--- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 2647)
+++ branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dommsk.F90	(revision 2648)
@@ -12,5 +12,5 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O manager
+   USE lib_mpp
 
    IMPLICIT NONE
@@ -19,10 +19,5 @@
    PUBLIC   dom_msk    ! routine called by inidom.F90
 
-#if defined key_degrad
-   !! ------------------------------------------------
-   !! Degradation method
-   !! --------------------------------------------------
-   REAL(wp), PUBLIC, DIMENSION (jpi,jpj,jpk) ::   facvol  !! volume for degraded regions
-#endif
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   facvol  !! volume for degraded regions
 
    !! * Substitutions
@@ -34,4 +29,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
    
    SUBROUTINE dom_msk
@@ -49,9 +45,16 @@
       !!               tpol     : ???
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: iwrk_in_use, iwrk_not_released
+      USE wrk_nemo, ONLY: imsk => iwrk_2d_1
       INTEGER  ::   ji, jk                   ! dummy loop indices
       INTEGER  ::   iif, iil, ijf, ijl       ! local integers
-      INTEGER, DIMENSION(jpi,jpj) ::  imsk   ! 2D workspace
       !!---------------------------------------------------------------------
       !
+      IF( iwrk_in_use(2, 1) ) THEN
+         CALL ctl_stop('dom_msk: ERROR: requested workspace arrays unavailable.')  ;  RETURN
+      END IF
+      !
+      CALL dom_msk_alloc
+
       ! Interior domain mask (used for global sum)
       ! --------------------
@@ -95,5 +98,21 @@
       ENDIF
       !
+      IF( iwrk_not_released(2, 1) ) CALL ctl_stop('dom_msk: failed to release workspace arrays.')
+      !
    END SUBROUTINE dom_msk
+
+   SUBROUTINE dom_msk_alloc()
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dom_msk_alloc  ***
+      !!---------------------------------------------------------------------
+#if defined key_degrad
+      INTEGER :: ierr 
+
+      ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr )
+      IF( ierr /= 0 )  &
+        &           CALL ctl_stop('STOP', 'dom_msk : unable to allocate facvol array')
+#endif
+
+   END SUBROUTINE dom_msk_alloc
 
    !!======================================================================
Index: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/domrea.F90
===================================================================
--- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 2647)
+++ branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/domrea.F90	(revision 2648)
@@ -16,5 +16,6 @@
    USE dommsk          ! domain: masks
    USE lbclnk          ! lateral boundary condition - MPP exchanges
-   USE in_out_manager  ! I/O manager
+   USE lib_mpp 
+   USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
 
    IMPLICIT NONE
@@ -53,9 +54,9 @@
       !!----------------------------------------------------------------------
       USE iom
+      USE wrk_nemo, ONLY: zmbk => wrk_2d_1, zprt => wrk_2d_2, zprw => wrk_2d_3
       !!
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
       INTEGER  ::   ik, inum0 , inum1 , inum2 , inum3 , inum4   ! local integers
       REAL(wp) ::   zrefdep         ! local real
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk, zprt, zprw   ! 2D workspace
       !!----------------------------------------------------------------------
 
@@ -63,4 +64,8 @@
       IF(lwp) WRITE(numout,*) 'dom_rea : read NetCDF mesh and mask information file(s)'
       IF(lwp) WRITE(numout,*) '~~~~~~~'
+
+      IF( wrk_in_use(2, 1,2,3)  ) THEN
+         CALL ctl_stop('dom_rea: ERROR: requested workspace arrays unavailable.') ; RETURN
+      END IF
 
       zmbk(:,:) = 0._wp
@@ -141,4 +146,6 @@
          CALL iom_get( inum3, jpdom_data, 'e2u', e2u )
          CALL iom_get( inum3, jpdom_data, 'e2v', e2v )
+
+         e1e2t(:,:) = e1t(:,:) * e2t(:,:)                              ! surface at T-points
 
          CALL iom_get( inum3, jpdom_data, 'ff', ff )
@@ -314,4 +321,6 @@
       END SELECT
       !
+      IF( wrk_not_released(2, 1,2,3)  ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
+      !
    END SUBROUTINE dom_rea
 
@@ -327,11 +336,17 @@
       !! ** Action  : - update mbathy: level bathymetry (in level index)
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: zmbk => wrk_2d_4
+      !
       INTEGER ::   ji, jj   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace 
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+
       !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) '    zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
       IF(lwp) WRITE(numout,*) '    ~~~~~~~~~~~~~'
+      !
+      IF( wrk_in_use(2, 4) ) THEN
+         CALL ctl_stop('dom_rea: ERROR: requested workspace arrays unavailable.')  ;  RETURN
+      END IF
       !
       mbkt(:,:) = MAX( mbathy(:,:) , 1 )    ! bottom k-index of T-level (=1 over land)
@@ -347,4 +362,6 @@
       zmbk(:,:) = REAL( mbkv(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
       !
+      IF( wrk_not_released(2, 4) ) CALL ctl_stop('dom_rea:failed to release workspace arrays.')
+      !
    END SUBROUTINE zgr_bot_level
 
Index: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90
===================================================================
--- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 2647)
+++ branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/dtadyn.F90	(revision 2648)
@@ -63,32 +63,32 @@
    INTEGER ::   numfl_t, numfl_u, numfl_v, numfl_w
 
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: tdta       ! temperature at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: sdta       ! salinity at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: udta       ! zonal velocity at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: vdta       ! meridional velocity at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wdta       ! vertical velocity at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: avtdta     ! vertical diffusivity coefficient
-
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: hmlddta    ! mixed layer depth at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: wspddta    ! wind speed at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: frlddta    ! sea-ice fraction at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: empdta     ! E-P at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: qsrdta     ! short wave heat flux at two consecutive times
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: bblxdta    ! frequency of bbl in the x direction at 2 consecutive times 
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: bblydta    ! frequency of bbl in the y direction at 2 consecutive times 
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: tdta       ! temperature at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: sdta       ! salinity at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: udta       ! zonal velocity at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vdta       ! meridional velocity at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wdta       ! vertical velocity at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: avtdta     ! vertical diffusivity coefficient
+
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: hmlddta    ! mixed layer depth at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: wspddta    ! wind speed at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: frlddta    ! sea-ice fraction at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: empdta     ! E-P at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: qsrdta     ! short wave heat flux at two consecutive times
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: bblxdta    ! frequency of bbl in the x direction at 2 consecutive times 
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: bblydta    ! frequency of bbl in the y direction at 2 consecutive times 
    LOGICAL :: l_offbbl
 #if defined key_ldfslp
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: uslpdta    ! zonal isopycnal slopes
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: vslpdta    ! meridional isopycnal slopes
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wslpidta   ! zonal diapycnal slopes
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: wslpjdta   ! meridional diapycnal slopes
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: uslpdta    ! zonal isopycnal slopes
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: vslpdta    ! meridional isopycnal slopes
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpidta   ! zonal diapycnal slopes
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: wslpjdta   ! meridional diapycnal slopes
 #endif
 #if ! defined key_degrad &&  defined key_traldf_c2d && defined key_traldf_eiv 
-   REAL(wp), DIMENSION(jpi,jpj    ,2) :: aeiwdta    ! G&M coefficient
+   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:,:) :: aeiwdta    ! G&M coefficient
 #endif
 #if defined key_degrad
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: ahtudta, ahtvdta, ahtwdta   ! Lateral diffusivity
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: ahtudta, ahtvdta, ahtwdta   ! Lateral diffusivity
 # if defined key_traldf_eiv
-   REAL(wp), DIMENSION(jpi,jpj,jpk,2) :: aeiudta, aeivdta, aeiwdta   ! G&M coefficient
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: aeiudta, aeivdta, aeiwdta   ! G&M coefficient
 # endif
 #endif
@@ -296,4 +296,35 @@
    END SUBROUTINE dta_dyn
 
+   INTEGER FUNCTION dta_dyn_alloc()
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE dta_dyn_alloc  ***
+      !!---------------------------------------------------------------------
+
+      ALLOCATE( tdta    (jpi,jpj,jpk,2), sdta    (jpi,jpj,jpk,2),    &
+       &        udta    (jpi,jpj,jpk,2), vdta    (jpi,jpj,jpk,2),    &
+       &        wdta    (jpi,jpj,jpk,2), avtdta  (jpi,jpj,jpk,2),    &
+#if defined key_ldfslp 
+       &        uslpdta (jpi,jpj,jpk,2), vslpdta (jpi,jpj,jpk,2),    &
+       &        wslpidta(jpi,jpj,jpk,2), wslpjdta(jpi,jpj,jpk,2),    &
+#endif
+#if defined key_degrad
+       &        ahtudta (jpi,jpj,jpk,2), ahtvdta (jpi,jpj,jpk,2),    &
+       &        ahtwdta (jpi,jpj,jpk,2),                             &
+# if defined key_traldf_eiv
+       &        aeiudta (jpi,jpj,jpk,2), aeivdta (jpi,jpj,jpk,2),    &
+       &        aeiwdta (jpi,jpj,jpk,2),                             &
+# endif
+#endif
+#if ! defined key_degrad &&  defined key_traldf_c2d && defined key_traldf_eiv
+       &        aeiwdta (jpi,jpj,    2),                             &
+#endif
+
+       &        hmlddta (jpi,jpj,    2), wspddta (jpi,jpj,    2),    &
+       &        frlddta (jpi,jpj,    2), qsrdta  (jpi,jpj,    2),    &
+       &        empdta  (jpi,jpj,    2),                         STAT=dta_dyn_alloc ) 
+
+      IF( dta_dyn_alloc /= 0 ) CALL ctl_warn('dta_dyn_alloc: failed to allocate facvol array.')
+
+   END FUNCTION dta_dyn_alloc
 
    SUBROUTINE dynrea( kt, kenr )
@@ -305,23 +336,27 @@
       !! ** Method : READ the kenr records of DATA and store in udta(...,2), ....  
       !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY: zu    => wrk_3d_1 , zv    => wrk_3d_2 , zw    => wrk_3d_3
+      USE wrk_nemo, ONLY: zt    => wrk_3d_4 , zs    => wrk_3d_5
+      USE wrk_nemo, ONLY: zavt  => wrk_3d_6 , zhdiv => wrk_3d_7
+      USE wrk_nemo, ONLY: zahtu => wrk_3d_8 , zahtv => wrk_3d_9 , zahtw => wrk_3d_10
+      USE wrk_nemo, ONLY: zaeiu => wrk_3d_11, zaeiv => wrk_3d_12, zaeiw => wrk_3d_13
+      !
+      USE wrk_nemo, ONLY: zemp  => wrk_2d_1 , zqsr  => wrk_2d_2 , zmld  => wrk_2d_3
+      USE wrk_nemo, ONLY: zice  => wrk_2d_4 , zwspd => wrk_2d_5 
+      USE wrk_nemo, ONLY: ztaux => wrk_2d_6 , ztauy => wrk_2d_7
+      USE wrk_nemo, ONLY: zbblx => wrk_2d_8 , zbbly => wrk_2d_9
+      USE wrk_nemo, ONLY: zaeiw2d => wrk_2d_10
+      !
       INTEGER, INTENT(in) ::   kt, kenr   ! time index
       !!
       INTEGER ::  jkenr
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::  zu, zv, zw, zt, zs, zavt , zhdiv              ! 3D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zemp, zqsr, zmld, zice, zwspd, ztaux, ztauy   ! 2D workspace
-      REAL(wp), DIMENSION(jpi,jpj)     ::  zbblx, zbbly
-
-#if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
-      REAL(wp), DIMENSION(jpi,jpj) :: zaeiw 
-#endif
-#if defined key_degrad
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zahtu, zahtv, zahtw  !  Lateral diffusivity
-# if defined key_traldf_eiv
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zaeiu, zaeiv, zaeiw  ! G&M coefficient
-# endif
-#endif
-      !!----------------------------------------------------------------------
-
-      ! 0. Initialization
+      !!----------------------------------------------------------------------
+
+      ! 0. Memory allocation
+      IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13) .OR. &
+          wrk_in_use(2, 1,2,3,4,5,6,7,8,9,10)               ) THEN
+         CALL ctl_stop('domrea/dta_dyn: requested workspace arrays unavailable')  ;  RETURN
+      END IF
       
       ! cas d'un fichier non periodique : on utilise deux fois le premier et
@@ -390,5 +425,5 @@
 
 #if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv 
-      CALL iom_get( numfl_w, jpdom_data, 'soleaeiw', zaeiw (:,: ), jkenr )
+      CALL iom_get( numfl_w, jpdom_data, 'soleaeiw', zaeiw2d(:,: ), jkenr )
 #endif
 
@@ -413,5 +448,5 @@
 
 #if ! defined key_degrad && defined key_traldf_c2d && defined key_traldf_eiv
-      aeiwdta(:,:,2)  = zaeiw(:,:) * tmask(:,:,1)
+      aeiwdta(:,:,2)  = zaeiw2d(:,:) * tmask(:,:,1)
 #endif
 
@@ -451,4 +486,9 @@
       ENDIF
       !      
+      IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13) .OR. &
+          wrk_not_released(2, 1,2,3,4,5,6,7,8,9,10)               ) THEN
+         CALL ctl_stop('domrea/dta_dyn: failed to release workspace arrays.')
+      END IF
+      !
    END SUBROUTINE dynrea
 
@@ -462,9 +502,13 @@
       !! ** Method :
       !!----------------------------------------------------------------------
-      REAL(wp) ::   znspyr   !: number of time step per year
+      REAL(wp) :: znspyr   !: number of time step per year
+      INTEGER  :: ierr
       !!
       NAMELIST/namdyn/ ndtadyn, ndtatot, nsptint, lperdyn,  &
       &                cfile_grid_T, cfile_grid_U, cfile_grid_V, cfile_grid_W
       !!----------------------------------------------------------------------
+
+      ierr = dta_dyn_alloc()
+      IF( ierr /= 0 )  CALL ctl_stop( 'STOP', 'dta_dyn_alloc : unable to allocate standard ocean arrays' )
 
       !  Define the dynamical input parameters
Index: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90
===================================================================
--- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 2647)
+++ branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/nemogcm.F90	(revision 2648)
@@ -5,4 +5,5 @@
    !!======================================================================
    !! History :  3.3  ! 2010-05  (C. Ethe)  Full reorganization of the off-line: phasing with the on-line
+   !!            4.0  ! 2011-01  (C. Ethe, A. R. Porter, STFC Daresbury) dynamical allocation
    !!----------------------------------------------------------------------
 
@@ -26,4 +27,5 @@
    USE trabbl          ! bottom boundary layer          (tra_bbl_init routine)
    USE zdfini          ! vertical physics: initialization
+   USE sbcmod          ! surface boundary condition       (sbc_init     routine)
    USE phycst          ! physical constant                  (par_cst routine)
    USE dtadyn          ! Lecture and Interpolation of the dynamical fields
@@ -137,12 +139,37 @@
       !                             !--------------------------------------------!
 #if defined key_iomput
-      CALL init_ioclient( ilocal_comm )       ! nemo local communicator (used or not) given by the io_server
-      narea = mynode( cltxt, ilocal_comm )    ! Nodes selection
+      CALL  init_ioclient( ilocal_comm )                 ! exchange io_server nemo local communicator with the io_server
+      narea = mynode( cltxt, numnam, nstop, ilocal_comm )   ! Nodes selection
 #else
-      narea = mynode( cltxt )                 ! Nodes selection (control print return in cltxt)
+      ilocal_comm = 0
+      narea = mynode( cltxt, numnam, nstop )                 ! Nodes selection (control print return in cltxt)
 #endif
+
       narea = narea + 1                       ! mynode return the rank of proc (0 --> jpnij -1 )
 
       lwp = (narea == 1) .OR. ln_ctl          ! control of all listing output print
+
+      ! 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 nemo_partition(mppsize)
+#else
+      jpni = 1
+      jpnj = 1
+      jpnij = jpni*jpnj
+#endif
+      ! Calculate domain dimensions given calculated jpni and jpnj
+      ! This used to be done in par_oce.F90 when they were parameters rather
+      ! than variables
+      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 nemo_alloc()
 
       IF(lwp) THEN                            ! open listing units
@@ -182,4 +209,5 @@
 
       !                                     ! Ocean physics
+                            CALL     sbc_init   ! Forcings : surface module
 #if ! defined key_degrad
                             CALL ldf_tra_init   ! Lateral ocean tracer physics
@@ -307,4 +335,157 @@
    END SUBROUTINE nemo_closefile
 
+   SUBROUTINE nemo_alloc
+     !!----------------------------------------------------------------------
+     !!                     ***  ROUTINE nemo_alloc  ***
+     !!
+     !! ** Purpose :   Allocate all the dynamic arrays of the OPA modules
+     !!
+     !! ** Method  :
+     !!----------------------------------------------------------------------
+     USE diawri,       ONLY: dia_wri_alloc
+     USE dom_oce,      ONLY: dom_oce_alloc
+     USE zdf_oce,      ONLY: zdf_oce_alloc
+     USE zdfmxl,       ONLY: zdf_mxl_alloc
+     USE ldftra_oce,   ONLY: ldftra_oce_alloc
+     USE trc_oce,      ONLY: trc_oce_alloc
+
+      USE wrk_nemo,    ONLY: wrk_alloc
+
+      INTEGER :: ierr
+      !!----------------------------------------------------------------------
+
+      ierr =        oce_alloc       ()          ! ocean 
+      ierr = ierr + dia_wri_alloc   ()
+      ierr = ierr + dom_oce_alloc   ()          ! ocean domain
+      ierr = ierr + ldftra_oce_alloc()          ! ocean lateral  physics : tracers
+      ierr = ierr + zdf_oce_alloc   ()          ! ocean vertical physics
+      ierr = ierr + zdf_mxl_alloc   ()          ! ocean vertical physics
+      !
+      ierr = ierr + lib_mpp_alloc   (numout)    ! mpp exchanges
+      ierr = ierr + trc_oce_alloc   ()          ! shared TRC / TRA arrays
+      ierr = ierr + wrk_alloc(numout, lwp)
+
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'nemo_alloc : unable to allocate standard ocean arrays' )
+      !
+   END SUBROUTINE nemo_alloc
+
+   SUBROUTINE nemo_partition( num_pes )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE nemo_partition  ***
+      !!
+      !! ** Purpose :   
+      !!
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      USE par_oce
+      INTEGER, INTENT(in) :: num_pes ! The number of MPI processes we have
+      ! Local variables
+      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
+
+      CALL factorise( ifact, nfactmax, nfact, num_pes, ierr )
+
+      IF( nfact <= 1 ) THEN
+         WRITE (numout, *) 'WARNING: factorisation of number of PEs failed'
+         WRITE (numout, *) '       : using grid of ',num_pes,' x 1'
+         jpnj = 1
+         jpni = num_pes
+      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
+      !
+   END SUBROUTINE nemo_partition
+
+
+   SUBROUTINE factorise( ifax, maxfax, nfax, n, ierr )
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE factorise  ***
+      !!
+      !! ** Purpose :   return the prime factors of n.
+      !!                nfax factors are returned in array ifax which is of 
+      !!                maximum dimension maxfax.
+      !! ** Method  :
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in)  :: n, maxfax
+      INTEGER, INTENT(Out) :: ierr, nfax
+      INTEGER, INTENT(out) :: ifax(maxfax)
+      ! Local variables.
+      INTEGER :: i, ifac, l, nu
+      INTEGER, PARAMETER :: ntest = 14
+      INTEGER :: lfax(ntest)
+
+      ! lfax contains the set of allowed factors.
+      data (lfax(i),i=1,ntest) / 16384, 8192, 4096, 2048, 1024, 512, 256,  &
+         &                         128,   64,   32,   16,    8,   4,   2  /
+      !!----------------------------------------------------------------------
+
+      ! Clear the error flag and initialise output vars
+      ierr = 0
+      ifax = 1
+      nfax = 0
+
+      ! Find the factors of n.
+      IF( n == 1 ) GOTO 20
+
+      ! nu holds the unfactorised part of the number.
+      ! nfax holds the number of factors found.
+      ! l points to the allowed factor list.
+      ! ifac holds the current factor.
+
+      nu   = n
+      nfax = 0
+
+      DO l = ntest, 1, -1
+         !
+         ifac = lfax(l)
+         IF(ifac > nu)CYCLE
+
+         ! Test whether the factor will divide.
+
+         IF( MOD(nu,ifac) == 0 ) THEN
+            !
+            nfax = nfax+1            ! Add the factor to the list
+            IF( nfax > maxfax ) THEN
+               ierr = 6
+               write (*,*) 'FACTOR: insufficient space in factor array ',nfax
+               return
+            ENDIF
+            ifax(nfax) = ifac
+            ! Store the other factor that goes with this one
+            nfax = nfax + 1
+            ifax(nfax) = nu / ifac
+            !WRITE (*,*) 'ARPDBG, factors ',nfax-1,' & ',nfax,' are ', &
+            !            ifax(nfax-1),' and ',ifax(nfax)
+         ENDIF
+         !
+      END DO
+
+   20 CONTINUE      ! Label 20 is the exit point from the factor search loop.
+      !
+      RETURN
+      !
+   END SUBROUTINE factorise
    !!======================================================================
 END MODULE nemogcm
Index: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/zdfmxl.F90
===================================================================
--- branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/zdfmxl.F90	(revision 2648)
+++ branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OFF_SRC/zdfmxl.F90	(revision 2648)
@@ -0,0 +1,119 @@
+MODULE zdfmxl
+   !!======================================================================
+   !!                       ***  MODULE  zdfmxl  ***
+   !! Ocean physics: mixed layer depth 
+   !!======================================================================
+   !! History :  1.0  ! 2003-08  (G. Madec)  original code
+   !!            3.2  ! 2009-07  (S. Masson, G. Madec)  IOM + merge of DO-loop
+   !!----------------------------------------------------------------------
+   !!   zdf_mxl      : Compute the turbocline and mixed layer depths.
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE zdf_oce         ! ocean vertical physics
+   USE in_out_manager  ! I/O manager
+   USE prtctl          ! Print control
+   USE iom             ! I/O library
+   USE lib_mpp         ! MPP library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_mxl       ! called by step.F90
+   PUBLIC   zdf_mxl_alloc ! called by nemogcm.F90
+
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   nmln    !: number of level in the mixed layer (used by TOP)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmld    !: mixing layer depth (turbocline)      [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: mixed layer depth at t-points        [m]
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: zdfmxl.F90 2636 2011-03-01 19:04:06Z gm $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_mxl_alloc()
+      !!----------------------------------------------------------------------
+      !!               ***  FUNCTION zdf_mxl_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT=zdf_mxl_alloc)
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( zdf_mxl_alloc )
+      IF( zdf_mxl_alloc /= 0 )   CALL ctl_warn('zdf_mxl_alloc: failed to allocate arrays.')
+   END FUNCTION zdf_mxl_alloc
+
+
+   SUBROUTINE zdf_mxl( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdfmxl  ***
+      !!                   
+      !! ** Purpose :   Compute the turbocline depth and the mixed layer depth
+      !!              with density criteria.
+      !!
+      !! ** Method  :   The mixed layer depth is the shallowest W depth with 
+      !!      the density of the corresponding T point (just bellow) bellow a
+      !!      given value defined locally as rho(10m) + zrho_c
+      !!               The turbocline depth is the depth at which the vertical
+      !!      eddy diffusivity coefficient (resulting from the vertical physics
+      !!      alone, not the isopycnal part, see trazdf.F) fall below a given
+      !!      value defined locally (avt_c here taken equal to 5 cm/s2)
+      !!
+      !! ** Action  :   nmln, hmld, hmlp, hmlpt
+      !!----------------------------------------------------------------------
+      USE wrk_nemo, ONLY: iwrk_in_use, iwrk_not_released
+      USE wrk_nemo, ONLY: imld => iwrk_2d_1    ! 2D integer workspace
+      !!
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !!
+      INTEGER                     ::   ji, jj, jk          ! dummy loop indices
+      INTEGER                     ::   iikn, iiki          ! temporary integer within a do loop
+      REAL(wp)                    ::   zrho_c = 0.01_wp    ! density criterion for mixed layer depth
+      REAL(wp)                    ::   zavt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth
+      !!----------------------------------------------------------------------
+
+      IF( iwrk_in_use(2, 1) )THEN
+         CALL ctl_stop('zdf_mxl : requested workspace array unavailable')   ;   RETURN
+      END IF
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+
+      ! w-level of the mixing and mixed layers
+      nmln(:,:) = mbkt(:,:) + 1        ! Initialization to the number of w ocean point
+      imld(:,:) = mbkt(:,:) + 1
+      DO jk = jpkm1, nlb10, -1         ! from the bottom to nlb10
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + zrho_c )   nmln(ji,jj) = jk      ! Mixed layer
+               IF( avt (ji,jj,jk) < zavt_c                     )   imld(ji,jj) = jk      ! Turbocline 
+            END DO
+         END DO
+      END DO
+      ! depth of the mixing and mixed layers
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            iiki = imld(ji,jj)
+            iikn = nmln(ji,jj)
+            hmld (ji,jj) = fsdepw(ji,jj,iiki  ) * tmask(ji,jj,1)    ! Turbocline depth 
+            hmlp (ji,jj) = fsdepw(ji,jj,iikn  ) * tmask(ji,jj,1)    ! Mixed layer depth
+            hmlpt(ji,jj) = fsdept(ji,jj,iikn-1)                     ! depth of the last T-point inside the mixed layer
+         END DO
+      END DO
+      CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
+      CALL iom_put( "mldkz5"  , hmld )   ! turbocline depth
+      
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 )
+      !
+      IF( iwrk_not_released(2, 1) )   CALL ctl_stop('zdf_mxl : failed to release workspace array')
+      !
+   END SUBROUTINE zdf_mxl
+
+   !!======================================================================
+END MODULE zdfmxl
