Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 2620)
@@ -47,7 +47,6 @@
    !! NEMO/OPA 3.3 , 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
 
@@ -66,6 +65,5 @@
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY: wrk_use, wrk_release, llwrk_use, llwrk_release
-      USE wrk_nemo, ONLY: llwrk_3d_1 
-      USE wrk_nemo, ONLY: wrk_3d_6, wrk_3d_7
+      USE wrk_nemo, ONLY: wrk_3d_6, wrk_3d_7    ! 3D workspace
       !!
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -86,5 +84,4 @@
       REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2
       !!
-      LOGICAL,  POINTER, DIMENSION(:,:,:) ::   llmask    ! Pointer to sub-array of workspace array
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness
@@ -93,18 +90,10 @@
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
 
-      IF( (.not. llwrk_use(3,1)) .OR. (.not. wrk_use(3, 6,7)) )THEN
-         CALL ctl_stop('albedo_ice: requested workspace arrays are unavailable.')
-         RETURN
-      ELSE IF(ijpl > jpk)THEN
-         ! 3D workspace arrays have extent jpk in 3rd dimension - check that 
-         ! ijpl doesn't exceed it.
-         CALL ctl_stop('albedo_ice: 3rd dimension of standard workspace arrays too small for them to be used here.')
-         RETURN
-      ELSE
-         ! Associate pointers with sub-arrays of workspace arrays
-         llmask  => llwrk_3d_1(:,:,1:ijpl)
-         zalbfz  =>   wrk_3d_6(:,:,1:ijpl)
-         zficeth =>   wrk_3d_7(:,:,1:ijpl)
-      END IF
+      IF( (.not. llwrk_use(3,1)) .OR. (.not. wrk_use(3, 6,7)) ) THEN
+         CALL ctl_stop('albedo_ice: requested workspace arrays are unavailable')   ;   RETURN
+      ENDIF
+      ! Associate pointers with sub-arrays of workspace arrays
+      zalbfz  =>   wrk_3d_6(:,:,1:ijpl)
+      zficeth =>   wrk_3d_7(:,:,1:ijpl)
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
@@ -113,27 +102,36 @@
       !  Computation of  zficeth
       !---------------------------
-      llmask(:,:,1:ijpl) = ( ph_snw == 0.e0 ) .AND. ( pt_ice >= rt0_ice )
       ! ice free of snow and melts
-      WHERE( llmask(:,:,1:ijpl) )   ;   zalbfz = rn_albice
-      ELSEWHERE                     ;   zalbfz = rn_alphdi
+      WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
+      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
+      END  WHERE
+
+      WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
+      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
+      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
+         &                                                                 - 0.8608 * ph_ice * ph_ice     &
+         &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
       END WHERE
 
-      DO jl = 1, ijpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
-                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
-               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
-                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
-               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
-                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
-                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
-                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
-               ELSE
-                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
-               ENDIF
-            END DO
-         END DO
-      END DO
+!!gm old code
+!      DO jl = 1, ijpl
+!         DO jj = 1, jpj
+!            DO ji = 1, jpi
+!               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
+!                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
+!               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
+!                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
+!               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
+!                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
+!                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
+!                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
+!               ELSE
+!                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
+!               ENDIF
+!            END DO
+!         END DO
+!      END DO
+!!gm end old code
       
       !----------------------------------------------- 
@@ -174,7 +172,5 @@
       pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction
       !
-      IF( (.not. llwrk_release(3, 1)) .OR. (.not. wrk_release(3, 6,7)) )THEN
-         CALL ctl_stop('albedo_ice: failed to release workspace arrays.')
-      END IF
+      IF( .not. wrk_release(3, 6,7) )   CALL ctl_stop('albedo_ice: failed to release workspace arrays')
       !
    END SUBROUTINE albedo_ice
@@ -186,6 +182,4 @@
       !! 
       !! ** Purpose :   Computation of the albedo of the ocean
-      !!
-      !! ** Method  :   ....
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_os   !  albedo of ocean under overcast sky
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis3.F90	(revision 2620)
@@ -6,6 +6,6 @@
    !!=====================================================================
    !! History :   
-   !!   9.0  !  04-06  (R. Redler, NEC Laboratories Europe, St Augustin, Germany) Original code
-   !!   " "  !  04-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Kiel, Germany) revision
+   !!   9.0  !  04-06  (R. Redler, NEC Laboratories Europe, Germany) Original code
+   !!   " "  !  04-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Germany) revision
    !!   " "  !  04-11  (V. Gayler, MPI M&D) Grid writing
    !!   " "  !  05-08  (R. Redler, W. Park) frld initialization, paral(2) revision
@@ -17,5 +17,4 @@
    !!----------------------------------------------------------------------
    !!   'key_oasis3'                    coupled Ocean/Atmosphere via OASIS3
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   cpl_prism_init     : initialization of coupled mode communication
@@ -34,10 +33,18 @@
    USE in_out_manager               ! I/O manager
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+
    IMPLICIT NONE
    PRIVATE
-!
-   LOGICAL, PUBLIC, PARAMETER :: lk_cpl = .TRUE.   ! coupled flag
-   INTEGER, PUBLIC            :: OASIS_Rcv  = 1    ! return code if received field
-   INTEGER, PUBLIC            :: OASIS_idle = 0    ! return code if nothing done by oasis
+
+   PUBLIC cpl_prism_init
+   PUBLIC cpl_prism_define
+   PUBLIC cpl_prism_snd
+   PUBLIC cpl_prism_rcv
+   PUBLIC cpl_prism_freq
+   PUBLIC cpl_prism_finalize
+
+   LOGICAL, PUBLIC, PARAMETER :: lk_cpl = .TRUE.   !: coupled flag
+   INTEGER, PUBLIC            :: OASIS_Rcv  = 1    !: return code if received field
+   INTEGER, PUBLIC            :: OASIS_idle = 0    !: return code if nothing done by oasis
    INTEGER                    :: ncomp_id          ! id returned by prism_init_comp
    INTEGER                    :: nerror            ! return error code
@@ -45,5 +52,5 @@
    INTEGER, PARAMETER :: nmaxfld=40    ! Maximum number of coupling fields
    
-   TYPE, PUBLIC ::   FLD_CPL            ! Type for coupling field information
+   TYPE, PUBLIC ::   FLD_CPL            !: Type for coupling field information
       LOGICAL            ::   laction   ! To be coupled or not
       CHARACTER(len = 8) ::   clname    ! Name of the coupling field   
@@ -53,26 +60,16 @@
    END TYPE FLD_CPL
 
-   TYPE(FLD_CPL), DIMENSION(nmaxfld), PUBLIC :: srcv, ssnd   ! Coupling fields
+   TYPE(FLD_CPL), DIMENSION(nmaxfld), PUBLIC :: srcv, ssnd   !: Coupling fields
 
    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: exfld  ! Temporary buffer for receiving
-
-   !! Routine accessibility
-   PUBLIC cpl_prism_init
-   PUBLIC cpl_prism_define
-   PUBLIC cpl_prism_snd
-   PUBLIC cpl_prism_rcv
-   PUBLIC cpl_prism_freq
-   PUBLIC cpl_prism_finalize
 
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , 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
 
-   SUBROUTINE cpl_prism_init (kl_comm) 
-
+   SUBROUTINE cpl_prism_init( kl_comm )
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE cpl_prism_init  ***
@@ -83,5 +80,5 @@
       !! ** Method  :   OASIS3 MPI communication 
       !!--------------------------------------------------------------------
-      INTEGER, INTENT(   OUT )   :: kl_comm       ! local communicator of the model
+      INTEGER, INTENT(out) ::   kl_comm   ! local communicator of the model
       !!--------------------------------------------------------------------
 
@@ -103,10 +100,9 @@
       IF ( nerror /= PRISM_Ok ) &
          CALL prism_abort_proto (ncomp_id, 'cpl_prism_init','Failure in prism_get_localcomm_proto' )
-
+      !
    END SUBROUTINE cpl_prism_init
 
 
-   SUBROUTINE cpl_prism_define (krcv, ksnd)
-
+   SUBROUTINE cpl_prism_define( krcv, ksnd )
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE cpl_prism_define  ***
@@ -117,11 +113,10 @@
       !! ** Method  :   OASIS3 MPI communication 
       !!--------------------------------------------------------------------
-      INTEGER, INTENT( IN    )   :: krcv, ksnd     ! Number of received and sent coupling fields
-      !
-      INTEGER                    :: id_part
-      INTEGER                    :: paral(5)       ! OASIS3 box partition
-      INTEGER                    :: ishape(2,2)    ! shape of arrays passed to PSMILe
-      INTEGER                    :: ji             ! local loop indicees
-      !!
+      INTEGER, INTENT(in) ::   krcv, ksnd     ! Number of received and sent coupling fields
+      !
+      INTEGER :: id_part
+      INTEGER :: paral(5)       ! OASIS3 box partition
+      INTEGER :: ishape(2,2)    ! shape of arrays passed to PSMILe
+      INTEGER :: ji             ! local loop indicees
       !!--------------------------------------------------------------------
 
@@ -142,7 +137,6 @@
       !
       ALLOCATE(exfld(nlei-nldi+1, nlej-nldj+1), stat = nerror)
-      IF (nerror > 0) THEN
-         CALL prism_abort_proto ( ncomp_id, 'cpl_prism_define', 'Failure in allocating exfld')
-         RETURN
+      IF( nerror > 0 ) THEN
+         CALL prism_abort_proto ( ncomp_id, 'cpl_prism_define', 'Failure in allocating exfld')   ;   RETURN
       ENDIF
       !
@@ -197,11 +191,10 @@
       
       CALL prism_enddef_proto(nerror)
-      IF ( nerror /= PRISM_Ok )   CALL prism_abort_proto ( ncomp_id, 'cpl_prism_define', 'Failure in prism_enddef')
-      
+      IF( nerror /= PRISM_Ok )   CALL prism_abort_proto ( ncomp_id, 'cpl_prism_define', 'Failure in prism_enddef')
+      !
    END SUBROUTINE cpl_prism_define
    
    
    SUBROUTINE cpl_prism_snd( kid, kstep, pdata, kinfo )
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_snd  ***
@@ -210,12 +203,8 @@
       !!      like sst or ice cover to the coupler or remote application.
       !!----------------------------------------------------------------------
-      !! * Arguments
-      !!
-      INTEGER,                      INTENT( IN    )   :: kid       ! variable index in the array
-      INTEGER,                      INTENT(   OUT )   :: kinfo     ! OASIS3 info argument
-      INTEGER,                      INTENT( IN    )   :: kstep     ! ocean time-step in seconds
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN    )   :: pdata
-      !!
-      !!
+      INTEGER                 , INTENT(in   ) ::   kid       ! variable index in the array
+      INTEGER                 , INTENT(  out) ::   kinfo     ! OASIS3 info argument
+      INTEGER                 , INTENT(in   ) ::   kstep     ! ocean time-step in seconds
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdata
       !!--------------------------------------------------------------------
       !
@@ -236,11 +225,11 @@
             WRITE(numout,*) '     -     Sum value is ', SUM(pdata)
             WRITE(numout,*) '****************'
-        ENDIF
-     ENDIF
+         ENDIF
+      ENDIF
+      !
     END SUBROUTINE cpl_prism_snd
 
 
    SUBROUTINE cpl_prism_rcv( kid, kstep, pdata, kinfo )
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_rcv  ***
@@ -249,8 +238,8 @@
       !!      like stresses and fluxes from the coupler or remote application.
       !!----------------------------------------------------------------------
-      INTEGER,                      INTENT( IN    )   :: kid       ! variable index in the array
-      INTEGER,                      INTENT( IN    )   :: kstep     ! ocean time-step in seconds
-      REAL(wp),     DIMENSION(:,:), INTENT( INOUT )   :: pdata     ! IN to keep the value if nothing is done
-      INTEGER,                      INTENT(   OUT )   :: kinfo     ! OASIS3 info argument
+      INTEGER                 , INTENT(in   ) ::   kid       ! variable index in the array
+      INTEGER                 , INTENT(in   ) ::   kstep     ! ocean time-step in seconds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdata     ! IN to keep the value if nothing is done
+      INTEGER                 , INTENT(  out) ::   kinfo     ! OASIS3 info argument
       !!
       LOGICAL                :: llaction
@@ -291,10 +280,9 @@
          kinfo = OASIS_idle     
       ENDIF
-
+      !
    END SUBROUTINE cpl_prism_rcv
 
 
    FUNCTION cpl_prism_freq( kid )  
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_freq  ***
@@ -304,11 +292,11 @@
       INTEGER,INTENT( IN )   :: kid              ! variable index 
       INTEGER                :: cpl_prism_freq   ! coupling frequency
+      !!----------------------------------------------------------------------
       cpl_prism_freq = ig_def_freq( kid )
-
+      !
    END FUNCTION cpl_prism_freq
 
 
    SUBROUTINE cpl_prism_finalize
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_finalize  ***
@@ -318,16 +306,13 @@
       !!      MPI communication.
       !!----------------------------------------------------------------------
-
+      !
       DEALLOCATE(exfld)
       CALL prism_terminate_proto ( nerror )         
-
+      !
    END SUBROUTINE cpl_prism_finalize
 
 #else
-
-   !!----------------------------------------------------------------------
-   !!   Default case                                Forced Ocean/Atmosphere
-   !!----------------------------------------------------------------------
-   !!   Empty module
+   !!----------------------------------------------------------------------
+   !!   Default case          Dummy module          Forced Ocean/Atmosphere
    !!----------------------------------------------------------------------
    USE in_out_manager               ! I/O manager
@@ -335,18 +320,15 @@
    PUBLIC cpl_prism_init
    PUBLIC cpl_prism_finalize
-
 CONTAINS
-
    SUBROUTINE cpl_prism_init (kl_comm) 
-      INTEGER, INTENT(   OUT )   :: kl_comm       ! local communicator of the model
+      INTEGER, INTENT(out)   :: kl_comm       ! local communicator of the model
       kl_comm = -1
       WRITE(numout,*) 'cpl_prism_init: Error you sould not be there...'
    END SUBROUTINE cpl_prism_init
-
    SUBROUTINE cpl_prism_finalize
       WRITE(numout,*) 'cpl_prism_finalize: Error you sould not be there...'
    END SUBROUTINE cpl_prism_finalize
-
 #endif
 
+   !!=====================================================================
 END MODULE cpl_oasis3
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis4.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis4.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/cpl_oasis4.F90	(revision 2620)
@@ -5,17 +5,16 @@
    !!=====================================================================
    !! History :   
-   !!   9.0  !  04-06  (R. Redler, NEC Laboratories Europe, St Augustin, Germany) Original code
-   !!   " "  !  04-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Kiel, Germany) revision
-   !!   " "  !  04-11  (V. Gayler, MPI M&D) Grid writing
-   !!   " "  !  05-08  (R. Redler, W. Park) frld initialization, paral(2) revision
-   !!   " "  !  05-09  (R. Redler) extended to allow for communication over root only
-   !!   " "  !  06-01  (W. Park) modification of physical part
-   !!   " "  !  06-02  (R. Redler, W. Park) buffer array fix for root exchange
-   !!   " "  !  2010   (E. Maisonnave and S. Masson) complete rewrite
+   !!   9.0  !  2004-06  (R. Redler, NEC Laboratories Europe, St Augustin, Germany) Original code
+   !!    -   !  2004-11  (R. Redler, NEC Laboratories Europe; N. Keenlyside, W. Park, IFM-GEOMAR, Kiel, Germany) revision
+   !!    -   !  2004-11  (V. Gayler, MPI M&D) Grid writing
+   !!    -   !  2005-08  (R. Redler, W. Park) frld initialization, paral(2) revision
+   !!    -   !  2005-09  (R. Redler) extended to allow for communication over root only
+   !!    -   !  2006-01  (W. Park) modification of physical part
+   !!    -   !  2006-02  (R. Redler, W. Park) buffer array fix for root exchange
+   !!    -   !  2010-10  (E. Maisonnave and S. Masson) complete rewrite
    !!----------------------------------------------------------------------
 #if defined key_oasis4
    !!----------------------------------------------------------------------
    !!   'key_oasis4'                    coupled Ocean/Atmosphere via OASIS4
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   cpl_prism_init     : initialization of coupled mode communication
@@ -35,5 +34,12 @@
    IMPLICIT NONE
    PRIVATE
-!
+
+   PUBLIC cpl_prism_init
+   PUBLIC cpl_prism_define
+   PUBLIC cpl_prism_snd
+   PUBLIC cpl_prism_rcv
+   PUBLIC cpl_prism_update_time
+   PUBLIC cpl_prism_finalize
+   
 !   LOGICAL, PUBLIC, PARAMETER :: lk_cpl = .TRUE.    ! coupled flag
    INTEGER                    :: ncomp_id           ! id returned by prism_init_comp
