Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 2715)
@@ -4,4 +4,8 @@
    !! LIM ice model : horizontal diffusion of sea-ice quantities
    !!======================================================================
+   !! History :  LIM  !  2000-01 (LIM) Original code
+   !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm
+   !!            1.0  !  2002-08 (C. Ethe)  F90, free form
+   !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
@@ -10,31 +14,27 @@
    !!   lim_hdf  : diffusion trend on sea-ice variable
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE dom_oce
-   USE in_out_manager
-   USE ice
-   USE lbclnk
-   USE lib_mpp
-   USE prtctl          ! Print control
+   USE dom_oce          ! ocean domain
+   USE ice              ! LIM-3: ice variables
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE prtctl           ! Print control
+   USE in_out_manager   ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_hdf    ! called by lim_tra
+   PUBLIC   lim_hdf     ! called by lim_tra
 
-   !! * Module variables
-   LOGICAL  ::   linit = .TRUE.              ! ???
+   LOGICAL  ::   linit = .TRUE.              ! initialization flag (set to flase after the 1st call)
    REAL(wp) ::   epsi04 = 1e-04              ! constant
-   REAL(wp), DIMENSION(jpi,jpj) ::   zfact   ! ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   efact   ! metric coefficient
 
    !! * Substitution 
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -43,79 +43,63 @@
       !!                  ***  ROUTINE lim_hdf  ***
       !!
-      !! ** purpose :   Compute and add the diffusive trend on sea-ice
-      !!      variables
+      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
       !!
       !! ** method  :   Second order diffusive operator evaluated using a
-      !!      Cranck-Nicholson time Scheme.
+      !!              Cranck-Nicholson time Scheme.
       !!
       !! ** Action  :    update ptab with the diffusive contribution
-      !!
-      !! History :
-      !!        !  00-01 (LIM) Original code
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!        !  02-08 (C. Ethe)  F90, free form
       !!-------------------------------------------------------------------
-      ! * Arguments
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   &
-         ptab                 ! Field on which the diffusion is applied  
-      REAL(wp), DIMENSION(jpi,jpj) ::   &
-         ptab0                ! ???
-
-      ! * Local variables
-      INTEGER ::  ji, jj      ! dummy loop indices
-      INTEGER ::  &
-         its, iter            ! temporary integers
-      CHARACTER (len=55) :: charout
-      REAL(wp) ::  &
-         zalfa, zrlxint, zconv, zeps   ! temporary scalars
-      REAL(wp), DIMENSION(jpi,jpj) ::  & 
-         zrlx, zflu, zflv, &  ! temporary workspaces
-         zdiv0, zdiv          !    "           "
+      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
+      USE wrk_nemo, ONLY:   zflu => wrk_2d_11, zdiv  => wrk_2d_13, zrlx  => wrk_2d_15 
+      USE wrk_nemo, ONLY:   zflv => wrk_2d_12, zdiv0 => wrk_2d_14, ztab0 => wrk_2d_16
+      !
+      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
+      !
+      INTEGER  ::  ji, jj            ! dummy loop indices
+      INTEGER  ::  its, iter, ierr   ! local integers
+      REAL(wp) ::   zalfa, zrlxint, zconv, zeps   ! local scalars
+      CHARACTER(lc) ::   charout   ! local character
       !!-------------------------------------------------------------------
-
-      ! Initialisation
-      ! ---------------   
-      ! Time integration parameters
-      zalfa = 0.5       ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
-      its   = 100       ! Maximum number of iteration
-      zeps  =  2. * epsi04
-
-      ! Arrays initialization
-      ptab0 (:, : ) = ptab(:,:)
-      !bug  zflu (:,jpj) = 0.e0
-      !bug  zflv (:,jpj) = 0.e0
-      zdiv0(:, 1 ) = 0.e0
-      zdiv0(:,jpj) = 0.e0
-      IF( .NOT.lk_vopt_loop ) THEN
-         zflu (jpi,:) = 0.e0   
-         zflv (jpi,:) = 0.e0
-         zdiv0(1,  :) = 0.e0
-         zdiv0(jpi,:) = 0.e0
+      
+      IF( wrk_in_use(2, 11,12,13,14,15,16) ) THEN
+         CALL ctl_stop( 'lim_hdf: requested workspace arrays unavailable' )   ;   RETURN
       ENDIF
 