@@ -59,23 +65,12 @@
    TYPE(PRISM_Time_struct), PUBLIC    :: date_bound(2)   ! date info for send operation
 
-
-   !! Routine accessibility
-   PUBLIC cpl_prism_init
-   PUBLIC cpl_prism_define
-   PUBLIC cpl_prism_snd
-   PUBLIC cpl_prism_rcv
-   PUBLIC cpl_prism_update_time
-   PUBLIC cpl_prism_finalize
-
-   !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
-   !! $Header$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE cpl_prism_init (kl_comm) 
-
+   SUBROUTINE cpl_prism_init( kl_comm ) 
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE cpl_prism_init  ***
@@ -86,6 +81,6 @@
       !! ** Method  :   OASIS4 MPI communication 
       !!--------------------------------------------------------------------
-      INTEGER, INTENT(   OUT )   :: kl_comm       ! local communicator of the model
-      !
+      INTEGER, INTENT(out) ::   kl_comm   ! local communicator of the model
+      !!--------------------------------------------------------------------
       
       CALL prism_init( 'nemo', nerror )
@@ -102,11 +97,9 @@
       CALL prism_get_localcomm( ncomp_id, kl_comm, nerror )
       IF( nerror /= PRISM_Success )   CALL prism_abort( ncomp_id, 'cpl_prism_init', 'Failure in prism_get_localcomm' )
-
-
+      !
    END SUBROUTINE cpl_prism_init
 
 
-   SUBROUTINE cpl_prism_define (krcv, ksnd)
-
+   SUBROUTINE cpl_prism_define( krcv, ksnd )
       !!-------------------------------------------------------------------
       !!             ***  ROUTINE cpl_prism_define  ***
@@ -120,23 +113,20 @@
       USE wrk_nemo, ONLY: zclo => wrk_3d_1, zcla => wrk_3d_2
       USE wrk_nemo, ONLY: zlon => wrk_2d_1, zlat => wrk_2d_2
-      !!
-      INTEGER, INTENT( IN    )   :: krcv, ksnd     ! Number of received and sent coupling fields
+      !
+      INTEGER, INTENT(in) :: krcv, ksnd     ! Number of received and sent coupling fields
       !
       INTEGER, DIMENSION(4)      :: igrid     ! ids returned by prism_def_grid
       INTEGER, DIMENSION(4)      :: iptid     ! ids returned by prism_set_points
-
-      INTEGER, DIMENSION(4)      :: imskid  ! ids returned by prism_set_mask
-      INTEGER, DIMENSION(4)      :: iishift     ! 
-      INTEGER, DIMENSION(4)      :: ijshift    ! 
+      INTEGER, DIMENSION(4)      :: imskid    ! ids returned by prism_set_mask
+      INTEGER, DIMENSION(4)      :: iishift   ! 
+      INTEGER, DIMENSION(4)      :: ijshift   ! 
       INTEGER, DIMENSION(4)      :: iioff     ! 
-      INTEGER, DIMENSION(4)      :: ijoff    ! 
-      INTEGER, DIMENSION(4)      :: itmp    ! 
-      INTEGER, DIMENSION(1,3)    :: iextent    ! 
-      INTEGER, DIMENSION(1,3)    :: ioffset    ! 
-
-
-      INTEGER                    :: ishape(2,3)     ! shape of arrays passed to PSMILe
+      INTEGER, DIMENSION(4)      :: ijoff     ! 
+      INTEGER, DIMENSION(4)      :: itmp      ! 
+      INTEGER, DIMENSION(1,3)    :: iextent   ! 
+      INTEGER, DIMENSION(1,3)    :: ioffset   ! 
+
+      INTEGER                    :: ishape(2,3)    ! shape of arrays passed to PSMILe
       INTEGER                    :: data_type      ! data type of transients
-
 
       LOGICAL                    :: new_points
@@ -144,7 +134,7 @@
       LOGICAL, ALLOCATABLE, SAVE :: llmask(:,:,:) ! jpi,jpj,1
 
-      INTEGER                    :: ji, jj, jg, jc     ! local loop indicees
-      INTEGER                    :: ii,ij     ! index
-      INTEGER, DIMENSION(1)      :: ind     ! index
+      INTEGER                    :: ji, jj, jg, jc   ! local loop indicees
+      INTEGER                    :: ii, ij           ! index
+      INTEGER, DIMENSION(1)      :: ind              ! index
 
       CHARACTER(len=32)          :: clpt_name     ! name of the grid points
@@ -154,11 +144,9 @@
       TYPE(PRISM_Time_struct)    :: tmpdate
       INTEGER                    :: idate_incr      ! date increment
-      !!
-      !!--------------------------------------------------------------------
-
-      IF( (.not. wrk_use(3, 1,2)) .OR. (.not. wrk_use(2, 1,2)) )THEN
-         CALL ctl_stop('cpl_prism_define: ERROR: requested workspace arrays are unavailable.')
-         RETURN
-      END IF
+      !!--------------------------------------------------------------------
+
+      IF( .NOT. wrk_use(3, 1,2) .OR. .NOT. wrk_use(2, 1,2) )THEN
+         CALL ctl_stop('cpl_prism_define: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
+      ENDIF
 
       IF(lwp) WRITE(numout,*)
@@ -333,13 +321,11 @@
       IF ( nerror /= PRISM_Success )   CALL prism_abort ( ncomp_id, 'cpl_prism_define', 'Failure in prism_enddef')
       
-      IF( (.not. wrk_release(3, 1,2)) .OR. (.not. wrk_release(2, 1,2)) )THEN
-         CALL ctl_stop('cpl_prism_define: ERROR: failed to release workspace arrays.')
-      END IF
-
+      IF( .not. wrk_release(3, 1,2) .OR.   &
+          .not. wrk_release(2, 1,2)   )   CALL ctl_stop('cpl_prism_define: failed to release workspace arrays')
+      !
    END SUBROUTINE cpl_prism_define
    
    
    SUBROUTINE cpl_prism_snd( kid, kstep, pdata, kinfo )
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_snd  ***
@@ -348,12 +334,8 @@
       !!      like sst or ice cover to the coupler or remote application.
       !!----------------------------------------------------------------------
-      !! * Arguments
-      !!
-      INTEGER,                      INTENT( IN    )   :: kid       ! variable intex in the array
-      INTEGER,                      INTENT(   OUT )   :: kinfo     ! OASIS4 info argument
-      INTEGER,                      INTENT( IN    )   :: kstep     ! ocean time-step in seconds
-      REAL(wp),     DIMENSION(:,:), INTENT( IN    )   :: pdata
-      !!
-      !!
+      INTEGER                 , INTENT(in   ) ::   kid     ! variable intex in the array
+      INTEGER                 , INTENT(  out) ::   kinfo   ! OASIS4 info argument
+      INTEGER                 , INTENT(in   ) ::   kstep   ! ocean time-step in seconds
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdata
       !!--------------------------------------------------------------------
       !
@@ -365,5 +347,5 @@
          &                                               'Failure in prism_put for '//TRIM(ssnd(kid)%clname) )
 
-      IF ( ln_ctl ) THEN        
+      IF( ln_ctl ) THEN        
          IF ( kinfo >= PRISM_Cpl     .OR. kinfo == PRISM_Rst .OR.   &
             & kinfo == PRISM_RstTimeop ) THEN
@@ -377,11 +359,11 @@
             WRITE(numout,*) '     -     Sum value is ', SUM(pdata)
             WRITE(numout,*) '****************'
-        ENDIF
-     ENDIF
-    END SUBROUTINE cpl_prism_snd
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE cpl_prism_snd
 
 
    SUBROUTINE cpl_prism_rcv( kid, kstep, pdata, kinfo )
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_rcv  ***
@@ -390,9 +372,9 @@
       !!      like stresses and fluxes from the coupler or remote application.
       !!----------------------------------------------------------------------
-      INTEGER,                      INTENT( IN    )   :: kid       ! variable intex in the array
-      INTEGER,                      INTENT( IN    )   :: kstep     ! ocean time-step in seconds
-      REAL(wp),     DIMENSION(:,:), INTENT( INOUT )   :: pdata     ! IN to keep the value if nothing is done
-      INTEGER,                      INTENT(   OUT )   :: kinfo     ! OASIS4 info argument
-      !!
+      INTEGER                 , INTENT(in   ) ::   kid     ! variable intex in the array
+      INTEGER                 , INTENT(in   ) ::   kstep   ! ocean time-step in seconds
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pdata   ! IN to keep the value if nothing is done
+      INTEGER                 , INTENT(  out) ::   kinfo   ! OASIS4 info argument
+      !
       LOGICAL                :: llaction
       !!--------------------------------------------------------------------
@@ -435,11 +417,9 @@
          kinfo = OASIS_idle     
       ENDIF
-
-
+      !
    END SUBROUTINE cpl_prism_rcv
 
 
    SUBROUTINE cpl_prism_finalize
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_finalize  ***
@@ -449,33 +429,30 @@
       !!      MPI communication.
       !!----------------------------------------------------------------------
-
+      !
       DEALLOCATE(exfld)
       CALL prism_terminate ( nerror )         
-
+      !
    END SUBROUTINE cpl_prism_finalize
 
+
    SUBROUTINE cpl_prism_update_time(kt)
-
       !!---------------------------------------------------------------------
       !!              ***  ROUTINE cpl_prism_update_time  ***
       !!
       !! ** Purpose : - Increment date with model timestep
-      !!      called explicitly at the end of each timestep
+      !!                called explicitly at the end of each timestep
       !!----------------------------------------------------------------------
-
-      INTEGER, INTENT(in) ::   kt       ! ocean model time step index
-
-      TYPE(PRISM_Time_struct)    :: tmpdate
-      INTEGER                    :: idate_incr     ! date increment
-
-
-      IF( kt == nit000 ) THEN
-      !
-      ! Define the actual date 
-      !
-      ! date is determined by adding days since beginning of the run to the corresponding initial date.
-      ! Note that OPA internal info about the start date of the experiment is bypassed.
-      ! Instead we rely sololy on the info provided by the SCC.xml file.
-      !
+      INTEGER, INTENT(in) ::   kt   ! ocean model time step index
+
+      TYPE(PRISM_Time_struct) ::   tmpdate
+      INTEGER                 ::   idate_incr   ! date increment
+      !!----------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN      ! Define the actual date 
+         !
+         ! date is determined by adding days since beginning of the run to the corresponding initial date.
+         ! Note that OPA internal info about the start date of the experiment is bypassed.
+         ! Instead we rely sololy on the info provided by the SCC.xml file.
+         !
          date = PRISM_Jobstart_date
          !
@@ -486,9 +463,7 @@
          tmpdate = date   ;   CALL PRISM_calc_newdate ( tmpdate, -idate_incr, nerror )   ;   date_bound(1) = tmpdate
          tmpdate = date   ;   CALL PRISM_calc_newdate ( tmpdate,  idate_incr, nerror )   ;   date_bound(2) = tmpdate
-
-      ELSE
-      !
-      ! Date update
-      !
+         !
+      ELSE      ! Date update
+         !
          idate_incr  = rdttra(1)
          CALL PRISM_calc_newdate( date, idate_incr, nerror )
@@ -497,10 +472,11 @@
          CALL PRISM_calc_newdate( tmpdate, idate_incr, nerror )
          date_bound(2) = tmpdate
-
+         !
       END IF
-
+      !
    END SUBROUTINE cpl_prism_update_time
 
 #endif
 
+   !!=====================================================================
 END MODULE cpl_oasis4
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/fldread.F90	(revision 2620)
@@ -4,7 +4,7 @@
    !! Ocean forcing:  read input field for surface boundary condition
    !!=====================================================================
-   !! History :  9.0  !  06-06  (G. Madec) Original code
-   !!                 !  05-08  (S. Alderson) Modified for Interpolation in memory
-   !!                 !         from input grid to model grid
+   !! History :  2.0  !  06-2006  (S. Masson, G. Madec) Original code
+   !!                 !  05-2008  (S. Alderson) Modified for Interpolation in memory
+   !!                 !                         from input grid to model grid
    !!----------------------------------------------------------------------
 
@@ -33,6 +33,6 @@
       CHARACTER(len = 34)  ::   wname       ! generic name of a NetCDF weights file to be used, blank if not
       CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation
-                                            ! a string starting with "U" or "V" for each component   
-                                            ! chars 2 onwards identify which components go together  
+      !                                     ! a string starting with "U" or "V" for each component   
+      !                                     ! chars 2 onwards identify which components go together  
    END TYPE FLD_N
 
@@ -51,5 +51,5 @@
       REAL(wp) , ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta       ! 2 consecutive record of input fields
       CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key
-                                                        ! into the WGTLIST structure
+      !                                                 ! into the WGTLIST structure
       CHARACTER(len = 34)             ::   vcomp        ! symbolic name for a vector component that needs rotation
       LOGICAL                         ::   rotn         ! flag to indicate whether field has been rotated
@@ -65,5 +65,5 @@
       INTEGER , DIMENSION(2)                  ::   ddims        ! shape of input grid
       INTEGER , DIMENSION(2)                  ::   botleft      ! top left corner of box in input grid containing 
-                                                                ! current processor grid
+      !                                                         ! current processor grid
       INTEGER , DIMENSION(2)                  ::   topright     ! top right corner of box 
       INTEGER                                 ::   jpiwgt       ! width of box on input grid
@@ -72,6 +72,6 @@
       INTEGER                                 ::   nestid       ! for agrif, keep track of nest we're in
       INTEGER                                 ::   overlap      ! =0 when cyclic grid has no overlapping EW columns
-                                                                ! =>1 when they have one or more overlapping columns      
-                                                                ! =-1 not cyclic
+      !                                                         ! =>1 when they have one or more overlapping columns      
+      !                                                         ! =-1 not cyclic
       LOGICAL                                 ::   cyclic       ! east-west cyclic or not
       INTEGER,  DIMENSION(:,:,:), POINTER     ::   data_jpi     ! array of source integers
@@ -93,7 +93,6 @@
    !! NEMO/OPA 3.3 , 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
 
@@ -259,6 +258,4 @@
       !! ** Purpose :  - if time interpolation, read before data 
       !!               - open current year file
-      !!
-      !! ** Method  :   
       !!----------------------------------------------------------------------
       INTEGER  , INTENT(in   ) ::   kn_fsbc   ! sbc computation period (in time step) 
@@ -408,6 +405,4 @@
       !!                  nrec_a(1): record number
       !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record (for print only)
-      !!
-      !! ** Method  :   
       !!----------------------------------------------------------------------
       INTEGER  , INTENT(in   )           ::   kn_fsbc   ! sbc computation period (in time step) 
@@ -555,11 +550,9 @@
       !!
       !! ** Purpose :   read the data
-      !!
-      !! ** Method  :   
-      !!----------------------------------------------------------------------
-      TYPE(FLD), INTENT(inout)   ::   sdjf   ! input field related variables
-      !!
-      INTEGER                    ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
-      INTEGER                    ::   iw     ! index into wgts array
+      !!----------------------------------------------------------------------
+      TYPE(FLD), INTENT(inout) ::   sdjf   ! input field related variables
+      !!
+      INTEGER                  ::   ipk    ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk )
+      INTEGER                  ::   iw     ! index into wgts array
       !!---------------------------------------------------------------------
             
@@ -593,9 +586,7 @@
       !!
       !! ** Purpose :   Vector fields may need to be rotated onto the local grid direction
-      !!
-      !! ** Method  :   
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY: wrk_use, wrk_release
-      USE wrk_nemo, ONLY: utmp => wrk_2d_4, vtmp => wrk_2d_5
+      USE wrk_nemo, ONLY: utmp => wrk_2d_4, vtmp => wrk_2d_5      ! 2D workspace
       !!
       INTEGER  , INTENT(in   )               ::   kt        ! ocean time step
@@ -609,7 +600,6 @@
       !!---------------------------------------------------------------------
 
-      IF(.not. wrk_use(2, 4,5))THEN
-         CALL ctl_stop('fld_rot: ERROR: requested workspace arrays are unavailable.')
-         RETURN
+      IF(.not. wrk_use(2, 4,5) ) THEN
+         CALL ctl_stop('fld_rot: ERROR: requested workspace arrays are unavailable.')   ;   RETURN
       END IF
 
@@ -646,9 +636,7 @@
           ENDIF
        END DO
-
-      IF(.not. wrk_release(2, 4,5))THEN
-         CALL ctl_stop('fld_rot: ERROR: failed to release workspace arrays.')
-      END IF
-
+      !
+      IF(.not. wrk_release(2, 4,5) )    CALL ctl_stop('fld_rot: ERROR: failed to release workspace arrays.')
+      !
    END SUBROUTINE fld_rot
 
@@ -659,12 +647,11 @@
       !!
       !! ** Purpose :   update the file name and open the file
-      !!
-      !! ** Method  :   
-      !!----------------------------------------------------------------------
-      TYPE(FLD), INTENT(inout)           ::   sdjf                  ! input field related variables
-      INTEGER  , INTENT(in   )           ::   kyear                 ! year value
-      INTEGER  , INTENT(in   )           ::   kmonth                ! month value
-      INTEGER  , INTENT(in   )           ::   kday                  ! day value
-      LOGICAL  , INTENT(in   ), OPTIONAL ::   ldstop                ! stop if open to read a non-existing file (default = .TRUE.)
+      !!----------------------------------------------------------------------
+      TYPE(FLD)        , INTENT(inout) ::   sdjf     ! input field related variables
+      INTEGER          , INTENT(in   ) ::   kyear    ! year value
+      INTEGER          , INTENT(in   ) ::   kmonth   ! month value
+      INTEGER          , INTENT(in   ) ::   kday     ! day value
+      LOGICAL, OPTIONAL, INTENT(in   ) ::   ldstop   ! stop if open to read a non-existing file (default = .TRUE.)
+      !!----------------------------------------------------------------------
 
       IF( sdjf%num /= 0 )   CALL iom_close( sdjf%num )   ! close file if already open
@@ -693,6 +680,4 @@
       !!
       !! ** Purpose :   fill sdf with sdf_n and control print
-      !!
-      !! ** Method  :   
       !!----------------------------------------------------------------------
       TYPE(FLD)  , DIMENSION(:), INTENT(inout) ::   sdf        ! structure of input fields (file informations, fields read)
@@ -748,13 +733,10 @@
       !!                if it is a new entry, the weights data is read in and
       !!                restructured (fld_weight)
-      !!
-      !! ** Method  :   
-      !!----------------------------------------------------------------------
-      TYPE( FLD ),      INTENT(in)    ::   sd        ! field with name of weights file
-      INTEGER,      INTENT(inout)     ::   kwgt      ! index of weights
-      !!
-      INTEGER                         ::   kw
-      INTEGER                         ::   nestid
-      LOGICAL                         ::   found
+      !!----------------------------------------------------------------------
+      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file
+      INTEGER    , INTENT(inout) ::   kwgt      ! index of weights
+      !!
+      INTEGER ::   kw, nestid   ! local integer
+      LOGICAL ::   found        ! local logical
       !!----------------------------------------------------------------------
       !
@@ -782,7 +764,8 @@
          CALL fld_weight( sd )
       ENDIF
-
+      !
    END SUBROUTINE wgt_list
 
+
    SUBROUTINE wgt_print( )
       !!---------------------------------------------------------------------
@@ -790,12 +773,8 @@
       !!
       !! ** Purpose :   print the list of known weights
-      !!
-      !! ** Method  :   
-      !!----------------------------------------------------------------------
-      !!
-      INTEGER                         ::   kw
-      !!----------------------------------------------------------------------
-      !
-
+      !!----------------------------------------------------------------------
+      INTEGER ::   kw   !
+      !!----------------------------------------------------------------------
+      !
       DO kw = 1, nxt_wgt-1
          WRITE(numout,*) 'weight file:  ',TRIM(ref_wgts(kw)%wgtname)
@@ -814,6 +793,7 @@
          IF( ASSOCIATED(ref_wgts(kw)%data_wgt) )  WRITE(numout,*) '       allocated'
       END DO
-
+      !
    END SUBROUTINE wgt_print
+
 
    SUBROUTINE fld_weight( sd )
@@ -823,6 +803,4 @@
       !! ** Purpose :   create a new WGT structure and fill in data from  
       !!                file, restructuring as required
-      !!
-      !! ** Method  :   
       !!----------------------------------------------------------------------
       USE wrk_nemo, ONLY: wrk_use, wrk_release, iwrk_use, iwrk_release
@@ -830,19 +808,18 @@
       USE wrk_nemo, ONLY: data_src => iwrk_2d_1
       !!
-      TYPE( FLD ),      INTENT(in)            ::   sd            ! field with name of weights file
-      !!
-      INTEGER                                 ::   jn            ! dummy loop indices
-      INTEGER                                 ::   inum          ! temporary logical unit
-      INTEGER                                 ::   id            ! temporary variable id
-      INTEGER                                 ::   ipk           ! temporary vertical dimension
-      CHARACTER (len=5)                       ::   aname
-      INTEGER , DIMENSION(3)                  ::   ddims
-      LOGICAL                                 ::   cyclical
-      INTEGER                                 ::   zwrap         ! temporary integer
-      !!----------------------------------------------------------------------
-      !
-      IF( (.NOT. wrk_use(2, 1)) .OR. (.NOT. iwrk_use(2,1)) )THEN
-         CALL ctl_stop('fld_weights: requested workspace arrays are unavailable.')
-         RETURN
+      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file
+      !!
+      INTEGER                ::   jn            ! dummy loop indices
+      INTEGER                ::   inum          ! temporary logical unit
+      INTEGER                ::   id            ! temporary variable id
+      INTEGER                ::   ipk           ! temporary vertical dimension
+      CHARACTER (len=5)      ::   aname
+      INTEGER , DIMENSION(3) ::   ddims
+      LOGICAL                ::   cyclical
+      INTEGER                ::   zwrap      ! local integer
+      !!----------------------------------------------------------------------
+      !
+      IF(  .NOT. wrk_use(2, 1)  .OR.  .NOT. iwrk_use(2,1)  ) THEN
+         CALL ctl_stop('fld_weights: requested workspace arrays are unavailable.')   ;   RETURN
       END IF
       !
@@ -957,11 +934,11 @@
       ENDIF
 
-      IF( (.NOT. wrk_release(2, 1)) .OR. (.NOT. iwrk_release(2,1)) )THEN
-         CALL ctl_stop('fld_weights: failed to release workspace arrays.')
-      END IF
-
+      IF( .NOT. wrk_release(2, 1) .OR.    &
+          .NOT.iwrk_release(2, 1)   )   CALL ctl_stop('fld_weights: failed to release workspace arrays')
+      !
    END SUBROUTINE fld_weight
 
-   SUBROUTINE fld_interp(num, clvar, kw, kk, dta, nrec)
+
+   SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec )
       !!---------------------------------------------------------------------
       !!                    ***  ROUTINE fld_interp  ***
@@ -969,23 +946,20 @@
       !! ** Purpose :   apply weights to input gridded data to create data
       !!                on model grid
-      !!
-      !! ** Method  :   
-      !!----------------------------------------------------------------------
-      INTEGER,          INTENT(in)                        ::   num                 ! stream number
-      CHARACTER(LEN=*), INTENT(in)                        ::   clvar               ! variable name
-      INTEGER,          INTENT(in)                        ::   kw                  ! weights number
-      INTEGER,          INTENT(in)                        ::   kk                  ! vertical dimension of kk
-      REAL(wp),         INTENT(inout), DIMENSION(:,:,:)   ::   dta              ! output field on model grid
-      INTEGER,          INTENT(in)                        ::   nrec                ! record number to read (ie time slice)
+      !!----------------------------------------------------------------------
+      INTEGER                   , INTENT(in   ) ::   num     ! stream number
+      CHARACTER(LEN=*)          , INTENT(in   ) ::   clvar   ! variable name
+      INTEGER                   , INTENT(in   ) ::   kw      ! weights number
+      INTEGER                   , INTENT(in   ) ::   kk      ! vertical dimension of kk
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   dta     ! output field on model grid
+      INTEGER                   , INTENT(in   ) ::   nrec    ! record number to read (ie time slice)
       !! 
-      INTEGER, DIMENSION(3)                               ::   rec1,recn           ! temporary arrays for start and length
-      INTEGER                                             ::  jk, jn, jm           ! loop counters
-      INTEGER                                             ::  ni, nj               ! lengths
-      INTEGER                                             ::  jpimin,jpiwid        ! temporary indices
-      INTEGER                                             ::  jpjmin,jpjwid        ! temporary indices
-      INTEGER                                             ::  jpi1,jpi2,jpj1,jpj2  ! temporary indices
-      !!----------------------------------------------------------------------
-      !
-
+      INTEGER, DIMENSION(3) ::   rec1,recn   ! temporary arrays for start and length
+      INTEGER ::  jk, jn, jm           ! loop counters
+      INTEGER ::  ni, nj               ! lengths
+      INTEGER ::  jpimin,jpiwid        ! temporary indices
+      INTEGER ::  jpjmin,jpjwid        ! temporary indices
+      INTEGER ::  jpi1,jpi2,jpj1,jpj2  ! temporary indices
+      !!----------------------------------------------------------------------
+      !
       !! for weighted interpolation we have weights at four corners of a box surrounding 
       !! a model grid point, each weight is multiplied by a grid value (bilinear case)
@@ -1107,19 +1081,19 @@
         END DO
 
-        ! gradient in the ij direction
-        DO jk = 1,4
-          DO jn = 1, jpj
-            DO jm = 1,jpi
-              ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
-              nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
-              dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
+         ! gradient in the ij direction
+         DO jk = 1,4
+            DO jn = 1, jpj
+               DO jm = 1,jpi
+                  ni = ref_wgts(kw)%data_jpi(jm,jn,jk)
+                  nj = ref_wgts(kw)%data_jpj(jm,jn,jk)
+                  dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( &
                                (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   &
                                (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:)))
+               END DO
             END DO
-          END DO
-        END DO
-
+         END DO
+         !
       END IF
-
+      !
    END SUBROUTINE fld_interp
 
@@ -1130,6 +1104,4 @@
       !!
       !! ** Purpose :  
-      !!
-      !! ** Method  :
       !!---------------------------------------------------------------------
       CHARACTER(len=*), INTENT(in)   ::   cdday   !3 first letters of the first day of the weekly file
@@ -1143,5 +1115,5 @@
       DO ijul = 1, 7
          IF( cl_week(ijul) == TRIM(cdday) ) EXIT
-      ENDDO
+      END DO
       IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) )
       !
@@ -1153,4 +1125,4 @@
    END FUNCTION ksec_week
 
-
+   !!======================================================================
 END MODULE fldread
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/geo2ocean.F90	(revision 2620)
@@ -4,81 +4,342 @@
    !! Ocean mesh    :  ???
    !!======================================================================
-   !! History :  OPA  !  07-1996  (O. Marti)  Original code
-   !!   NEMO     1.0  !  02-2008  (G. Madec)  F90: Free form
-   !!            3.0  !  
+   !! History :  OPA  ! 1989-07  (O. Marti)  Original code
+   !!   NEMO     1.0  ! 2002-08  (G. Madec)  F90: Free form + opt.
+   !!             -   ! 2006-02  (A. Caubel)  oce2geo - Original code
+   !!            2.0  ! 2007-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary
+   !!            3.0  ! 2008-07  (G. Madec)  geo2oce suppress lon/lat agruments
+   !!            3.3  ! 2010-10  (K. Mogensen) add obs_rot
+   !!            4.0  ! 2011-02  (G. Madec)  dynamical allocation + addition of angle_geo
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
-   !!   repcmo      : 
-   !!   angle       :
-   !!   geo2oce     :
-   !!   repere      :   old routine suppress it ???
+   !!   angle_msh_geo    : local sin/cos of the angle between model grid lines and the North direction
+   !!   angle_geo        : local sin/cos of the latitude and longitude at each mesh grid point
+   !!   geo2oce, oce2geo : ?
+   !!   obs_rot          : provide to OBS operator the sin &cos at u- and v-points
+   !!
+   !!   repcmo, repere, rot_rep :   old routines  ==> to be suppressed
    !!----------------------------------------------------------------------