-      ! Metric coefficient (compute at the first call and saved in
-      IF( linit ) THEN
+      !                       !==  Initialisation  ==!
+      !
+      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
+         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
+         IF( lk_mpp    )   CALL mpp_sum( ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
          DO jj = 2, jpjm1  
             DO ji = fs_2 , fs_jpim1   ! vector opt.
-               zfact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj  ) + e1v(ji,jj) + e1v(ji,jj-1) ) &
-                  &          / ( e1t(ji,jj) * e2t(ji,jj) )
+               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
             END DO
          END DO
          linit = .FALSE.
       ENDIF
+      !                             ! Time integration parameters
+      zalfa = 0.5_wp                      ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
+      its   = 100                         ! Maximum number of iteration
+      zeps  =  2._wp * epsi04
+      !
+      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
+      zdiv0(:, 1 ) = 0._wp
+      zdiv0(:,jpj) = 0._wp
+      IF( .NOT.lk_vopt_loop ) THEN
+         zflu (jpi,:) = 0._wp   
+         zflv (jpi,:) = 0._wp
+         zdiv0(1,  :) = 0._wp
+         zdiv0(jpi,:) = 0._wp
+      ENDIF
 
-
-      ! Sub-time step loop
-      zconv = 1.e0
+      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
       iter  = 0
-
-      !                                                   !===================
-      DO WHILE ( ( zconv > zeps ) .AND. (iter <= its) )   ! Sub-time step loop
-         !                                                !===================
-         ! incrementation of the sub-time step number
-         iter = iter + 1
-
-         ! diffusive fluxes in U- and V- direction
-         DO jj = 1, jpjm1
+      !
+      DO WHILE( zconv > zeps .AND. iter <= its )   ! Sub-time step loop
+         !
+         iter = iter + 1                                 ! incrementation of the sub-time step number
+         !
+         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
             DO ji = 1 , fs_jpim1   ! vector opt.
                zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
@@ -123,7 +107,6 @@
             END DO
          END DO
-
-         ! diffusive trend : divergence of the fluxes
-         DO jj= 2, jpjm1
+         !
+         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
             DO ji = fs_2 , fs_jpim1   ! vector opt. 
                zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
@@ -131,24 +114,19 @@
             END DO
          END DO
-
-         ! save the first evaluation of the diffusive trend in zdiv0
-         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)       
-
-         ! XXXX iterative evaluation?????
-         DO jj = 2, jpjm1
+         !
+         IF( iter == 1 )   zdiv0(:,:) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
+         !
+         DO jj = 2, jpjm1                                ! iterative evaluation
             DO ji = fs_2 , fs_jpim1   ! vector opt.
-               zrlxint = (   ptab0(ji,jj)    &
-                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + zfact(ji,jj) * ptab(ji,jj) )   &
+               zrlxint = (   ztab0(ji,jj)    &
+                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
                   &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
-                  &    / ( 1.0 + zalfa * rdt_ice * zfact(ji,jj) )
+                  &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
                zrlx(ji,jj) = ptab(ji,jj) + om * ( zrlxint - ptab(ji,jj) )
             END DO
          END DO
-
-         ! lateral boundary condition on zrlx
-         CALL lbc_lnk( zrlx, 'T', 1. )
-
-         ! convergence test
-         zconv = 0.e0
+         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
+         !
+         zconv = 0._wp                                   ! convergence test
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1
@@ -156,23 +134,18 @@
             END DO
          END DO
-         IF( lk_mpp )   CALL mpp_max( zconv )   ! max over the global domain
-
-         DO jj = 1, jpj
-            DO ji = 1 , jpi
-               ptab(ji,jj) = zrlx(ji,jj)
-            END DO
-         END DO
-
-         !                                         !==========================
+         IF( lk_mpp )   CALL mpp_max( zconv )            ! max over the global domain
+         !
+         ptab(:,:) = zrlx(:,:)
+         !
       END DO                                       ! end of sub-time step loop
-      !                                            !==========================
 
       IF(ln_ctl)   THEN
-         zrlx(:,:) = ptab(:,:) - ptab0(:,:)
+         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
          WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
-         CALL prt_ctl(tab2d_1=zrlx, clinfo1=charout)
+         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
       ENDIF
-
-
+      !
+      IF( wrk_not_released(2, 11,12,13,14,15,16) )   CALL ctl_stop('lim_hdf: failed to release workspace arrays')
+      !
    END SUBROUTINE lim_hdf
 