-   USE dom_oce         ! mesh and scale factors
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE dom_oce        ! mesh and scale factors
+   USE phycst         ! physical constants
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   rot_rep, repcmo, repere, geo2oce, oce2geo   ! only rot_rep should be used
-                                             ! repcmo and repere are keep only for compatibility.
-                                             ! they are only a useless overlay of rot_rep
-
-   PUBLIC   obs_rot
-   PUBLIC   geo2oce_alloc ! Called in nemogcm.F90
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   &
-      gsint, gcost,   &  ! cos/sin between model grid lines and NP direction at T point
-      gsinu, gcosu,   &  ! cos/sin between model grid lines and NP direction at U point
-      gsinv, gcosv,   &  ! cos/sin between model grid lines and NP direction at V point
-      gsinf, gcosf       ! cos/sin between model grid lines and NP direction at F point
-
-   LOGICAL ::   lmust_init = .TRUE.        !: used to initialize the cos/sin variables (se above)
-
-   ! Local 'saved' arrays - one set for geo2oce and one set for oce2geo.
-   ! Declared here so can be allocated in ge2oce_alloc().
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zsinlon_o2g, zcoslon_o2g, zsinlat_o2g, zcoslat_o2g
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zsinlon_g2o, zcoslon_g2o, zsinlat_g2o, zcoslat_g2o
+   PUBLIC   geo2oce, oce2geo   !    
+   PUBLIC   obs_rot            !
+   !
+   PUBLIC   repcmo, repere, rot_rep    ! CAUTION: these routines are kept only for compatibility.
+   !                                   !          They are only a useless overlay of rot_rep
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsin, gcos         ! sinus & cosinus at T,U,V,F points
+   !                                                                     ! between mesh and NP direction
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gsinlon, gsinlat   ! sinus   of lon & lat at T,U,V,F points
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   gcoslon, gcoslat   ! cosinus of lon & lat at T,U,V,F points
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
-   FUNCTION geo2oce_alloc()
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE geo2oce_alloc  ***
-      !!----------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: geo2oce_alloc
-
-      ALLOCATE(gsint(jpi,jpj), gcost(jpi,jpj),   &
-               gsinu(jpi,jpj), gcosu(jpi,jpj),   &
-               gsinv(jpi,jpj), gcosv(jpi,jpj),   &
-               gsinf(jpi,jpj), gcosf(jpi,jpj),   &
-               !
-               zsinlon_o2g(jpi,jpj,4), zcoslon_o2g(jpi,jpj,4), & 
-               zsinlat_o2g(jpi,jpj,4), zcoslat_o2g(jpi,jpj,4), &
-               zsinlon_g2o(jpi,jpj,4), zcoslon_g2o(jpi,jpj,4), &
-               zsinlat_g2o(jpi,jpj,4), zcoslat_g2o(jpi,jpj,4), &
-               Stat=geo2oce_alloc)
-
-   END FUNCTION geo2oce_alloc
-
-
-   SUBROUTINE repcmo ( pxu1, pyu1, pxv1, pyv1,   &
-                       px2 , py2 )
+   SUBROUTINE angle_msh_geo
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE angle_msh_geo  ***
+      !! 
+      !! ** Purpose :   Compute angles between model grid lines and the North direction
+      !!
+      !! ** Action  :   allocate and compute (gsint, gcost, gsinu, gcosu, gsinv, 
+      !!              gcosv, gsinf, gcosf) arrays: sinus and cosinus of the angle 
+      !!              between the north-south axe and the j-direction at t, u, v and f-points
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      INTEGER  ::   ierr     ! local integer
+      REAL(wp) ::   zpi_4    ! local scalar (pi/4)
+      REAL(wp) ::   zlam, zphi, zlan, zphh   ! local scalars
+      REAL(wp) ::   zxnpt, zynpt, znnpt      ! x,y components & norm of the vector: T point to North Pole
+      REAL(wp) ::   zxnpu, zynpu, znnpu      ! x,y components & norm of the vector: U point to North Pole
+      REAL(wp) ::   zxnpv, zynpv, znnpv      ! x,y components & norm of the vector: V point to North Pole
+      REAL(wp) ::   zxnpf, zynpf, znnpf      ! x,y components & norm of the vector: F point to North Pole
+      REAL(wp) ::   zxvvt, zyvvt, znvvt      ! x,y components & norm of the vector: between V pts below and above a T pt
+      REAL(wp) ::   zxffu, zyffu, znffu      ! x,y components & norm of the vector: between F pts below and above a U pt
+      REAL(wp) ::   zxffv, zyffv, znffv      ! x,y components & norm of the vector: between F pts left  and right a V pt
+      REAL(wp) ::   zxuuf, zyuuf, znuuf      ! x,y components & norm of the vector: between U pts below and above a F pt
+      !!----------------------------------------------------------------------
+      !
+      ! already allocated & initialized  ==> return
+      ! -------------------------------
+      IF( ALLOCATED( gsin ) .AND. ALLOCATED( gcos )  )   RETURN
+
+      ! allocate cos & sin arrays 
+      ! -------------------------
+      ALLOCATE( gsin(jpi,jpj,4) , gcos(jpi,jpj,4) , STAT=ierr )
+      IF(lk_mpp)   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop('STOP', 'angle_msh_geo: unable to allocate arrays' )
+      !
+      ! initialize cos & sin arrays
+      ! ----------------------------
+      ! (computation done on the north stereographic polar plane)
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' angle_msh_geo : geographic <--> model mesh , cos/sin initialization '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~'
+
+      zpi_4 = rpi / 4._wp
+      !
+      DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+         DO ji = fs_2, jpi   ! vector opt.
+            !
+            zlam = glamt(ji,jj)            ! north pole direction & modulous (at t-point)
+            zphi = gphit(ji,jj)
+            zxnpt = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            zynpt = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            znnpt = zxnpt*zxnpt + zynpt*zynpt
+
+            zlam = glamu(ji,jj)            ! north pole direction & modulous (at u-point)
+            zphi = gphiu(ji,jj)
+            zxnpu = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            zynpu = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            znnpu = zxnpu*zxnpu + zynpu*zynpu
+
+            zlam = glamv(ji,jj)            ! north pole direction & modulous (at v-point)
+            zphi = gphiv(ji,jj)
+            zxnpv = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            zynpv = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            znnpv = zxnpv*zxnpv + zynpv*zynpv
+            
+            zlam = glamf(ji,jj)            ! north pole direction & modulous (at f-point)
+            zphi = gphif(ji,jj)
+            zxnpf = 0._wp - 2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            zynpf = 0._wp - 2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )
+            znnpf = zxnpf*zxnpf + zynpf*zynpf
+
+            zlam = glamv(ji,jj  )            ! j-direction: v-point segment direction (around t-point)
+            zphi = gphiv(ji,jj  )
+            zlan = glamv(ji,jj-1)
+            zphh = gphiv(ji,jj-1)
+            zxvvt =  2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            zyvvt =  2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
+            znvvt = MAX( znvvt, 1.e-14 )
+
+            zlam = glamf(ji,jj  )            ! j-direction: f-point segment direction (around u-point)
+            zphi = gphif(ji,jj  )
+            zlan = glamf(ji,jj-1)
+            zphh = gphif(ji,jj-1)
+            zxffu =  2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            zyffu =  2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
+            znffu = MAX( znffu, 1.e-14 )
+
+            zlam = glamf(ji  ,jj)            ! i-direction: f-point segment direction (around v-point)
+            zphi = gphif(ji  ,jj)
+            zlan = glamf(ji-1,jj)
+            zphh = gphif(ji-1,jj)
+            zxffv =  2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            zyffv =  2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
+            znffv = MAX( znffv, 1.e-14 )
+
+            zlam = glamu(ji,jj+1)            ! j-direction: u-point segment direction (around f-point)
+            zphi = gphiu(ji,jj+1)
+            zlan = glamu(ji,jj  )
+            zphh = gphiu(ji,jj  )
+            zxuuf =  2._wp * COS( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * COS( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            zyuuf =  2._wp * SIN( rad*zlam ) * TAN( zpi_4 - rad*zphi*0.5 )   &
+               &  -  2._wp * SIN( rad*zlan ) * TAN( zpi_4 - rad*zphh*0.5 )
+            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
+            znuuf = MAX( znuuf, 1.e-14 )
+
+            ! cosinus and sinus using scalar and vectorial products
+            gsin(ji,jj,1) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt      ! T-point
+            gcos(ji,jj,1) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
+
+            gsin(ji,jj,2) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu      ! U-point
+            gcos(ji,jj,2) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
+            !
+            gsin(ji,jj,3) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv        ! F-point
+            gcos(ji,jj,3) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv            ! (caution, rotation of 90 degres)
+            !
+            gsin(ji,jj,4) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf      ! V-point
+            gcos(ji,jj,4) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
+            !
+         END DO
+      END DO
+      !
+      ! Geographic mesh case
+      ! --------------------
+      DO jj = 2, jpjm1
+         DO ji = fs_2, jpi   ! vector opt.
+            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
+               gsin(ji,jj,1) = 0.
+               gcos(ji,jj,1) = 1.
+            ENDIF
+            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
+               gsin(ji,jj,2) = 0.
+               gcos(ji,jj,2) = 1.
+            ENDIF
+            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN
+               gsin(ji,jj,3) = 0.
+               gcos(ji,jj,3) = 1.
+            ENDIF
+            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
+               gsin(ji,jj,4) = 0.
+               gcos(ji,jj,4) = 1.
+            ENDIF
+         END DO
+      END DO
+      !
+      CALL lbc_lnk( gcos(:,:,1), 'T', -1. )   ;   CALL lbc_lnk( gsin(:,:,1), 'T', -1. )      ! lateral boundary cond.
+      CALL lbc_lnk( gcos(:,:,2), 'U', -1. )   ;   CALL lbc_lnk( gsin(:,:,2), 'U', -1. )
+      CALL lbc_lnk( gcos(:,:,3), 'V', -1. )   ;   CALL lbc_lnk( gsin(:,:,3), 'V', -1. )
+      CALL lbc_lnk( gcos(:,:,4), 'F', -1. )   ;   CALL lbc_lnk( gsin(:,:,4), 'F', -1. )
+      !
+   END SUBROUTINE angle_msh_geo
+
+
+   SUBROUTINE angle_geo
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE angle_geo  ***
+      !!      
+      !! ** Purpose :   compute one for all, and at each mesh grid points, 
+      !!              the local cos/sin of latitude/longitude.
+      !!----------------------------------------------------------------------
+      INTEGER ierr
+      !!----------------------------------------------------------------------
+      !
+      IF( ALLOCATED( gsinlon ) .AND. ALLOCATED( gcoslon ) .AND.   &   !==  already allocated & initialized
+          ALLOCATED( gsinlat ) .AND. ALLOCATED( gcoslat )         )   RETURN
+
+      ALLOCATE( gsinlon(jpi,jpj,4) , gcoslon(jpi,jpj,4) ,     &       !==  allocate the arrays
+         &      gsinlat(jpi,jpj,4) , gcoslat(jpi,jpj,4) , STAT=ierr )
+      IF( lk_mpp    )   CALL mpp_sum( ierr )
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'angle_geo: unable to allocate arrays')
+         !
+      gsinlon(:,:,1) = SIN( rad * glamt(:,:) )     ! T-point
+      gcoslon(:,:,1) = COS( rad * glamt(:,:) )
+      gsinlat(:,:,1) = SIN( rad * gphit(:,:) )
+      gcoslat(:,:,1) = COS( rad * gphit(:,:) )
+      !
+      gsinlon(:,:,2) = SIN( rad * glamu(:,:) )     ! U-point
+      gcoslon(:,:,2) = COS( rad * glamu(:,:) )
+      gsinlat(:,:,2) = SIN( rad * gphiu(:,:) )
+      gcoslat(:,:,2) = COS( rad * gphiu(:,:) )
+      !
+      gsinlon(:,:,3) = SIN( rad * glamv(:,:) )     ! V-point
+      gcoslon(:,:,3) = COS( rad * glamv(:,:) )
+      gsinlat(:,:,3) = SIN( rad * gphiv(:,:) )
+      gcoslat(:,:,3) = COS( rad * gphiv(:,:) )
+      !
+      gsinlon(:,:,4) = SIN( rad * glamf(:,:) )     ! T-point
+      gcoslon(:,:,4) = COS( rad * glamf(:,:) )
+      gsinlat(:,:,4) = SIN( rad * gphif(:,:) )
+      gcoslat(:,:,4) = COS( rad * gphif(:,:) )
+      !
+   END SUBROUTINE angle_geo
+
+
+   SUBROUTINE geo2oce( pxx, pyy, pzz, cgrid, pte, ptn )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE geo2oce  ***
+      !!      
+      !! ** Purpose :   Change a vector from geocentric to east/north
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxx, pyy, pzz   !
+      CHARACTER(len=1)            , INTENT(in   ) ::   cgrid           ! type of grid point
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pte, ptn        !
+      !!
+      INTEGER ::   ipt     ! local integer
+      !!----------------------------------------------------------------------
+
+      CALL angle_geo             ! initialization of coc & sin (just a return if not the 1st call)
+      !
+      SELECT CASE (cgrid)        ! type of point
+      CASE ('T')   ;   ipt = 1
+      CASE ('U')   ;   ipt = 2
+      CASE ('V')   ;   ipt = 3
+      CASE ('F')   ;   ipt = 4
+      CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
+      END SELECT
+  
+      !                          ! transformation
+      pte(:,:) = - gsinlon(:,:,ipt)                    * pxx(:,:) + gcoslon(:,:,ipt) * pyy(:,:)
+      ptn(:,:) = - gsinlat(:,:,ipt) * gsinlat(:,:,ipt) * pxx(:,:)    &
+         &       - gsinlon(:,:,ipt) * gsinlat(:,:,ipt) * pyy(:,:) + gcoslat(:,:,ipt) * pzz(:,:)
+!!$   ptv(:,:) =   gsinlat(:,:,ipt) * gcoslat(:,:,ipt) * pxx(:,:)    &
+!!$              + zsinlon(:,:,ipt) * gcoslat(:,:,ipt) * pyy(:,:) + gsinlat(:,:,ipt) * pzz(:,:)
+      !
+   END SUBROUTINE geo2oce
+
+
+   SUBROUTINE oce2geo ( pte, ptn, cgrid, pxx , pyy , pzz )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE oce2geo  ***
+      !!      
+      !! ** Purpose :   Change vector from east/north to geocentric
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::  pte, ptn          ! 
+      CHARACTER(len=1)        , INTENT(in   ) ::  cgrid             !
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::  pxx , pyy , pzz   !
+      !!
+      INTEGER ::   ipt     !
+      !!----------------------------------------------------------------------
+      !
+      CALL angle_geo             ! initialization of coc & sin (just a return if not the 1st call)
+      !
+      SELECT CASE (cgrid)        ! type of point
+      CASE ('T')   ;   ipt = 1
+      CASE ('U')   ;   ipt = 2
+      CASE ('V')   ;   ipt = 3
+      CASE ('F')   ;   ipt = 4
+      CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
+      END SELECT
+      !
+      !                          ! transformation
+      pxx(:,:) = - gsinlon(:,:,ipt) * pte(:,:) - gcoslon(:,:,ipt) * gsinlat(:,:,ipt) * ptn(:,:) 
+      pyy(:,:) =   gcoslon(:,:,ipt) * pte(:,:) - gsinlon(:,:,ipt) * gsinlat(:,:,ipt) * ptn(:,:)
+      pzz(:,:) =   gcoslat(:,:,ipt) * ptn(:,:)
+      !
+   END SUBROUTINE oce2geo
+
+
+   SUBROUTINE obs_rot( psinu, pcosu, psinv, pcosv )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE obs_rot  ***
+      !!
+      !! ** Purpose :   provide to OBS operator the sin &cos at u- and v-points
+      !!              in order to rotate currents at observation points
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(out)::   psinu, pcosu, psinv, pcosv   ! copy of data
+      !!----------------------------------------------------------------------
+      !
+      CALL angle_msh_geo      ! initialization of coc & sin (just a return if not the 1st call)
+      !
+      psinu(:,:) = gsin(:,:,2)      ! U-point
+      pcosu(:,:) = gcos(:,:,2)
+      psinv(:,:) = gsin(:,:,3)      ! V-point
+      pcosv(:,:) = gcos(:,:,3)
+      !
+   END SUBROUTINE obs_rot
+
+
+   SUBROUTINE repcmo( pxu1, pyu1, pxv1, pyv1, px2 , py2 )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE repcmo  ***
       !!
       !! ** Purpose :   Change vector componantes from a geographic grid to a
-      !!      stretched coordinates grid.
+      !!              stretched coordinates grid.
       !!
       !! ** Method  :   Initialization of arrays at the first call.
@@ -92,446 +353,27 @@
       REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   py2          ! j-componante (defined at v-point)
       !!----------------------------------------------------------------------
-      
-      ! Change from geographic to stretched coordinate
-      ! ----------------------------------------------
-      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )
+      !      
+      CALL rot_rep( pxu1, pyu1, 'U', 'en->i',px2 )    ! Change from geographic to stretched coordinate
       CALL rot_rep( pxv1, pyv1, 'V', 'en->j',py2 )
-      
+      !
    END SUBROUTINE repcmo
 
 
-   SUBROUTINE rot_rep ( pxin, pyin, cd_type, cdtodo, prot )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE rot_rep  ***
-      !!
-      !! ** Purpose :   Rotate the Repere: Change vector componantes between
-      !!                geographic grid <--> stretched coordinates grid.
-      !!
-      !! History :
-      !!   9.2  !  07-04  (S. Masson)  
-      !!                  (O. Marti ) Original code (repere and repcmo)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( IN ) ::   pxin, pyin   ! vector componantes
-      CHARACTER(len=1),             INTENT( IN ) ::   cd_type      ! define the nature of pt2d array grid-points
-      CHARACTER(len=5),             INTENT( IN ) ::   cdtodo       ! specify the work to do:
-      !!                                                           ! 'en->i' east-north componantes to model i componante
-      !!                                                           ! 'en->j' east-north componantes to model j componante
-      !!                                                           ! 'ij->e' model i-j componantes to east componante
-      !!                                                           ! 'ij->n' model i-j componantes to east componante
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) ::   prot      
-
-      !!----------------------------------------------------------------------
-
-      ! Initialization of gsin* and gcos* at first call
-      ! -----------------------------------------------
-
-      IF( lmust_init ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' rot_rep : geographic <--> stretched'
-         IF(lwp) WRITE(numout,*) ' ~~~~~    coordinate transformation'
-
-         CALL angle       ! initialization of the transformation
-         lmust_init = .FALSE.
-
-      ENDIF
-      
-      SELECT CASE (cdtodo)
-      CASE ('en->i')      ! 'en->i' est-north componantes to model i componante
-         SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) + pyin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) + pyin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) + pyin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) + pyin(:,:) * gsinf(:,:)
-         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
-         END SELECT
-      CASE ('en->j')      ! 'en->j' est-north componantes to model j componante
-         SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) - pxin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) - pxin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) - pxin(:,:) * gsinv(:,:)   
-         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) - pxin(:,:) * gsinf(:,:)   
-         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
-         END SELECT
-      CASE ('ij->e')      ! 'ij->e' model i-j componantes to est componante
-         SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pxin(:,:) * gcost(:,:) - pyin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pxin(:,:) * gcosu(:,:) - pyin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pxin(:,:) * gcosv(:,:) - pyin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pxin(:,:) * gcosf(:,:) - pyin(:,:) * gsinf(:,:)
-         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
-         END SELECT
-      CASE ('ij->n')      ! 'ij->n' model i-j componantes to est componante
-         SELECT CASE (cd_type)
-         CASE ('T')   ;   prot(:,:) = pyin(:,:) * gcost(:,:) + pxin(:,:) * gsint(:,:)
-         CASE ('U')   ;   prot(:,:) = pyin(:,:) * gcosu(:,:) + pxin(:,:) * gsinu(:,:)
-         CASE ('V')   ;   prot(:,:) = pyin(:,:) * gcosv(:,:) + pxin(:,:) * gsinv(:,:)
-         CASE ('F')   ;   prot(:,:) = pyin(:,:) * gcosf(:,:) + pxin(:,:) * gsinf(:,:)
-         CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
-         END SELECT
-      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
-      END SELECT
-      
-   END SUBROUTINE rot_rep
-
-
-   SUBROUTINE angle
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE angle  ***
-      !! 
-      !! ** Purpose :   Compute angles between model grid lines and the North direction
-      !!
-      !! ** Method  :
-      !!
-      !! ** Action  :   Compute (gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf) arrays:
-      !!      sinus and cosinus of the angle between the north-south axe and the 
-      !!      j-direction at t, u, v and f-points
-      !!
-      !! History :
-      !!   7.0  !  96-07  (O. Marti )  Original code
-      !!   8.0  !  98-06  (G. Madec )
-      !!   8.5  !  98-06  (G. Madec )  Free form, F90 + opt.
-      !!   9.2  !  07-04  (S. Masson)  Add T, F points and bugfix in cos lateral boundary
-      !!----------------------------------------------------------------------
-      INTEGER ::   ji, jj      ! dummy loop indices
-      !!
-      REAL(wp) ::   &
-         zlam, zphi,            &  ! temporary scalars
-         zlan, zphh,            &  !    "         "
-         zxnpt, zynpt, znnpt,   &  ! x,y components and norm of the vector: T point to North Pole
-         zxnpu, zynpu, znnpu,   &  ! x,y components and norm of the vector: U point to North Pole
-         zxnpv, zynpv, znnpv,   &  ! x,y components and norm of the vector: V point to North Pole
-         zxnpf, zynpf, znnpf,   &  ! x,y components and norm of the vector: F point to North Pole
-         zxvvt, zyvvt, znvvt,   &  ! x,y components and norm of the vector: between V points below and above a T point
-         zxffu, zyffu, znffu,   &  ! x,y components and norm of the vector: between F points below and above a U point
-         zxffv, zyffv, znffv,   &  ! x,y components and norm of the vector: between F points left  and right a V point
-         zxuuf, zyuuf, znuuf       ! x,y components and norm of the vector: between U points below and above a F point
-      !!----------------------------------------------------------------------
-
-      ! ============================= !
-      ! Compute the cosinus and sinus !
-      ! ============================= !
-      ! (computation done on the north stereographic polar plane)
-
-      DO jj = 2, jpjm1
-!CDIR NOVERRCHK
-         DO ji = fs_2, jpi   ! vector opt.
-
-            ! north pole direction & modulous (at t-point)
-            zlam = glamt(ji,jj)
-            zphi = gphit(ji,jj)
-            zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            znnpt = zxnpt*zxnpt + zynpt*zynpt
-
-            ! north pole direction & modulous (at u-point)
-            zlam = glamu(ji,jj)
-            zphi = gphiu(ji,jj)
-            zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            znnpu = zxnpu*zxnpu + zynpu*zynpu
-
-            ! north pole direction & modulous (at v-point)
-            zlam = glamv(ji,jj)
-            zphi = gphiv(ji,jj)
-            zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            znnpv = zxnpv*zxnpv + zynpv*zynpv
-
-            ! north pole direction & modulous (at f-point)
-            zlam = glamf(ji,jj)
-            zphi = gphif(ji,jj)
-            zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
-            znnpf = zxnpf*zxnpf + zynpf*zynpf
-
-            ! j-direction: v-point segment direction (around t-point)
-            zlam = glamv(ji,jj  )
-            zphi = gphiv(ji,jj  )
-            zlan = glamv(ji,jj-1)
-            zphh = gphiv(ji,jj-1)
-            zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            zyvvt =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            znvvt = SQRT( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  )
-            znvvt = MAX( znvvt, 1.e-14 )
-
-            ! j-direction: f-point segment direction (around u-point)
-            zlam = glamf(ji,jj  )
-            zphi = gphif(ji,jj  )
-            zlan = glamf(ji,jj-1)
-            zphh = gphif(ji,jj-1)
-            zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            zyffu =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            znffu = SQRT( znnpu * ( zxffu*zxffu + zyffu*zyffu )  )
-            znffu = MAX( znffu, 1.e-14 )
-
-            ! i-direction: f-point segment direction (around v-point)
-            zlam = glamf(ji  ,jj)
-            zphi = gphif(ji  ,jj)
-            zlan = glamf(ji-1,jj)
-            zphh = gphif(ji-1,jj)
-            zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            zyffv =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            znffv = SQRT( znnpv * ( zxffv*zxffv + zyffv*zyffv )  )
-            znffv = MAX( znffv, 1.e-14 )
-
-            ! j-direction: u-point segment direction (around f-point)
-            zlam = glamu(ji,jj+1)
-            zphi = gphiu(ji,jj+1)
-            zlan = glamu(ji,jj  )
-            zphh = gphiu(ji,jj  )
-            zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            zyuuf =  2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
-               &  -  2. * SIN( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
-            znuuf = SQRT( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  )
-            znuuf = MAX( znuuf, 1.e-14 )
-
-            ! cosinus and sinus using scalar and vectorial products
-            gsint(ji,jj) = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt
-            gcost(ji,jj) = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt
-
-            gsinu(ji,jj) = ( zxnpu*zyffu - zynpu*zxffu ) / znffu
-            gcosu(ji,jj) = ( zxnpu*zxffu + zynpu*zyffu ) / znffu
-
-            gsinf(ji,jj) = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf
-            gcosf(ji,jj) = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf
-
-            ! (caution, rotation of 90 degres)
-            gsinv(ji,jj) = ( zxnpv*zxffv + zynpv*zyffv ) / znffv
-            gcosv(ji,jj) =-( zxnpv*zyffv - zynpv*zxffv ) / znffv
-
-         END DO
-      END DO
-
-      ! =============== !
-      ! Geographic mesh !
-      ! =============== !
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, jpi   ! vector opt.
-            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
-               gsint(ji,jj) = 0.
-               gcost(ji,jj) = 1.
-            ENDIF
-            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
-               gsinu(ji,jj) = 0.
-               gcosu(ji,jj) = 1.
-            ENDIF
-            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN
-               gsinv(ji,jj) = 0.
-               gcosv(ji,jj) = 1.
-            ENDIF
-            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
-               gsinf(ji,jj) = 0.
-               gcosf(ji,jj) = 1.
-            ENDIF
-         END DO
-      END DO
-
-      ! =========================== !
-      ! Lateral boundary conditions !
-      ! =========================== !
-
-      ! lateral boundary cond.: T-, U-, V-, F-pts, sgn
-      CALL lbc_lnk( gcost, 'T', -1. )   ;   CALL lbc_lnk( gsint, 'T', -1. )
-      CALL lbc_lnk( gcosu, 'U', -1. )   ;   CALL lbc_lnk( gsinu, 'U', -1. )
-      CALL lbc_lnk( gcosv, 'V', -1. )   ;   CALL lbc_lnk( gsinv, 'V', -1. )
-      CALL lbc_lnk( gcosf, 'F', -1. )   ;   CALL lbc_lnk( gsinf, 'F', -1. )
-
-   END SUBROUTINE angle
-
-
-   SUBROUTINE geo2oce ( pxx, pyy, pzz, cgrid,     &
-                        pte, ptn )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE geo2oce  ***
-      !!      
-      !! ** Purpose :
-      !!
-      !! ** Method  :   Change wind stress from geocentric to east/north
-      !!
-      !! History :
-      !!        !         (O. Marti)  Original code
-      !!        !  91-03  (G. Madec)
-      !!        !  92-07  (M. Imbard)
-      !!        !  99-11  (M. Imbard) NetCDF format with IOIPSL
-      !!        !  00-08  (D. Ludicone) Reduced section at Bab el Mandeb
-      !!   8.5  !  02-06  (G. Madec)  F90: Free form
-      !!   3.0  !  07-08  (G. Madec)  geo2oce suppress lon/lat agruments
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::  pxx, pyy, pzz
-      CHARACTER(len=1)            , INTENT(in   ) ::  cgrid
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::  pte, ptn
-      !!
-      REAL(wp), PARAMETER :: rpi = 3.141592653E0
-      REAL(wp), PARAMETER :: rad = rpi / 180.e0
-      INTEGER ::   ig     !
-      !! * Local save
-      LOGICAL , SAVE, DIMENSION(4)         ::   linit = .FALSE.
-      !!----------------------------------------------------------------------
-
-      SELECT CASE( cgrid)
-         CASE ( 'T' )   
-            ig = 1
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_g2o(:,:,ig) = SIN( rad * glamt(:,:) )
-               zcoslon_g2o(:,:,ig) = COS( rad * glamt(:,:) )
-               zsinlat_g2o(:,:,ig) = SIN( rad * gphit(:,:) )
-               zcoslat_g2o(:,:,ig) = COS( rad * gphit(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'U' )   
-            ig = 2
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_g2o(:,:,ig) = SIN( rad * glamu(:,:) )
-               zcoslon_g2o(:,:,ig) = COS( rad * glamu(:,:) )
-               zsinlat_g2o(:,:,ig) = SIN( rad * gphiu(:,:) )
-               zcoslat_g2o(:,:,ig) = COS( rad * gphiu(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'V' )   
-            ig = 3
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_g2o(:,:,ig) = SIN( rad * glamv(:,:) )
-               zcoslon_g2o(:,:,ig) = COS( rad * glamv(:,:) )
-               zsinlat_g2o(:,:,ig) = SIN( rad * gphiv(:,:) )
-               zcoslat_g2o(:,:,ig) = COS( rad * gphiv(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'F' )   
-            ig = 4
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_g2o(:,:,ig) = SIN( rad * glamf(:,:) )
-               zcoslon_g2o(:,:,ig) = COS( rad * glamf(:,:) )
-               zsinlat_g2o(:,:,ig) = SIN( rad * gphif(:,:) )
-               zcoslat_g2o(:,:,ig) = COS( rad * gphif(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE default   
-            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
-            CALL ctl_stop( ctmp1 )
-      END SELECT
-      
-      pte = - zsinlon_g2o(:,:,ig) * pxx + zcoslon_g2o(:,:,ig) * pyy
-      ptn = - zcoslon_g2o(:,:,ig) * zsinlat_g2o(:,:,ig) * pxx    &
-            - zsinlon_g2o(:,:,ig) * zsinlat_g2o(:,:,ig) * pyy    &
-            + zcoslat_g2o(:,:,ig) * pzz
-!!$   ptv =   zcoslon(:,:,ig) * zcoslat(:,:,ig) * pxx    &
-!!$         + zsinlon(:,:,ig) * zcoslat(:,:,ig) * pyy    &
-!!$         + zsinlat(:,:,ig) * pzz
-      !
-   END SUBROUTINE geo2oce
-
-   SUBROUTINE oce2geo ( pte, ptn, cgrid,     &
-                        pxx , pyy , pzz )
-      !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE oce2geo  ***
-      !!      
-      !! ** Purpose :
-      !!
-      !! ** Method  :   Change vector from east/north to geocentric
-      !!
-      !! History :
-      !!        !         (A. Caubel)  oce2geo - Original code
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:,:), INTENT( IN    ) ::  pte, ptn
-      CHARACTER(len=1)        , INTENT( IN    ) ::  cgrid
-      REAL(wp), DIMENSION(:,:), INTENT(   OUT ) ::  pxx , pyy , pzz
-      !!
-      REAL(wp), PARAMETER :: rpi = 3.141592653E0
-      REAL(wp), PARAMETER :: rad = rpi / 180.e0
-      INTEGER ::   ig     !
-      !! * Local save
-      LOGICAL , SAVE, DIMENSION(4)         ::   linit = .FALSE.
-      !!----------------------------------------------------------------------
-
-      SELECT CASE( cgrid)
-         CASE ( 'T' )   
-            ig = 1
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_o2g(:,:,ig) = SIN( rad * glamt(:,:) )
-               zcoslon_o2g(:,:,ig) = COS( rad * glamt(:,:) )
-               zsinlat_o2g(:,:,ig) = SIN( rad * gphit(:,:) )
-               zcoslat_o2g(:,:,ig) = COS( rad * gphit(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'U' )   
-            ig = 2
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_o2g(:,:,ig) = SIN( rad * glamu(:,:) )
-               zcoslon_o2g(:,:,ig) = COS( rad * glamu(:,:) )
-               zsinlat_o2g(:,:,ig) = SIN( rad * gphiu(:,:) )
-               zcoslat_o2g(:,:,ig) = COS( rad * gphiu(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'V' )   
-            ig = 3
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_o2g(:,:,ig) = SIN( rad * glamv(:,:) )
-               zcoslon_o2g(:,:,ig) = COS( rad * glamv(:,:) )
-               zsinlat_o2g(:,:,ig) = SIN( rad * gphiv(:,:) )
-               zcoslat_o2g(:,:,ig) = COS( rad * gphiv(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE ( 'F' )   
-            ig = 4
-            IF( .NOT. linit(ig) ) THEN 
-               zsinlon_o2g(:,:,ig) = SIN( rad * glamf(:,:) )
-               zcoslon_o2g(:,:,ig) = COS( rad * glamf(:,:) )
-               zsinlat_o2g(:,:,ig) = SIN( rad * gphif(:,:) )
-               zcoslat_o2g(:,:,ig) = COS( rad * gphif(:,:) )
-               linit(ig) = .TRUE.
-            ENDIF
-         CASE default   
-            WRITE(ctmp1,*) 'geo2oce : bad grid argument : ', cgrid
-            CALL ctl_stop( ctmp1 )
-      END SELECT
-
-       pxx = - zsinlon_o2g(:,:,ig) * pte - zcoslon_o2g(:,:,ig) * zsinlat_o2g(:,:,ig) * ptn 
-       pyy =   zcoslon_o2g(:,:,ig) * pte - zsinlon_o2g(:,:,ig) * zsinlat_o2g(:,:,ig) * ptn
-       pzz =   zcoslat_o2g(:,:,ig) * ptn
-
-      
-   END SUBROUTINE oce2geo
-
-
-   SUBROUTINE repere ( px1, py1, px2, py2, kchoix, cd_type )
+   SUBROUTINE repere( px1, py1, px2, py2, kchoix, cd_type )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE repere  ***
       !!        
       !! ** Purpose :   Change vector componantes between a geopgraphic grid 
-      !!      and a stretched coordinates grid.
-      !!
-      !! ** Method  :   
-      !!
-      !! ** Action  :
-      !!
-      !! History :
-      !!        !  89-03  (O. Marti)  original code
-      !!        !  92-02  (M. Imbard)
-      !!        !  93-03  (M. Guyon)  symetrical conditions
-      !!        !  98-05  (B. Blanke)
-      !!   8.5  !  02-08  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), INTENT( IN   ), DIMENSION(:,:) ::   &
-         px1, py1          ! two horizontal components to be rotated
-      REAL(wp), INTENT( OUT  ), DIMENSION(:,:) ::   &
-         px2, py2          ! the two horizontal components in the model repere
-      INTEGER, INTENT( IN ) ::   &
-         kchoix   ! type of transformation
-                  ! = 1 change from geographic to model grid.
-                  ! =-1 change from model to geographic grid
-      CHARACTER(len=1), INTENT( IN ), OPTIONAL ::   cd_type      ! define the nature of pt2d array grid-points
+      !!              and a stretched coordinates grid.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   px1, py1   ! the 2 horizontal components to be rotated
+      REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   px2, py2   ! the 2 horizontal components in the model repere
+      INTEGER                   , INTENT(in   ) ::   kchoix     ! = 1 change from geographic to model grid.
+      !                                                         ! =-1 change from model to geographic grid
+      CHARACTER(len=1), OPTIONAL, INTENT(in   ) ::   cd_type    ! define the nature of pt2d array grid-points
       !
       CHARACTER(len=1) ::   cl_type      ! define the nature of pt2d array grid-points (T point by default)
       !!----------------------------------------------------------------------
-
+      !
       cl_type = 'T'
       IF( PRESENT(cd_type) )   cl_type = cd_type
@@ -546,44 +388,46 @@
       CASE DEFAULT   ;   CALL ctl_stop( 'repere: Syntax Error in the definition of kchoix (1 OR -1' )
       END SELECT
+      !
+   END SUBROUTINE repere
+
+
+   SUBROUTINE rot_rep( pxin, pyin, cd_type, cdtodo, prot )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE rot_rep  ***
+      !!
+      !! ** Purpose :   Rotate the Repere: Change vector componantes between
+      !!                geographic grid <--> stretched coordinates grid.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pxin, pyin   ! vector componantes
+      CHARACTER(len=1),             INTENT(in   ) ::   cd_type      ! define the nature of pt2d array grid-points
+      CHARACTER(len=5),             INTENT(in   ) ::   cdtodo       ! specify the work to do:
+      !!                                                            ! 'en->i' east-north componantes to model i-componante
+      !!                                                            ! 'en->j' east-north componantes to model j-componante
+      !!                                                            ! 'ij->e' model i-j componantes to east componante
+      !!                                                            ! 'ij->n' model i-j componantes to east componante
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   prot         !
+      !
+      INTEGER ::   ipt   !
+      !!----------------------------------------------------------------------
+
+      CALL angle_msh_geo   ! initialization of the transformation (just a return if not the 1st call)
       
-   END SUBROUTINE repere
-
-
-   SUBROUTINE obs_rot ( psinu, pcosu, psinv, pcosv )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE obs_rot  ***
-      !!
-      !! ** Purpose :   Copy gsinu, gcosu, gsinv and gsinv
-      !!                to input data for rotations of
-      !!                current at observation points
-      !!
-      !! History :
-      !!   9.2  !  09-02  (K. Mogensen)
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( OUT )::   &
-         & psinu, pcosu, psinv, pcosv! copy of data
-
-      !!----------------------------------------------------------------------
-
-      ! Initialization of gsin* and gcos* at first call
-      ! -----------------------------------------------
-
-      IF( lmust_init ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
-         IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation'
-
-         CALL angle       ! initialization of the transformation
-         lmust_init = .FALSE.
-
-      ENDIF
-
-      psinu(:,:) = gsinu(:,:)
-      pcosu(:,:) = gcosu(:,:)
-      psinv(:,:) = gsinv(:,:)
-      pcosv(:,:) = gcosv(:,:)
-
-   END SUBROUTINE obs_rot
-
+      SELECT CASE (cd_type)
+      CASE ('T')   ;   ipt = 1
+      CASE ('U')   ;   ipt = 2
+      CASE ('V')   ;   ipt = 3
+      CASE ('F')   ;   ipt = 4
+      CASE DEFAULT   ;   CALL ctl_stop( 'Only T, U, V and F grid points are coded' )
+      END SELECT
+      
+      SELECT CASE (cdtodo)
+      CASE ('en->i')   ;   prot(:,:) = pxin(:,:) * gcos(:,:,ipt) + pyin(:,:) * gsin(:,:,ipt)    ! east-north to model i-component
+      CASE ('en->j')   ;   prot(:,:) = pyin(:,:) * gcos(:,:,ipt) - pxin(:,:) * gsin(:,:,ipt)    ! east-north to model j-component
+      CASE ('ij->e')   ;   prot(:,:) = pxin(:,:) * gcos(:,:,ipt) - pyin(:,:) * gsin(:,:,ipt)    ! model i-j  to east    component
+      CASE ('ij->n')   ;   prot(:,:) = pyin(:,:) * gcos(:,:,ipt) + pxin(:,:) * gsin(:,:,ipt)    ! model i-j  to north   component
+      CASE DEFAULT   ;   CALL ctl_stop( 'rot_rep: Syntax Error in the definition of cdtodo' )
+      END SELECT
+      
+   END SUBROUTINE rot_rep
 
   !!======================================================================
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/oasis4_date.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/oasis4_date.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/oasis4_date.F90	(revision 2620)
@@ -1,33 +1,41 @@
 MODULE oasis4_date
-   !!---------------------------------------------------------------------
-   !!                     ***  oasis_date.h90  ***  
+   !!======================================================================
+   !!                     ***  MODULE oasis_date  ***  
    !!   Date and related information required to couple NEMO via OASIS4
    !!   Made separate from cpl_oasis4 module to allow wider use.
-   !!=====================================================================
-   !! History :   
-   !!   9.0  !  05-12  (R. Hill, Met. Office) Original code
+   !!======================================================================
+   !! History :  2.0  ! 2005-12  (R. Hill, Met. Office) Original code
    !!----------------------------------------------------------------------
 #if defined key_oasis4
+   !!----------------------------------------------------------------------
+   !!   'key_oasis4'                                   coupled with OASIS-4
+   !!----------------------------------------------------------------------
+   !##################### WARNING coupled mode ###############################
+   !##################### WARNING coupled mode ###############################
+   !   Following line must be enabled if coupling with OASIS
+   !   USE PRISM
+   !##################### WARNING coupled mode ###############################
+   !##################### WARNING coupled mode ###############################
+
+   PRIVATE
+   IMPLICIT NONE
+
+   INTEGER, PUBLIC :: date_err, date_info   !: 
+
+   TYPE(PRISM_Time_struct), PUBLIC ::   dates            !: date info for send operation
+   TYPE(PRISM_Time_struct), PUBLIC ::   dates_bound(2)   !: date info for send operation
+   TYPE(PRISM_Time_struct), PUBLIC ::   dater            !: date info for receive operation
+   TYPE(PRISM_Time_struct), PUBLIC ::   dater_bound(2)   !: date info for receive operation
+   TYPE(PRISM_Time_struct), PUBLIC ::   tmpdate          !:
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option            Dummy module                   NO OASIS-4
+   !!----------------------------------------------------------------------
+#endif
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-!##################### WARNING coupled mode ###############################
-!##################### WARNING coupled mode ###############################
-!   Following line must be enabled if coupling with OASIS
-!   USE PRISM
-!##################### WARNING coupled mode ###############################
-!##################### WARNING coupled mode ###############################
-
-   IMPLICIT NONE
-
-   INTEGER :: date_err, date_info
-   !!---------------------------------------------------------------------
-   TYPE(PRISM_Time_struct), PUBLIC    :: dates          ! date info for send operation
-   TYPE(PRISM_Time_struct), PUBLIC    :: dates_bound(2) ! date info for send operation
-   TYPE(PRISM_Time_struct), PUBLIC    :: dater          ! date info for receive operation
-   TYPE(PRISM_Time_struct), PUBLIC    :: dater_bound(2) ! date info for receive operation
-   TYPE(PRISM_Time_struct), PUBLIC    :: tmpdate
-#endif
+   !!======================================================================
 END MODULE oasis4_date
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 2620)
@@ -19,9 +19,11 @@
    USE par_ice_2        ! LIM-2 parameters
 # endif
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC sbc_ice_alloc ! called in nemogcm.F90
+   PUBLIC sbc_ice_alloc ! called in iceini(_2).F90
 
 # if defined  key_lim2
@@ -48,9 +50,9 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice   !: albedo of ice
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   utau_ice    !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau_ice    !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr1_i0      !: 1st Qsr fraction penetrating inside ice cover    [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fr2_i0      !: 2nd Qsr fraction penetrating inside ice cover    [-]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   emp_ice     !: sublimation-snow budget over ice             [kg/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice  !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice  !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0    !: 1st Qsr fraction penetrating inside ice cover    [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0    !: 2nd Qsr fraction penetrating inside ice cover    [-]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice   !: sublimation-snow budget over ice             [kg/m2]
 
 # if defined key_lim3
@@ -65,13 +67,8 @@
 CONTAINS
 
-   FUNCTION sbc_ice_alloc()
+   INTEGER FUNCTION sbc_ice_alloc()
       !!----------------------------------------------------------------------
-     !!                     ***  FUNCTION sbc_ice_alloc  ***
-     !!
-     !! ** Purpose :   Allocate all the dynamic arrays in the modules
+      !!                     ***  FUNCTION sbc_ice_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: sbc_ice_alloc   ! return value
-      !!----------------------------------------------------------------------
-      !
       ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     &
          &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     &
@@ -80,7 +77,9 @@
          &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     ,     &
          &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     &
-         &      emp_ice(jpi,jpj)                              , STAT=sbc_ice_alloc)
-      !
-  END FUNCTION sbc_ice_alloc
+         &      emp_ice(jpi,jpj)                              , STAT=sbc_ice_alloc )
+         !
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_ice_alloc )
+      IF( sbc_ice_alloc > 0 )   CALL ctl_warn('sbc_ice_alloc: allocation of arrays failed')
+   END FUNCTION sbc_ice_alloc
 
 #else
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2620)
@@ -10,15 +10,23 @@
    !!            3.3  ! 2010-10  (J. Chanut, C. Bricaud)  add the surface pressure forcing
    !!----------------------------------------------------------------------
-   USE par_oce          ! ocean parameters
+
+   !!----------------------------------------------------------------------
+   !!   sbc_oce_alloc : allocation of sbc arrays
+   !!   sbc_tau2wnd   : wind speed estimated from wind stress
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC sbc_oce_alloc ! routine called in nemogcm.F90
-
+   PUBLIC   sbc_oce_alloc   ! routine called in sbcmod.F90
+   PUBLIC   sbc_tau2wnd     ! routine called in several sbc modules
+   
    !!----------------------------------------------------------------------
    !!           Namelist for the Ocean Surface Boundary Condition
    !!----------------------------------------------------------------------
-   !                                             !! * namsbc namelist *
+   !                                            !!* namsbc namelist *
    LOGICAL , PUBLIC ::   ln_ana      = .FALSE.   !: analytical boundary condition flag
    LOGICAL , PUBLIC ::   ln_flx      = .FALSE.   !: flux      formulation
@@ -75,51 +83,43 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ssh_m     !: mean (nn_fsbc time-step) sea surface height                [m]
 
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!======================================================================
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   FUNCTION sbc_oce_alloc()
+   INTEGER FUNCTION sbc_oce_alloc()
       !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_oce_alloc  ***
+      !!                  ***  FUNCTION sbc_oce_alloc  ***
       !!---------------------------------------------------------------------
-      USE in_out_manager, ONLY: ctl_warn
-      IMPLICIT none
-      INTEGER :: sbc_oce_alloc
-      ! Local variables
       INTEGER :: ierr(4)
       !!---------------------------------------------------------------------
-
       ierr(:) = 0
-
-      ALLOCATE(utau(jpi,jpj),    utau_b(jpi,jpj),                      &
-               vtau(jpi,jpj),    vtau_b(jpi,jpj),                      &
-               taum(jpi,jpj),    wndm(jpi,jpj)  , Stat=ierr(1)) 
-
-      ALLOCATE(qsr(jpi,jpj),     qns(jpi,jpj),    qns_b(jpi,jpj),      &
-               qsr_tot(jpi,jpj), qns_tot(jpi,jpj),                     &
-               emp(jpi,jpj),     emp_b(jpi,jpj),                       &
-               emps(jpi,jpj),    emps_b(jpi,jpj), emp_tot(jpi,jpj),    &
-               Stat=ierr(2))
-
-      ALLOCATE(rnf(jpi,jpj),          rnf_b(jpi,jpj),                       &
-               sbc_tsc(jpi,jpj,jpts), sbc_tsc_b(jpi,jpj,jpts),         &  
-               qsr_hc(jpi,jpj,jpk) ,  qsr_hc_b(jpi,jpj,jpk),  Stat=ierr(3))
-
-      ALLOCATE(tprecip(jpi,jpj),      sprecip(jpi,jpj), fr_i(jpi,jpj), &
+      !
+      ALLOCATE( utau(jpi,jpj) , utau_b(jpi,jpj) , taum(jpi,jpj) ,     &
+         &      vtau(jpi,jpj) , vtau_b(jpi,jpj) , wndm(jpi,jpj) , STAT=ierr(1) ) 
+         !
+      ALLOCATE( qns_tot(jpi,jpj) , qns   (jpi,jpj) , qns_b(jpi,jpj),      &
+         &      qsr_tot(jpi,jpj) , qsr   (jpi,jpj) ,                     &
+         &      emp    (jpi,jpj) , emp_b (jpi,jpj) ,                       &
+         &      emps   (jpi,jpj) , emps_b(jpi,jpj) , emp_tot(jpi,jpj) , STAT=ierr(2))
+         !
+      ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
+         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
+         !
+      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
 #if defined key_cpl_carbon_cycle
-               atm_co2(jpi,jpj),                                       &
+         &      atm_co2(jpi,jpj) ,                                        &
 #endif
-               ssu_m(jpi,jpj), ssv_m(jpi,jpj), sst_m(jpi,jpj),         &
-               sss_m(jpi,jpj), ssh_m(jpi,jpj), Stat=ierr(4))
-
-      sbc_oce_alloc = MAXVAL(ierr)
-
-      IF(sbc_oce_alloc > 0)THEN
-         CALL ctl_warn('sbc_oce_alloc: allocation of arrays failed.')
-      END IF
-
+         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
+         &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+         !
+      sbc_oce_alloc = MAXVAL( ierr )
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_oce_alloc )
+      IF( sbc_oce_alloc > 0 )   CALL ctl_warn('sbc_oce_alloc: allocation of arrays failed')
+      !
    END FUNCTION sbc_oce_alloc
 
@@ -139,6 +139,4 @@
       REAL(wp) ::   ztx, zty, ztau, zcoef ! temporary variables
       INTEGER  ::   ji, jj                ! dummy indices
-      !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
       !!---------------------------------------------------------------------
       zcoef = 0.5 / ( zrhoa * zcdrag ) 
@@ -154,6 +152,7 @@
       END DO
       CALL lbc_lnk( wndm(:,:) , 'T', 1. )
-
+      !
    END SUBROUTINE sbc_tau2wnd
 
+   !!======================================================================
 END MODULE sbc_oce
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcana.F90	(revision 2620)
@@ -42,7 +42,6 @@
    !! NEMO/OPA 3.3 , 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
 
@@ -92,8 +91,8 @@
       ENDIF
 
-      qns   (:,:) = rn_qns0
-      qsr   (:,:) = rn_qsr0
-      emp   (:,:) = rn_emp0
-      emps  (:,:) = rn_emp0
+      qns (:,:) = rn_qns0
+      qsr (:,:) = rn_qsr0
+      emp (:,:) = rn_emp0
+      emps(:,:) = rn_emp0
    
       ! Increase the surface stress to its nominal value during the first nn_tau000 time-steps
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcapr.F90	(revision 2620)
@@ -19,6 +19,6 @@
    USE in_out_manager  ! I/O manager
    USE lib_fortran     ! distribued memory computing library
-   USE iom
-   USE restart
+   USE iom             ! IOM library
+   USE restart         ! ocean restart
 
    IMPLICIT NONE
@@ -32,21 +32,20 @@
    LOGICAL, PUBLIC ::   ln_ref_apr = .FALSE.  !: ref. pressure: global mean Patm (F) or a constant (F)
 
-   REAL(wp), ALLOCATABLE, PUBLIC, DIMENSION(:,:) ::   ssh_ib    ! Inverse barometer now    sea surface height     [m]
-   REAL(wp), ALLOCATABLE, PUBLIC, DIMENSION(:,:) ::   ssh_ibb   ! Inverse barometer before sea surface height     [m]
-   REAL(wp), ALLOCATABLE, PUBLIC, DIMENSION(:,:) ::   apr       ! atmospheric pressure at kt                   [N/m2]
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   ssh_ib    ! Inverse barometer now    sea surface height   [m]
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   ssh_ibb   ! Inverse barometer before sea surface height   [m]
+   REAL(wp), ALLOCATABLE, SAVE, PUBLIC, DIMENSION(:,:) ::   apr       ! atmospheric pressure at kt                 [N/m2]
    
    REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure          [N/m2]
    REAL(wp) ::   tarea                ! whole domain mean masked ocean surface
    REAL(wp) ::   r1_grau              ! = 1.e0 / (grav * rau0)
-
    
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:)   ::   sf_apr   ! structure of input fields (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE, SAVE, DIMENSION(:) ::   sf_apr   ! structure of input fields (file informations, fields read)
 
    !! * Substitutions
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010) 
-   !! $Id $
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -91,7 +90,6 @@
          !
          ALLOCATE( sf_apr(1), STAT=ierror )           !* allocate and fill sf_sst (forcing structure) with sn_sst
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'sbc_apr: unable to allocate sf_apr structure' )   ;   RETURN
-         ENDIF
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_apr: unable to allocate sf_apr structure' )
+         !
          CALL fld_fill( sf_apr, (/ sn_apr /), cn_dir, 'sbc_apr', 'Atmospheric pressure ', 'namsbc_apr' )
                                 ALLOCATE( sf_apr(1)%fnow(jpi,jpj,1)   )
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 2620)
@@ -43,5 +43,4 @@
    PUBLIC sbc_blk_clio        ! routine called by sbcmod.F90 
    PUBLIC blk_ice_clio        ! routine called by sbcice_lim.F90 
-   PUBLIC sbc_blk_clio_alloc  ! routine called by nemogcm.F90
 
    INTEGER , PARAMETER ::   jpfld   = 7           ! maximum number of files to read 
@@ -77,29 +76,14 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stauc        ! cloud optical depth 
    
-   REAL(wp)  ::   zeps    = 1.e-20                ! constant values
-   REAL(wp)  ::   zeps0   = 1.e-13  
+   REAL(wp) ::   eps20  = 1.e-20   ! constant values
    
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$ 
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
-
-   FUNCTION sbc_blk_clio_alloc()
-      !!---------------------------------------------------------------------
-      !!                 ***  ROUTINE sbc_blk_clio_alloc  ***
-      !!---------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: sbc_blk_clio_alloc
-      !!---------------------------------------------------------------------
-
-      ALLOCATE(sbudyko(jpi,jpj), &
-               stauc(jpi,jpj),   &
-               Stat=sbc_blk_clio_alloc)
-
-   END FUNCTION sbc_blk_clio_alloc
 
    SUBROUTINE sbc_blk_clio( kt )
@@ -134,5 +118,5 @@
       !!
       INTEGER  ::   ifpr, jfpr         ! dummy indices
-      INTEGER  ::   ierror             ! return error code
+      INTEGER  ::   ierr0, ierr1, ierr2, ierr3   ! return error code
       !!
       CHARACTER(len=100) ::  cn_dir                            !   Root directory for location of CLIO files
@@ -171,14 +155,17 @@
          
          ! set sf structure
-         ALLOCATE( sf(jpfld), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'sbc_blk_clio: unable to allocate sf structure' )   ;   RETURN
-         ENDIF
+         ALLOCATE( sf(jpfld), STAT=ierr0 )
+         IF( ierr0 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate sf structure' )
          DO ifpr= 1, jpfld
-            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
-         END DO
+            ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) , STAT=ierr1)
+            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) , STAT=ierr2 )
+         END DO
+         IF( ierr1+ierr2 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate sf array structure' )
          ! fill sf with slf_i and control print
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_clio', 'flux formulation for ocean surface boundary condition', 'namsbc_clio' )
+         
+         ! allocate sbcblk clio arrays
+         ALLOCATE( sbudyko(jpi,jpj) , stauc(jpi,jpj), STAT=ierr3 )
+         IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate arrays' )
          !
       ENDIF
@@ -330,5 +317,5 @@
             zdeltaq = zqatm - zqsato
             ztvmoy  = ztatm * ( 1. + 2.2e-3 * ztatm * zqatm )
-            zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) * ztvmoy, zeps )
+            zdenum  = MAX( sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) * ztvmoy, eps20 )
             zdtetar = zdteta / zdenum
             ztvmoyr = ztvmoy * ztvmoy * zdeltaq / zdenum
@@ -352,5 +339,5 @@
             zpsil   = zpsih
             
-            zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) / grav, zeps )
+            zvatmg         = MAX( 0.032 * 1.5e-3 * sf(jp_wndm)%fnow(ji,jj,1) * sf(jp_wndm)%fnow(ji,jj,1) / grav, eps20 )
             zcmn           = vkarmn / LOG ( 10. / zvatmg )
             zchn           = 0.0327 * zcmn
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2620)
@@ -78,5 +78,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
 
    SUBROUTINE sbc_blk_core( kt )
@@ -167,10 +166,8 @@
          !
          ALLOCATE( sf(jfld), STAT=ierror )         ! set sf structure
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'sbc_blk_core: unable to allocate sf structure' )   ;   RETURN
-         ENDIF
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_core: unable to allocate sf structure' )
          DO ifpr= 1, jfld
             ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
-            IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
+            IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
          END DO
          !                                         ! fill sf with slf_i and control print
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 2620)
@@ -58,5 +58,4 @@
    PUBLIC   sbc_cpl_ice_tau    ! routine called by sbc_ice_lim(_2).F90
    PUBLIC   sbc_cpl_ice_flx    ! routine called by sbc_ice_lim(_2).F90
-   PUBLIC   sbc_cpl_init_alloc ! routine called by nemogcm.F90
 
    INTEGER, PARAMETER ::   jpr_otx1   =  1            ! 3 atmosphere-ocean stress components on grid 1
@@ -155,5 +154,5 @@
    INTEGER , ALLOCATABLE, SAVE, DIMENSION(    :) ::   nrcvinfo           ! OASIS info argument
 
-#if ! defined key_lim2 && ! defined key_lim3
+#if ! defined key_lim2   &&   ! defined key_lim3
    ! quick patch to be able to run the coupled model without sea-ice...
    INTEGER, PARAMETER               ::   jpl = 1 
@@ -173,34 +172,25 @@
 CONTAINS
   
-   FUNCTION sbc_cpl_init_alloc()
-      !!----------------------------------------------------------------------
-      !!             ***  ROUTINE sbc_cpl_init_alloc  ***
-      !!----------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: sbc_cpl_init_alloc
+   INTEGER FUNCTION sbc_cpl_alloc()
+      !!----------------------------------------------------------------------
+      !!             ***  FUNCTION sbc_cpl_alloc  ***
+      !!----------------------------------------------------------------------
       INTEGER :: ierr(2)
       !!----------------------------------------------------------------------
-
       ierr(:) = 0
-
-      ALLOCATE(albedo_oce_mix(jpi,jpj), &
-               frcv(jpi,jpj,jprcv),     &
-               nrcvinfo(jprcv),  Stat=Stat=ierr(1))
-
+      !
+      ALLOCATE( albedo_oce_mix(jpi,jpj), frcv(jpi,jpj,jprcv), nrcvinfo(jprcv),  STAT=ierr(1) )
+      !
 #if ! defined key_lim2 && ! defined key_lim3
       ! quick patch to be able to run the coupled model without sea-ice...
-      ALLOCATE(hicif(jpi,jpj), hsnif(jpi,jpj), u_ice(jpi,jpj),  &
-               v_ice(jpi,jpj), fr1_i0(jpi,jpj),fr2_i0(jpi,jpj), &
-               tn_ice(jpi,jpj,jpl), alb_ice(jpi,jpj,jpl),       &
-               Stat=ierr(2) )
+      ALLOCATE( hicif(jpi,jpj) , u_ice(jpi,jpj) , fr1_i0(jpi,jpj) , tn_ice (jpi,jpj,jpl) ,     &
+                hsnif(jpi,jpj) , v_ice(jpi,jpj) , fr2_i0(jpi,jpj) , alb_ice(jpi,jpj,jpl) , STAT=ierr(2) )
 #endif
-
-      sbc_cpl_init_alloc = MAXVAL(ierr)
-
-      IF(sbc_cpl_init_alloc > 0)THEN
-         CALL ctl_warn('sbc_cpl_init_alloc: allocation of arrays failed.')
-      END IF
-
-    END FUNCTION sbc_cpl_init_alloc
+      sbc_cpl_init_alloc = MAXVAL( ierr )
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_cpl_alloc )
+      IF( sbc_cpl_alloc > 0 )   CALL ctl_warn('sbc_cpl_alloc: allocation of arrays failed')
+      !
+   END FUNCTION sbc_cpl_alloc
+
 
    SUBROUTINE sbc_cpl_init( k_ice )     
@@ -233,7 +223,6 @@
       !!---------------------------------------------------------------------
 
-      IF(.not. wrk_use(2,1,2))THEN
-         CALL ctl_stop('sbc_cpl_init: requested workspace arrays unavailable.')
-         RETURN
+      IF(.NOT. wrk_use(2,1,2) ) THEN
+         CALL ctl_stop('sbc_cpl_init: requested workspace arrays unavailable.')   ;   RETURN
       END IF
 
@@ -273,5 +262,5 @@
 
 #if defined key_cpl_carbon_cycle 
-      REWIND( numnam )                    ! ... read namlist namsbc_cpl_co2
+      REWIND( numnam )                    ! read namlist namsbc_cpl_co2
       READ  ( numnam, namsbc_cpl_co2 )
       IF(lwp) THEN                        ! control print
@@ -291,4 +280,7 @@
       cn_rcv_tau(1) = TRIM( cn_rcv_tau_nature )   ;   cn_rcv_tau(2) = TRIM( cn_rcv_tau_refere )
       cn_rcv_tau(3) = TRIM( cn_rcv_tau_orient )   ;   cn_rcv_tau(4) = TRIM( cn_rcv_tau_grid   )
+
+      !                                   ! allocate zdfric arrays
+      IF( sbc_cpl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_cpl_alloc : unable to allocate arrays' )
      
       ! ================================ !
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcdcy.F90	(revision 2620)
@@ -21,10 +21,11 @@
    IMPLICIT NONE
    PRIVATE
-   INTEGER, PUBLIC              ::   nday_qsr                    ! day when parameters were computed
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! parameters used to compute the diurnal cycle
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !     -       -         -           -      -
+   
+   INTEGER, PUBLIC ::   nday_qsr   !: day when parameters were computed
+   
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   raa , rbb  , rcc  , rab     ! diurnal cycle parameters
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   rtmd, rdawn, rdusk, rscal   !    -      -       -
   
    PUBLIC   sbc_dcy        ! routine called by sbc
-   PUBLIC   sbc_dcy_alloc  ! routine called by nemogcm.F90
 
    !!----------------------------------------------------------------------
@@ -35,24 +36,19 @@
 CONTAINS
 
-      FUNCTION sbc_dcy_alloc()
+      INTEGER FUNCTION sbc_dcy_alloc()
          !!----------------------------------------------------------------------
-         !!                ***  ROUTINE sbc_dcy_alloc  ***
+         !!                ***  FUNCTION sbc_dcy_alloc  ***
          !!----------------------------------------------------------------------
-         IMPLICIT none
-         INTEGER :: sbc_dcy_alloc
-         !!----------------------------------------------------------------------
-
-         ALLOCATE(raa(jpi,jpj),  rbb(jpi,jpj),   rcc(jpi,jpj),   rab(jpi,jpj),   &
-                  rtmd(jpi,jpj), rdawn(jpi,jpj), rdusk(jpi,jpj), rscal(jpi,jpj), &
-                  Stat=sbc_dcy_alloc)
-
-         IF(sbc_dcy_alloc /= 0)THEN
-            CALL ctl_warn('sbc_dcy_alloc: failed to allocate arrays.')
-         END IF
-
+         !
+         ALLOCATE( raa (jpi,jpj) , rbb  (jpi,jpj) , rcc  (jpi,jpj) , rab  (jpi,jpj) ,     &
+            &      rtmd(jpi,jpj) , rdawn(jpi,jpj) , rdusk(jpi,jpj) , rscal(jpi,jpj) , STAT=sbc_dcy_alloc )
+            !
+         IF( lk_mpp             )   CALL mpp_sum ( sbc_dcy_alloc )
+         IF( sbc_dcy_alloc /= 0 )   CALL ctl_warn('sbc_dcy_alloc: failed to allocate arrays')
+         !
       END FUNCTION sbc_dcy_alloc
 
 
-      FUNCTION sbc_dcy( pqsrin ) RESULT( zqsrout )
+   FUNCTION sbc_dcy( pqsrin ) RESULT( zqsrout )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_dcy  ***
@@ -85,12 +81,11 @@
       ! Initialization
       ! --------------
-      ztwopi    = 2. * rpi
-      zinvtwopi = 1. / ztwopi
-      zconvrad  = ztwopi / 360.
+      ztwopi    = 2._wp * rpi
+      zinvtwopi = 1._wp / ztwopi
+      zconvrad  = ztwopi / 360._wp
 
       ! When are we during the day (from 0 to 1)
-      zlo = ( REAL(nsec_day, wp) - 0.5   * rdttra(1) ) / rday
-      zup = zlo + ( REAL(nn_fsbc, wp) * rdttra(1) ) / rday
-
+      zlo = ( REAL(nsec_day, wp) - 0.5_wp * rdttra(1) ) / rday
+      zup = zlo + ( REAL(nn_fsbc, wp)     * rdttra(1) ) / rday
       !                                          
       IF( nday_qsr == -1 ) THEN       ! first time step only               
@@ -101,4 +96,6 @@
             WRITE(numout,*)
          ENDIF
+         ! allocate sbcdcy arrays
+         IF( sbc_dcy_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_dcy_alloc : unable to allocate arrays' )
          ! Compute rcc needed to compute the time integral of the diurnal cycle
          rcc(:,:) = zconvrad * glamt(:,:) - rpi
@@ -156,6 +153,6 @@
              END DO  
          END DO  
-         rdawn(:,:) = MOD((rdawn(:,:) + 1.), 1.)
-         rdusk(:,:) = MOD((rdusk(:,:) + 1.), 1.)
+         rdawn(:,:) = MOD( (rdawn(:,:) + 1._wp), 1._wp )
+         rdusk(:,:) = MOD( (rdusk(:,:) + 1._wp), 1._wp )
 
          !     2.2 Compute the scalling function:
@@ -185,5 +182,5 @@
          ztmp = rday / ( rdttra(1) * REAL(nn_fsbc, wp) )
          rscal(:,:) = rscal(:,:) * ztmp
-
+         !
       ENDIF 
 
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcflx.F90	(revision 2620)
@@ -42,5 +42,5 @@
    !! NEMO/OPA 3.3 , 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
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcfwb.F90	(revision 2620)
@@ -28,13 +28,10 @@
    PRIVATE
 
-   PUBLIC   sbc_fwb       ! routine called by step
-   PUBLIC   sbc_fwb_alloc ! routine called in nemogcm.F90
-
-   REAL(wp) ::   a_fwb_b            ! annual domain averaged freshwater budget
-   REAL(wp) ::   a_fwb              ! for 2 year before (_b) and before year.
-   REAL(wp) ::   fwfold             ! fwfold to be suppressed
-   REAL(wp) ::   area               ! global mean ocean surface (interior domain)
-
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e1e2    ! area of the interior domain (e1t*e2t)
+   PUBLIC   sbc_fwb    ! routine called by step
+
+   REAL(wp) ::   a_fwb_b   ! annual domain averaged freshwater budget
+   REAL(wp) ::   a_fwb     ! for 2 year before (_b) and before year.
+   REAL(wp) ::   fwfold    ! fwfold to be suppressed
+   REAL(wp) ::   area      ! global mean ocean surface (interior domain)
 
    !! * Substitutions
@@ -47,21 +44,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   FUNCTION sbc_fwb_alloc()
-      !!---------------------------------------------------------------------
-      !!                ***  ROUTINE sbc_fwb_alloc  ***
-      !!---------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: sbc_fwb_alloc
-      !!---------------------------------------------------------------------
-
-     ALLOCATE(e1e2(jpi,jpj), Stat=sbc_fwb_alloc)
-
-     IF(sbc_fwb_alloc /= 0)THEN
-        CALL ctl_warn('sbc_fwb_alloc: failed to allocate array.')
-     END IF
-
-   END FUNCTION sbc_fwb_alloc
-
 
    SUBROUTINE sbc_fwb( kt, kn_fwb, kn_fsbc )
@@ -93,7 +73,6 @@
       !!----------------------------------------------------------------------
       !
-      IF( .NOT. wrk_use(2, 1,2,3,4,5))THEN
-         CALL ctl_stop('sbc_fwb: requested workspace arrays are unavailable.')
-         RETURN
+      IF( .NOT. wrk_use(2, 1,2,3,4,5) ) THEN
+         CALL ctl_stop('sbc_fwb: requested workspace arrays are unavailable')   ;   RETURN
       END IF
       !
@@ -110,6 +89,5 @@
          IF( kn_fwb == 3 .AND. nn_sssr /= 2 )   CALL ctl_stop( 'sbc_fwb: nn_fwb = 3 requires nn_sssr = 2, we stop ' )
          !
-         e1e2(:,:) = e1t(:,:) * e2t(:,:) 
-         area = glob_sum( e1e2(:,:) )           ! interior global domain surface
+         area = glob_sum( e1e2t(:,:) )           ! interior global domain surface
       ENDIF
       
@@ -120,5 +98,5 @@
          !
          IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
-            z_fwf = glob_sum( e1e2(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area   ! sum over the global domain
             emp (:,:) = emp (:,:) - z_fwf 
             emps(:,:) = emps(:,:) - z_fwf 
@@ -143,5 +121,5 @@
          IF( MOD( kt, ikty ) == 0 ) THEN
             a_fwb_b = a_fwb
-            a_fwb   = glob_sum( e1e2(:,:) * sshn(:,:) )   ! sum over the global domain
+            a_fwb   = glob_sum( e1e2t(:,:) * sshn(:,:) )   ! sum over the global domain
             a_fwb   = a_fwb * 1.e+3 / ( area * 86400. * 365. )     ! convert in Kg/m3/s = mm/s
 !!gm        !                                                      !!bug 365d year 
@@ -168,8 +146,8 @@
             ztmsk_neg(:,:) = tmask_i(:,:) - ztmsk_pos(:,:)
             !
-            zsurf_neg = glob_sum( e1e2(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp 
-            zsurf_pos = glob_sum( e1e2(:,:)*ztmsk_pos(:,:) )
+            zsurf_neg = glob_sum( e1e2t(:,:)*ztmsk_neg(:,:) )   ! Area filled by <0 and >0 erp 
+            zsurf_pos = glob_sum( e1e2t(:,:)*ztmsk_pos(:,:) )
             !                                                  ! fwf global mean 
-            z_fwf = glob_sum( e1e2(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
+            z_fwf = glob_sum( e1e2t(:,:) * ( emp(:,:) - rnf(:,:) ) ) / area
             !            
             IF( z_fwf < 0._wp ) THEN         ! spread out over >0 erp area to increase evaporation
@@ -181,9 +159,9 @@
             ENDIF
             !
-            zsum_fwf   = glob_sum( e1e2(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
+            zsum_fwf   = glob_sum( e1e2t(:,:) * z_fwf )         ! fwf global mean over <0 or >0 erp area
 !!gm :  zsum_fwf   = z_fwf * area   ???  it is right?  I think so....
             z_fwf_nsrf =  zsum_fwf / ( zsurf_tospread + rsmall )
             !                                                  ! weight to respect erp field 2D structure 
-            zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2(:,:) )
+            zsum_erp = glob_sum( ztmsk_tospread(:,:) * erp(:,:) * e1e2t(:,:) )
             z_wgt(:,:) = ztmsk_tospread(:,:) * erp(:,:) / ( zsum_erp + rsmall )
             !                                                  ! final correction term to apply
@@ -200,10 +178,10 @@
                IF( z_fwf < 0._wp ) THEN
                   WRITE(numout,*)'   z_fwf < 0'
-                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp+)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
                ELSE
                   WRITE(numout,*)'   z_fwf >= 0'
-                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2(:,:) )*1.e-9,' Sv'
+                  WRITE(numout,*)'   SUM(erp-)     = ', SUM( ztmsk_tospread(:,:)*erp(:,:)*e1e2t(:,:) )*1.e-9,' Sv'
                ENDIF
-               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2(:,:) )*1.e-9,' Sv'
+               WRITE(numout,*)'   SUM(empG)     = ', SUM( z_fwf*e1e2t(:,:) )*1.e-9,' Sv'
                WRITE(numout,*)'   z_fwf         = ', z_fwf      ,' Kg/m2/s'
                WRITE(numout,*)'   z_fwf_nsrf    = ', z_fwf_nsrf ,' Kg/m2/s'
@@ -218,7 +196,5 @@
       END SELECT
       !
-      IF( .NOT. wrk_release(2, 1,2,3,4,5))THEN
-         CALL ctl_stop('sbc_fwb: failed to release workspace arrays.')
-      END IF
+      IF( .NOT. wrk_release(2, 1,2,3,4,5) )   CALL ctl_stop('sbc_fwb: failed to release workspace arrays')
       !
    END SUBROUTINE sbc_fwb
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 2620)
@@ -32,7 +32,6 @@
    !! NEMO/OPA 3.3 , 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
 
@@ -78,10 +77,7 @@
 
          ALLOCATE( sf_ice(1), STAT=ierror )
-         IF( ierror > 0 ) THEN
-            CALL ctl_stop( 'sbc_ice_if: unable to allocate sf_ice structure' )   ;   RETURN
-         ENDIF
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ice_if: unable to allocate sf_ice structure' )
          ALLOCATE( sf_ice(1)%fnow(jpi,jpj,1) )
-         IF( sn_ice%ln_tint ) ALLOCATE( sf_ice(1)%fdta(jpi,jpj,1,2) )
-
+         IF( sn_ice%ln_tint )   ALLOCATE( sf_ice(1)%fdta(jpi,jpj,1,2) )
 
          ! fill sf_ice with sn_ice and control print
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 2620)
@@ -102,6 +102,4 @@
       IF( .NOT. wrk_use(3, 1,2) ) THEN
          CALL ctl_stop( 'sbc_ice_lim: requested workspace arrays are unavailable.' )   ;   RETURN
-      ELSEIF( jpl > jpk ) THEN
-         CALL ctl_stop( 'sbc_ice_lim: extent of 3rd dimension of workspace arrays needs to exceed jpk.' )   ;   RETURN
       ENDIF
 
@@ -254,7 +252,5 @@
 !!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
-      IF( .NOT. wrk_release(3, 1,2) ) THEN
-         CALL ctl_stop( 'sbc_ice_lim: failed to release workspace arrays.' )
-      END IF
+      IF( .NOT. wrk_release(3, 1,2) )   CALL ctl_stop( 'sbc_ice_lim: failed to release workspace arrays.' )
       !
    END SUBROUTINE sbc_ice_lim
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2620)
@@ -84,5 +84,5 @@
       !!---------------------------------------------------------------------
       USE wrk_nemo, ONLY: wrk_use, wrk_release
-      USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2, wrk_3d_3
+      USE wrk_nemo, ONLY: wrk_3d_1, wrk_3d_2, wrk_3d_3   ! 3D workspace
       !!
       INTEGER, INTENT(in) ::   kt      ! ocean time step
@@ -96,12 +96,11 @@
       !!----------------------------------------------------------------------
 
-      IF(.NOT. wrk_use(3, 1,2,3))THEN
-         CALL ctl_stop('sbc_ice_lim_2: requested workspace arrays are unavailable.')
-         RETURN
+      IF(.NOT. wrk_use(3, 1,2,3) ) THEN
+         CALL ctl_stop('sbc_ice_lim_2: requested workspace arrays are unavailable.')   ;   RETURN
       END IF
       ! Use pointers to access only sub-arrays of workspaces
       zalb_ice_os => wrk_3d_1(:,:,1:1)
       zalb_ice_cs => wrk_3d_2(:,:,1:1)
-            zsist => wrk_3d_3(:,:,1:1)
+      zsist       => wrk_3d_3(:,:,1:1)
 
       IF( kt == nit000 ) THEN
@@ -204,8 +203,8 @@
          !                                             ! Ice surface fluxes in coupled mode 
          IF( ksbc == 5 )   CALL sbc_cpl_ice_flx( reshape( frld, (/jpi,jpj,1/) ),                 &
-      &                                                   qns_tot, qns_ice, qsr_tot , qsr_ice,   &
-      &                                                   emp_tot, emp_ice, dqns_ice, sprecip,   &
-      !                                      optional arguments, used only in 'mixed oce-ice' case
-      &                                                   palbi = zalb_ice_cs, psst = sst_m, pist = sist )
+            &                                             qns_tot, qns_ice, qsr_tot , qsr_ice,   &
+            &                                             emp_tot, emp_ice, dqns_ice, sprecip,   &
+            !                                optional arguments, used only in 'mixed oce-ice' case
+            &                                             palbi = zalb_ice_cs, psst = sst_m, pist = sist )
 #endif
                            CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
@@ -229,7 +228,5 @@
       IF( ln_limdyn    )   CALL lim_sbc_tau_2( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
       !
-      IF(.NOT. wrk_release(3, 1,2,3))THEN
-         CALL ctl_stop('sbc_ice_lim_2: failed to release workspace arrays.')
-      END IF
+      IF(.NOT. wrk_release(3, 1,2,3) )   CALL ctl_stop('sbc_ice_lim_2: failed to release workspace arrays')
       !
    END SUBROUTINE sbc_ice_lim_2
@@ -241,6 +238,5 @@
 CONTAINS
    SUBROUTINE sbc_ice_lim_2 ( kt, ksbc )     ! Dummy routine
-      INTEGER, INTENT(in) ::   kt      
-      INTEGER, INTENT(in) ::   ksbc    
+      INTEGER, INTENT(in) ::   kt, ksbc    
       WRITE(*,*) 'sbc_ice_lim_2: You should not have seen this print! error?', kt, ksbc
    END SUBROUTINE sbc_ice_lim_2
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2620)
@@ -57,7 +57,7 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO-consortium (2010) 
+   !! NEMO/OPA 4.0 , NEMO-consortium (2011) 
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -76,6 +76,6 @@
       INTEGER ::   icpt      ! temporary integer
       !!
-      NAMELIST/namsbc/ nn_fsbc, ln_ana  , ln_flx, ln_blk_clio, ln_blk_core, ln_cpl    ,   &
-         &             ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana , ln_flx  , ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_apr_dyn, nn_ice , ln_dm2dc, ln_rnf     , ln_ssr     , nn_fwb
       !!----------------------------------------------------------------------
 
@@ -117,4 +117,7 @@
          WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
       ENDIF
+
+      !                              ! allocate sbc arrays
+      IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
 
       !                          ! Checks:
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 2620)
@@ -30,7 +30,6 @@
    PUBLIC   sbc_rnf       ! routine call in sbcmod module
    PUBLIC   sbc_rnf_div   ! routine called in sshwzv module
-   PUBLIC   sbc_rnf_alloc ! routine called in nemogcm module
-
-   !                                                      !!* namsbc_rnf namelist *
+
+   !                                                     !!* namsbc_rnf namelist *
    CHARACTER(len=100), PUBLIC ::   cn_dir       = './'    !: Root directory for location of ssr files
    LOGICAL           , PUBLIC ::   ln_rnf_depth = .false. !: depth       river runoffs attribute specified in a file
@@ -48,10 +47,10 @@
    REAL(wp)          , PUBLIC ::   rn_rfact     = 1._wp   !: multiplicative factor for runoff
 
-   INTEGER , PUBLIC                          ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk            !: river mouth mask (hori.)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z          !: river mouth mask (vert.)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   h_rnf             !: depth of runoff in m
-   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nk_rnf            !: depth of runoff in model levels
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rnf_tsc_b, rnf_tsc  !: before and now T & S contents of runoffs  [K.m/s & PSU.m/s]
+   INTEGER , PUBLIC  ::   nkrnf = 0         !: nb of levels over which Kz is increased at river mouths
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rnfmsk              !: river mouth mask (hori.)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)     ::   rnfmsk_z            !: river mouth mask (vert.)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   h_rnf               !: depth of runoff in m
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   nk_rnf              !: depth of runoff in model levels
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   rnf_tsc_b, rnf_tsc  !: before and now T & S runoff contents   [K.m/s & PSU.m/s]
    
    REAL(wp) ::   r1_rau0   ! = 1 / rau0 
@@ -70,21 +69,14 @@
 CONTAINS
 
-   FUNCTION sbc_rnf_alloc()
+   INTEGER FUNCTION sbc_rnf_alloc()
       !!----------------------------------------------------------------------
       !!                ***  ROUTINE sbc_rnf_alloc  ***
       !!----------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: sbc_rnf_alloc
-      !!----------------------------------------------------------------------
-
-      ALLOCATE(rnfmsk(jpi,jpj),         rnfmsk_z(jpk),         &
-               h_rnf(jpi,jpj),          nk_rnf(jpi,jpj),       &
-               rnf_tsc_b(jpi,jpj,jpts), rnf_tsc(jpi,jpj,jpts), &
-               Stat=sbc_rnf_alloc)
-
-      IF(sbc_rnf_alloc > 0)THEN
-         CALL ctl_warn('sbc_rnf_alloc: allocation of arrays failed.')
-      END IF
-
+      ALLOCATE( rnfmsk(jpi,jpj)         , rnfmsk_z(jpk)          ,     &
+         &      h_rnf (jpi,jpj)         , nk_rnf  (jpi,jpj)      ,     &
+         &      rnf_tsc_b(jpi,jpj,jpts) , rnf_tsc (jpi,jpj,jpts) , STAT=sbc_rnf_alloc)
+         !
+      IF( lk_mpp            )   CALL mpp_sum ( sbc_rnf_alloc )
+      IF( sbc_rnf_alloc > 0 )   CALL ctl_warn('sbc_rnf_alloc: allocation of arrays failed.')
    END FUNCTION sbc_rnf_alloc
 
@@ -296,4 +288,6 @@
       !                                   !   Type of runoff
       !                                   ! ==================
+      !                                         !==  allocate runoff arrays
+      IF( sbc_rnf_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_rnf_alloc : unable to allocate arrays' )
       !
       IF( ln_rnf_emp ) THEN                     !==  runoffs directly provided in the precipitations  ==!
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssm.F90	(revision 2620)
@@ -4,5 +4,5 @@
    !! Surface module :  provide time-mean ocean surface variables
    !!======================================================================
-   !! History :  9.0   !  06-07  (G. Madec)  Original code
+   !! History :  9.0  ! 2006-07  (G. Madec)  Original code
    !!            3.3  ! 2010-10  (C. Bricaud, G. Madec)  add the Patm forcing for sea-ice
    !!----------------------------------------------------------------------
@@ -32,7 +32,6 @@
    !! NEMO/OPA 3.3 , 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
 
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 2620)
@@ -25,9 +25,8 @@
    PRIVATE
 
-   PUBLIC   sbc_ssr       ! routine called in sbcmod
-   PUBLIC   sbc_ssr_alloc ! routine called in nemgcm
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp      !: evaporation damping   [kg/m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp      !: heat flux damping        [w/m2]
+   PUBLIC   sbc_ssr    ! routine called in sbcmod
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   erp   !: evaporation damping   [kg/m2/s]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qrp   !: heat flux damping        [w/m2]
 
    !                                           !!* Namelist namsbc_ssr *
@@ -46,26 +45,9 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   FUNCTION sbc_ssr_alloc()
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE sbc_ssr_alloc  ***
-      !!---------------------------------------------------------------------
-      IMPLICIT none
-      INTEGER :: sbc_ssr_alloc
-      !!---------------------------------------------------------------------
-
-      ALLOCATE(erp(jpi,jpj), qrp(jpi,jpj), Stat=sbc_ssr_alloc)
-
-      IF(sbc_ssr_alloc > 0)THEN
-         CALL ctl_warn('sbc_ssr_alloc: allocation of arrays failed.')
-      END IF
-
-   END FUNCTION sbc_ssr_alloc
 
    SUBROUTINE sbc_ssr( kt )
@@ -125,28 +107,34 @@
          ENDIF
 
-         IF( nn_sstr == 1 ) THEN      !* set sf_sst structure
+         IF( nn_sstr == 1 ) THEN      !* set sf_sst structure & allocate arrays
             !
             ALLOCATE( sf_sst(1), STAT=ierror )
-            IF( ierror > 0 ) THEN
-               CALL ctl_stop( 'sbc_ssr: unable to allocate sf_sst structure' )   ;   RETURN
-            ENDIF
-            ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1) )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst structure' )
+            ALLOCATE( sf_sst(1)%fnow(jpi,jpj,1), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst now array' )
             !
             ! fill sf_sst with sn_sst and control print
             CALL fld_fill( sf_sst, (/ sn_sst /), cn_dir, 'sbc_ssr', 'SST restoring term toward SST data', 'namsbc_ssr' )
-            IF( sf_sst(1)%ln_tint ) ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2) )
-         ENDIF
-         !
-         IF( nn_sssr >= 1 ) THEN      ! set sf_sss structure
+            IF( sf_sst(1)%ln_tint )   ALLOCATE( sf_sst(1)%fdta(jpi,jpj,1,2), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sst data array' )
+            !
+            ALLOCATE( qrp(jpi,jpj), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate qrp array' )
+         ENDIF
+         !
+         IF( nn_sssr >= 1 ) THEN      ! set sf_sss structure & allocate arrays
             !
             ALLOCATE( sf_sss(1), STAT=ierror )
-            IF( ierror > 0 ) THEN
-               CALL ctl_stop( 'sbc_ssr: unable to allocate sf_sss structure' )   ;   RETURN
-            ENDIF
-            ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1) )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss structure' )
+            ALLOCATE( sf_sss(1)%fnow(jpi,jpj,1), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss now array' )
             !
             ! fill sf_sss with sn_sss and control print
             CALL fld_fill( sf_sss, (/ sn_sss /), cn_dir, 'sbc_ssr', 'SSS restoring term toward SSS data', 'namsbc_ssr' )
-            IF( sf_sss(1)%ln_tint )ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2) )
+            IF( sf_sss(1)%ln_tint )   ALLOCATE( sf_sss(1)%fdta(jpi,jpj,1,2), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate sf_sss data array' )
+            !
+            ALLOCATE( erp(jpi,jpj), STAT=ierror )
+            IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_ssr: unable to allocate erp array' )
          ENDIF
          !
Index: /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2619)
+++ /branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 2620)
@@ -472,5 +472,4 @@
      
      USE dynzdf_exp,   ONLY: dyn_zdf_exp_alloc
-     USE geo2ocean,    ONLY: geo2oce_alloc
 #if   defined key_mpp_mpi  
      USE lib_mpp,      ONLY: lib_mpp_alloc
@@ -480,13 +479,5 @@
      USE obc_oce,      ONLY: obc_oce_alloc
 #endif
-     USE sbcblk_clio,  ONLY: sbc_blk_clio_alloc
-#if defined key_oasis3 || defined key_oasis4
-     USE sbccpl,       ONLY: sbc_cpl_init_alloc
-#endif
-     USE sbcdcy,       ONLY: sbc_dcy_alloc
-     USE sbcfwb,       ONLY: sbc_fwb_alloc
-     USE sbc_oce,      ONLY: sbc_oce_alloc
-     USE sbcrnf,       ONLY: sbc_rnf_alloc
-     USE sbcssr,       ONLY: sbc_ssr_alloc
+
      USE sol_oce,      ONLY: sol_oce_alloc
      USE solmat,       ONLY: sol_mat_alloc
@@ -565,5 +556,4 @@
 
       ierr = ierr + dyn_zdf_exp_alloc()
-      ierr = ierr + geo2oce_alloc()
 #if defined key_mpp_mpi 
       ierr = ierr + lib_mpp_alloc()
@@ -573,13 +563,5 @@
       ierr = ierr + obc_oce_alloc()
 #endif
-      ierr = ierr + sbc_blk_clio_alloc()
-#if defined key_oasis3 || defined key_oasis4
-      ierr = ierr + sbc_cpl_init_alloc()
-#endif
-      ierr = ierr + sbc_dcy_alloc()
-      ierr = ierr + sbc_fwb_alloc()
-      ierr = ierr + sbc_oce_alloc()
-      ierr = ierr + sbc_rnf_alloc()
-      ierr = ierr + sbc_ssr_alloc()
+
       ierr = ierr + sol_oce_alloc()
       ierr = ierr + sol_mat_alloc()
@@ -763,6 +745,6 @@
       RETURN
       !
-  END SUBROUTINE factorise
-
-  !!======================================================================
+   END SUBROUTINE factorise
+
+   !!======================================================================
 END MODULE nemogcm
