Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdy_par.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdy_par.F90	(revision 2185)
@@ -8,5 +8,5 @@
    !!            3.3  !  2010-09  (D. Storkey and E. O'Dea) update for Shelf configurations
    !!----------------------------------------------------------------------
-#if defined key_bdy
+#if defined   key_bdy
    !!----------------------------------------------------------------------
    !!   'key_bdy' :                    Unstructured Open Boundary Condition
@@ -16,10 +16,9 @@
    PUBLIC
 
-
-   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .TRUE.  !: Unstructured Ocean Boundary Condition flag
+   LOGICAL, PUBLIC, PARAMETER ::   lk_bdy  = .TRUE.   !: Unstructured Ocean Boundary Condition flag
    INTEGER, PUBLIC, PARAMETER ::   jpbdta  = 20000    !: Max length of bdy field in file
    INTEGER, PUBLIC, PARAMETER ::   jpbdim  = 20000    !: Max length of bdy field on a processor
-   INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000    !: Max number of time dumps per file
-   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 6	     !: Number of horizontal grid types used  (T, u, v, f)
+   INTEGER, PUBLIC, PARAMETER ::   jpbtime = 1000     !: Max number of time dumps per file
+   INTEGER, PUBLIC, PARAMETER ::   jpbgrd  = 6	      !: Number of horizontal grid types used  (T, u, v, f)
 #else
    !!----------------------------------------------------------------------
@@ -30,5 +29,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 2185)
@@ -16,6 +16,5 @@
    !!----------------------------------------------------------------------
    !!   bdy_dta    : read u, v, t, s data along open boundaries
-   !!   bdy_dta_bt : read depth-mean velocities and elevation along open
-   !!                boundaries        
+   !!   bdy_dta_bt : read depth-mean velocities and elevation along open boundaries        
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers
@@ -37,32 +36,31 @@
    PUBLIC   bdy_dta_bt 
 
-   INTEGER ::   numbdyt, numbdyu, numbdyv                      !: logical units for T-, U-, & V-points data file, resp.
-   INTEGER ::   ntimes_bdy                                     !: exact number of time dumps in data files
-   INTEGER ::   nbdy_b, nbdy_a                                 !: record of bdy data file for before and after model time step
-   INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             !: logical unit for T-, U- & V-points data file, resp.
-   INTEGER ::   ntimes_bdy_bt                                  !: exact number of time dumps in data files
-   INTEGER ::   nbdy_b_bt, nbdy_a_bt                           !: record of bdy data file for before and after model time step
-
-   INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           !: time array in seconds in each data file
-
-   REAL(wp) ::  zoffset                                        !: time offset between time origin in file & start time of model run
-
-   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   tbdydta, sbdydta     !: time interpolated values of T and S bdy data   
-   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   ubdydta, vbdydta     !: time interpolated values of U and V bdy data 
-   REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta !: Arrays used for time interpolation of bdy data   
-   REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            !: bdy data of ssh
+   INTEGER ::   numbdyt, numbdyu, numbdyv                      ! logical units for T-, U-, & V-points data file, resp.
+   INTEGER ::   ntimes_bdy                                     ! exact number of time dumps in data files
+   INTEGER ::   nbdy_b, nbdy_a                                 ! record of bdy data file for before and after time step
+   INTEGER ::   numbdyt_bt, numbdyu_bt, numbdyv_bt             ! logical unit for T-, U- & V-points data file, resp.
+   INTEGER ::   ntimes_bdy_bt                                  ! exact number of time dumps in data files
+   INTEGER ::   nbdy_b_bt, nbdy_a_bt                           ! record of bdy data file for before and after time step
+
+   INTEGER, DIMENSION (jpbtime) ::   istep, istep_bt           ! time array in seconds in each data file
+
+   REAL(wp) ::  zoffset                                        ! time offset between time origin in file & start time of model run
+
+   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   tbdydta, sbdydta     ! time interpolated values of T and S bdy data   
+   REAL(wp), DIMENSION(jpbdim,jpk,2) ::   ubdydta, vbdydta     ! time interpolated values of U and V bdy data 
+   REAL(wp), DIMENSION(jpbdim,2)     ::   ubtbdydta, vbtbdydta ! Arrays used for time interpolation of bdy data   
+   REAL(wp), DIMENSION(jpbdim,2)     ::   sshbdydta            ! bdy data of ssh
 
 #if defined key_lim2
-   REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          !: }
-   REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         !: } Arrays used for time interpolation of bdy data for ice variables
-   REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         !: }
+   REAL(wp), DIMENSION(jpbdim,2)     ::   frld_bdydta          ! }
+   REAL(wp), DIMENSION(jpbdim,2)     ::   hicif_bdydta         ! } Arrays used for time interp. of ice bdy data 
+   REAL(wp), DIMENSION(jpbdim,2)     ::   hsnif_bdydta         ! }
 #endif
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -78,5 +76,5 @@
       !!                the file. If so read it in. Time interpolate.
       !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt                             ! ocean time-step index (for timesplitting option, otherwise zero)
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index (for timesplitting option, otherwise zero)
       !!
       CHARACTER(LEN=80), DIMENSION(3) ::   clfile               ! names of input files
@@ -90,5 +88,5 @@
       INTEGER ::   itimer, totime
       INTEGER ::   ii, ij                                       ! array addresses
-      INTEGER ::   ipi, ipj, ipk, inum                          ! temporary integers (NetCDF read)
+      INTEGER ::   ipi, ipj, ipk, inum                          ! local integers (NetCDF read)
       INTEGER ::   iyear0, imonth0, iday0
       INTEGER ::   ihours0, iminutes0, isec0
@@ -102,6 +100,7 @@
       !!---------------------------------------------------------------------------
 
-      IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs ) THEN  ! If these are both false then this routine 
-                                                     ! does nothing.
+
+      IF( ln_bdy_dyn_frs .OR. ln_bdy_tra_frs    &
+         &               .OR. ln_bdy_ice_frs ) THEN  ! If these are both false then this routine does nothing
 
       ! -------------------- !
@@ -113,5 +112,6 @@
       ! Some time variables for monthly climatological forcing:
       ! *******************************************************
- !!gm  here  use directely daymod variables
+
+!!gm  here  use directely daymod calendar variables
  
       iman = INT( raamo )      ! Number of months in a year
@@ -132,6 +132,6 @@
          !                                             !-------------------!
          istep(:) = 0
-         nbdy_b    = 0
-         nbdy_a    = 0
+         nbdy_b   = 0
+         nbdy_a   = 0
 
          ! Get time information from bdy data file
@@ -162,9 +162,7 @@
             igrd_start = 1
             igrd_end   = 3
-            IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN
-               ! No T-grid file.
+            IF(.NOT. ln_bdy_tra_frs .AND. .NOT. ln_bdy_ice_frs) THEN       ! No T-grid file.
                igrd_start = 2
-            ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN
-               ! No U-grid or V-grid file.
+            ELSEIF ( .NOT. ln_bdy_dyn_frs ) THEN                           ! No U-grid or V-grid file.
                igrd_end   = 1         
             ENDIF
@@ -176,10 +174,7 @@
 
                SELECT CASE( igrd )
-                  CASE (1) 
-                     numbdyt = inum
-                  CASE (2) 
-                     numbdyu = inum
-                  CASE (3) 
-                     numbdyv = inum
+                  CASE (1)   ;   numbdyt = inum
+                  CASE (2)   ;   numbdyu = inum
+                  CASE (3)   ;   numbdyv = inum
                END SELECT
 
@@ -216,11 +211,10 @@
 
                ! Check that time array increases:
-        
                it = 1
-               DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 )
-                 it = it + 1
+               DO WHILE( zstepr(it+1) > zstepr(it) .AND. it /= ntimes_bdy - 1 ) 
+                  it = it + 1
                END DO
-
-               IF( it.NE.ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN
+               !
+               IF( it /= ntimes_bdy-1 .AND. ntimes_bdy > 1 ) THEN
                      WRITE(ctmp1,*) 'Check file: ', clfile(igrd)
                      CALL ctl_stop( 'Time array in unstructured boundary data files',   &
@@ -274,8 +268,6 @@
             ENDIF
 
-            IF( igrd_start == 1 ) THEN
-               istep(:) = istept(:)
-            ELSE
-               istep(:) = istepu(:)
+            IF( igrd_start == 1 ) THEN   ;   istep(:) = istept(:)
+            ELSE                         ;   istep(:) = istepu(:)
             ENDIF
 
@@ -302,5 +294,5 @@
             it = 1
             DO WHILE( istep(it+1) <= 0 .AND. it <= ntimes_bdy - 1 )
-              it = it + 1
+               it = it + 1
             END DO
             nbdy_b = it
@@ -315,43 +307,43 @@
          ! *****************************************************************
 
-         IF( nbdy_dta == 0) THEN      ! boundary data arrays are filled with initial conditions
+         IF( nbdy_dta == 0 ) THEN      ! boundary data arrays are filled with initial conditions
             !
             IF (ln_bdy_tra_frs) THEN
-              igrd = 1            ! T-points data 
-              DO ib = 1, nblen(igrd)
-                ii = nbi(ib,igrd)
-                ij = nbj(ib,igrd)
-                DO ik = 1, jpkm1
-                  tbdy(ib,ik) = tn(ii, ij, ik)
-                  sbdy(ib,ik) = sn(ii, ij, ik)
-                ENDDO
-              END DO
+               igrd = 1            ! T-points data 
+               DO ib = 1, nblen(igrd)
+                  ii = nbi(ib,igrd)
+                  ij = nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     tbdy(ib,ik) = tn(ii,ij,ik)
+                     sbdy(ib,ik) = sn(ii,ij,ik)
+                  END DO
+               END DO
             ENDIF
 
             IF(ln_bdy_dyn_frs) THEN
-              igrd = 2            ! U-points data 
-              DO ib = 1, nblen(igrd)
-                ii = nbi(ib,igrd)
-                ij = nbj(ib,igrd)
-                DO ik = 1, jpkm1
-                  ubdy(ib,ik) = un(ii, ij, ik)
-                ENDDO
-              END DO
-
-              igrd = 3            ! V-points data 
-              DO ib = 1, nblen(igrd)            
-                ii = nbi(ib,igrd)
-                ij = nbj(ib,igrd)
-                DO ik = 1, jpkm1
-                  vbdy(ib,ik) = vn(ii, ij, ik)
-                ENDDO
-              END DO
+               igrd = 2            ! U-points data 
+               DO ib = 1, nblen(igrd)
+                  ii = nbi(ib,igrd)
+                  ij = nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     ubdy(ib,ik) = un(ii, ij, ik)
+                  END DO
+               END DO
+               !
+               igrd = 3            ! V-points data 
+               DO ib = 1, nblen(igrd)            
+                  ii = nbi(ib,igrd)
+                  ij = nbj(ib,igrd)
+                  DO ik = 1, jpkm1
+                     vbdy(ib,ik) = vn(ii, ij, ik)
+                  END DO
+               END DO
             ENDIF
             !
 #if defined key_lim2
-            IF (ln_bdy_ice_frs) THEN
+            IF( ln_bdy_ice_frs ) THEN
                igrd = 1            ! T-points data
                DO ib = 1, nblen(igrd)
-                  frld_bdy(ib)  =  frld(nbi(ib,igrd), nbj(ib,igrd))
+                  frld_bdy (ib) =  frld(nbi(ib,igrd), nbj(ib,igrd))
                   hicif_bdy(ib) = hicif(nbi(ib,igrd), nbj(ib,igrd))
                   hsnif_bdy(ib) = hsnif(nbi(ib,igrd), nbj(ib,igrd))
@@ -377,4 +369,5 @@
 
             IF(ln_bdy_tra_frs) THEN
+               !
                igrd = 1                                           ! Temperature
                IF( nblendta(igrd) <=  0 ) THEN 
@@ -385,5 +378,5 @@
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyt, jpdom_unknown, 'votemper', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-
+               !
                DO ib = 1, nblen(igrd)
                   DO ik = 1, jpkm1
@@ -400,5 +393,5 @@
                ipi = nblendta(igrd)
                CALL iom_get ( numbdyt, jpdom_unknown, 'vosaline', zdta(1:ipi,1:ipj,1:ipk), nbdy_a )
-
+               !
                DO ib = 1, nblen(igrd)
                   DO ik = 1, jpkm1
@@ -408,6 +401,6 @@
             ENDIF  ! ln_bdy_tra_frs
  
-            IF(ln_bdy_dyn_frs) THEN
-
+            IF( ln_bdy_dyn_frs ) THEN
+               !
                igrd = 2                                           ! u-velocity
                IF ( nblendta(igrd) .le. 0 ) THEN 
@@ -440,5 +433,5 @@
 
 #if defined key_lim2
-            IF(ln_bdy_ice_frs) THEN
+            IF( ln_bdy_ice_frs ) THEN
               !
               igrd=1                                              ! leads fraction
@@ -468,109 +461,103 @@
 #endif
 
-            IF ((.NOT.ln_bdy_clim) .AND. (istep(1) > 0)) THEN
-               ! First data time is after start of run
-               ! Put first value in both time levels
+            IF( .NOT.ln_bdy_clim .AND. istep(1) > 0 ) THEN     ! First data time is after start of run
+               nbdy_b = nbdy_a                                 ! Put first value in both time levels
+               IF( ln_bdy_tra_frs ) THEN
+                 tbdydta(:,:,1) = tbdydta(:,:,2)
+                 sbdydta(:,:,1) = sbdydta(:,:,2)
+               ENDIF
+               IF( ln_bdy_dyn_frs ) THEN
+                 ubdydta(:,:,1) = ubdydta(:,:,2)
+                 vbdydta(:,:,1) = vbdydta(:,:,2)
+               ENDIF
+#if defined key_lim2
+               IF( ln_bdy_ice_frs ) THEN
+                  frld_bdydta (:,1) =  frld_bdydta(:,2)
+                  hicif_bdydta(:,1) = hicif_bdydta(:,2)
+                  hsnif_bdydta(:,1) = hsnif_bdydta(:,2)
+               ENDIF
+#endif
+            END IF
+            !
+         END IF   ! nbdy_dta == 0/1
+ 
+         ! In the case of constant boundary forcing fill bdy arrays once for all
+         IF( ln_bdy_clim .AND. ntimes_bdy == 1 ) THEN
+            IF( ln_bdy_tra_frs ) THEN
+               tbdy  (:,:) = tbdydta  (:,:,2)
+               sbdy  (:,:) = sbdydta  (:,:,2)
+            ENDIF
+            IF( ln_bdy_dyn_frs) THEN
+               ubdy  (:,:) = ubdydta  (:,:,2)
+               vbdy  (:,:) = vbdydta  (:,:,2)
+            ENDIF
+#if defined key_lim2
+            IF( ln_bdy_ice_frs ) THEN
+               frld_bdy (:) = frld_bdydta (:,2)
+               hicif_bdy(:) = hicif_bdydta(:,2)
+               hsnif_bdy(:) = hsnif_bdydta(:,2)
+            ENDIF
+#endif
+
+            IF( ln_bdy_tra_frs .OR. ln_bdy_ice_frs) CALL iom_close( numbdyt )
+            IF( ln_bdy_dyn_frs                    ) CALL iom_close( numbdyu )
+            IF( ln_bdy_dyn_frs                    ) CALL iom_close( numbdyv )
+         END IF
+         !
+      ENDIF                                            ! End if nit000
+
+
+      !                                                !---------------------!
+      IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN    !  at each time step  !
+         !                                             !---------------------!
+         ! Read one more record if necessary
+         !**********************************
+
+         IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000
+            nbdy_b = imois
+            nbdy_a = imois + 1
+            nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman
+            nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman
+            lect=.true.
+         ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN
+
+            IF( nbdy_a < ntimes_bdy ) THEN
+               nbdy_b = nbdy_a
+               nbdy_a = nbdy_a + 1
+               lect  =.true.
+            ELSE
+               ! We have reached the end of the file
+               ! put the last data time into both time levels
                nbdy_b = nbdy_a
                IF(ln_bdy_tra_frs) THEN
-                 tbdydta(:,:,1) = tbdydta(:,:,2)
-                 sbdydta(:,:,1) = sbdydta(:,:,2)
+                  tbdydta(:,:,1) =  tbdydta(:,:,2)
+                  sbdydta(:,:,1) =  sbdydta(:,:,2)
                ENDIF
                IF(ln_bdy_dyn_frs) THEN
-                 ubdydta(:,:,1) = ubdydta(:,:,2)
-                 vbdydta(:,:,1) = vbdydta(:,:,2)
+                  ubdydta(:,:,1) =  ubdydta(:,:,2)
+                  vbdydta(:,:,1) =  vbdydta(:,:,2)
                ENDIF
 #if defined key_lim2
-             IF( ln_bdy_ice_frs ) THEN
-               frld_bdydta (:,1) =  frld_bdydta(:,2)
-               hicif_bdydta(:,1) = hicif_bdydta(:,2)
-               hsnif_bdydta(:,1) = hsnif_bdydta(:,2)
-             ENDIF
-#endif
-            END IF
-
-         END IF ! nbdy_dta == 0/1
- 
-         ! In the case of constant boundary forcing fill bdy arrays once for all
-         IF ((ln_bdy_clim).AND.(ntimes_bdy==1)) THEN
-            IF(ln_bdy_tra_frs) THEN
-              tbdy  (:,:) = tbdydta  (:,:,2)
-              sbdy  (:,:) = sbdydta  (:,:,2)
-            ENDIF
-            IF(ln_bdy_dyn_frs) THEN
-              ubdy  (:,:) = ubdydta  (:,:,2)
-              vbdy  (:,:) = vbdydta  (:,:,2)
-            ENDIF
-#if defined key_lim2
-          IF(ln_bdy_ice_frs) THEN
-            frld_bdy (:) = frld_bdydta (:,2)
-            hicif_bdy(:) = hicif_bdydta(:,2)
-            hsnif_bdy(:) = hsnif_bdydta(:,2)
-          ENDIF
-#endif
-
-            IF(ln_bdy_tra_frs .or. ln_bdy_ice_frs) CALL iom_close( numbdyt )
-            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyu )
-            IF(ln_bdy_dyn_frs) CALL iom_close( numbdyv )
-         END IF
-
-      ENDIF                                            ! End if nit000
-
-
-      !                                                !---------------------!
-      !                                                !  at each time step  !
-      !                                                !---------------------!
-
-      IF( nbdy_dta == 1 .AND. ntimes_bdy > 1 ) THEN 
-         !
-         ! Read one more record if necessary
-         !**********************************
-
-        IF( ln_bdy_clim .AND. imois /= nbdy_b ) THEN      ! remember that nbdy_b=0 for kt=nit000
-           nbdy_b = imois
-           nbdy_a = imois + 1
-           nbdy_b = MOD( nbdy_b, iman )   ;   IF( nbdy_b == 0 ) nbdy_b = iman
-           nbdy_a = MOD( nbdy_a, iman )   ;   IF( nbdy_a == 0 ) nbdy_a = iman
-           lect=.true.
-        ELSEIF( .NOT.ln_bdy_clim .AND. itimer >= istep(nbdy_a) ) THEN
-
-           IF ( nbdy_a < ntimes_bdy ) THEN
-              nbdy_b = nbdy_a
-              nbdy_a = nbdy_a + 1
-              lect  =.true.
-           ELSE
-              ! We have reached the end of the file
-              ! put the last data time into both time levels
-              nbdy_b = nbdy_a
-              IF(ln_bdy_tra_frs) THEN
-                tbdydta(:,:,1) =  tbdydta(:,:,2)
-                sbdydta(:,:,1) =  sbdydta(:,:,2)
-              ENDIF
-              IF(ln_bdy_dyn_frs) THEN
-                ubdydta(:,:,1) =  ubdydta(:,:,2)
-                vbdydta(:,:,1) =  vbdydta(:,:,2)
-              ENDIF
-#if defined key_lim2
-              IF(ln_bdy_ice_frs) THEN
-                frld_bdydta (:,1) =  frld_bdydta (:,2)
-                hicif_bdydta(:,1) =  hicif_bdydta(:,2)
-                hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
-              ENDIF
+               IF(ln_bdy_ice_frs) THEN
+                  frld_bdydta (:,1) =  frld_bdydta (:,2)
+                  hicif_bdydta(:,1) =  hicif_bdydta(:,2)
+                  hsnif_bdydta(:,1) =  hsnif_bdydta(:,2)
+               ENDIF
 #endif
             END IF ! nbdy_a < ntimes_bdy
-
+            !
         END IF
          
-        IF( lect ) THEN
-           ! Swap arrays
-           IF(ln_bdy_tra_frs) THEN
+        IF( lect ) THEN           ! Swap arrays
+           IF( ln_bdy_tra_frs ) THEN
              tbdydta(:,:,1) =  tbdydta(:,:,2)
              sbdydta(:,:,1) =  sbdydta(:,:,2)
            ENDIF
-           IF(ln_bdy_dyn_frs) THEN
+           IF( ln_bdy_dyn_frs ) THEN
              ubdydta(:,:,1) =  ubdydta(:,:,2)
              vbdydta(:,:,1) =  vbdydta(:,:,2)
            ENDIF
 #if defined key_lim2
-           IF(ln_bdy_ice_frs) THEN
+           IF( ln_bdy_ice_frs ) THEN
              frld_bdydta (:,1) =  frld_bdydta (:,2)
              hicif_bdydta(:,1) =  hicif_bdydta(:,2)
@@ -582,5 +569,5 @@
            ipk  = jpk
 
-           IF(ln_bdy_tra_frs) THEN
+           IF( ln_bdy_tra_frs ) THEN
               ! 
               igrd = 1                                   ! temperature
@@ -720,5 +707,5 @@
       !
       ENDIF ! ln_bdy_dyn_frs .OR. ln_bdy_tra_frs
-
+      !
    END SUBROUTINE bdy_dta
 
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydyn.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdydyn.F90	(revision 2185)
@@ -36,9 +36,8 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -56,11 +55,11 @@
       INTEGER, INTENT( in ) ::   kt   ! Main time step counter
       !!
-      INTEGER  ::   ib, ik, igrd      ! dummy loop indices
-      INTEGER  ::   ii, ij            ! 2D addresses
-      REAL(wp) ::   zwgt              ! boundary weight
+      INTEGER  ::   jb, jk         ! dummy loop indices
+      INTEGER  ::   ii, ij, igrd   ! local integers
+      REAL(wp) ::   zwgt           ! boundary weight
       !!----------------------------------------------------------------------
       !
-      IF(ln_bdy_dyn_frs) THEN ! If this is false, then this routine does nothing. 
-
+      IF(ln_bdy_dyn_frs) THEN       ! If this is false, then this routine does nothing. 
+         !
          IF( kt == nit000 ) THEN
             IF(lwp) WRITE(numout,*)
@@ -70,33 +69,37 @@
          !
          igrd = 2                      ! Relaxation of zonal velocity
-         DO ib = 1, nblen(igrd)
-            DO ik = 1, jpkm1
-               ii = nbi(ib,igrd)
-               ij = nbj(ib,igrd)
-               zwgt = nbw(ib,igrd)
-               ua(ii,ij,ik) = ( ua(ii,ij,ik) * ( 1.- zwgt ) + ubdy(ib,ik) * zwgt ) * umask(ii,ij,ik)
+         DO jb = 1, nblen(igrd)
+            DO jk = 1, jpkm1
+               ii   = nbi(jb,igrd)
+               ij   = nbj(jb,igrd)
+               zwgt = nbw(jb,igrd)
+               ua(ii,ij,jk) = ( ua(ii,ij,jk) * ( 1.- zwgt ) + ubdy(jb,jk) * zwgt ) * umask(ii,ij,jk)
             END DO
          END DO
          !
          igrd = 3                      ! Relaxation of meridional velocity
-         DO ib = 1, nblen(igrd)
-            DO ik = 1, jpkm1
-               ii = nbi(ib,igrd)
-               ij = nbj(ib,igrd)
-               zwgt = nbw(ib,igrd)
-               va(ii,ij,ik) = ( va(ii,ij,ik) * ( 1.- zwgt ) + vbdy(ib,ik) * zwgt ) * vmask(ii,ij,ik)
+         DO jb = 1, nblen(igrd)
+            DO jk = 1, jpkm1
+               ii   = nbi(jb,igrd)
+               ij   = nbj(jb,igrd)
+               zwgt = nbw(jb,igrd)
+               va(ii,ij,jk) = ( va(ii,ij,jk) * ( 1.- zwgt ) + vbdy(jb,jk) * zwgt ) * vmask(ii,ij,jk)
             END DO
          END DO 
-         !
-         CALL lbc_lnk( ua, 'U', -1. )   ! Boundary points should be updated
-         CALL lbc_lnk( va, 'V', -1. )   !
+         CALL lbc_lnk( ua, 'U', -1. )   ;   CALL lbc_lnk( va, 'V', -1. )   ! Boundary points should be updated
          !
       ENDIF ! ln_bdy_dyn_frs
-
+      !
    END SUBROUTINE bdy_dyn_frs
 
 
-#if defined key_dynspg_exp || defined key_dynspg_ts
+# if defined   key_dynspg_exp   ||   defined key_dynspg_ts
+   !!----------------------------------------------------------------------
+   !!   'key_dynspg_exp'        OR              explicit sea surface height
+   !!   'key_dynspg_ts '                  split-explicit sea surface height
+   !!----------------------------------------------------------------------
+   
 !! Option to use Flather with dynspg_flt not coded yet...
+
    SUBROUTINE bdy_dyn_fla( pssh )
       !!----------------------------------------------------------------------
@@ -121,5 +124,5 @@
       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pssh
 
-      INTEGER  ::   ib, igrd                         ! dummy loop indices
+      INTEGER  ::   jb, igrd                         ! dummy loop indices
       INTEGER  ::   ii, ij, iim1, iip1, ijm1, ijp1   ! 2D addresses
       REAL(wp) ::   zcorr                            ! Flather correction
@@ -136,9 +139,9 @@
          igrd = 4
          spgu(:,:) = 0.0
-         DO ib = 1, nblenrim(igrd)
-            ii = nbi(ib,igrd)
-            ij = nbj(ib,igrd)
-            IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(ib)
-            IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(ib)
+         DO jb = 1, nblenrim(igrd)
+            ii = nbi(jb,igrd)
+            ij = nbj(jb,igrd)
+            IF( ln_bdy_dyn_fla ) spgu(ii, ij) = sshbdy(jb)
+            IF( ln_bdy_tides )   spgu(ii, ij) = spgu(ii, ij) + sshtide(jb)
          END DO
          !
@@ -146,12 +149,12 @@
          !             ! remember that flagu=-1 if normal velocity direction is outward
          !             ! I think we should rather use after ssh ?
-         DO ib = 1, nblenrim(igrd)
-            ii  = nbi(ib,igrd)
-            ij  = nbj(ib,igrd) 
-            iim1 = ii + MAX( 0, INT( flagu(ib) ) )   ! T pts i-indice inside the boundary
-            iip1 = ii - MIN( 0, INT( flagu(ib) ) )   ! T pts i-indice outside the boundary 
+         DO jb = 1, nblenrim(igrd)
+            ii  = nbi(jb,igrd)
+            ij  = nbj(jb,igrd) 
+            iim1 = ii + MAX( 0, INT( flagu(jb) ) )   ! T pts i-indice inside the boundary
+            iip1 = ii - MIN( 0, INT( flagu(jb) ) )   ! T pts i-indice outside the boundary 
             !
-            zcorr = - flagu(ib) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
-            zforc = ubtbdy(ib) + utide(ib)
+            zcorr = - flagu(jb) * SQRT( grav * hur_e(ii, ij) ) * ( pssh(iim1, ij) - spgu(iip1,ij) )
+            zforc = ubtbdy(jb) + utide(jb)
             ua_e(ii,ij) = zforc + zcorr * umask(ii,ij,1) 
          END DO
@@ -159,12 +162,12 @@
          igrd = 6      ! Flather bc on v-velocity
          !             ! remember that flagv=-1 if normal velocity direction is outward
-         DO ib = 1, nblenrim(igrd)
-            ii  = nbi(ib,igrd)
-            ij  = nbj(ib,igrd) 
-            ijm1 = ij + MAX( 0, INT( flagv(ib) ) )   ! T pts j-indice inside the boundary
-            ijp1 = ij - MIN( 0, INT( flagv(ib) ) )   ! T pts j-indice outside the boundary 
+         DO jb = 1, nblenrim(igrd)
+            ii  = nbi(jb,igrd)
+            ij  = nbj(jb,igrd) 
+            ijm1 = ij + MAX( 0, INT( flagv(jb) ) )   ! T pts j-indice inside the boundary
+            ijp1 = ij - MIN( 0, INT( flagv(jb) ) )   ! T pts j-indice outside the boundary 
             !
-            zcorr = - flagv(ib) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
-            zforc = vbtbdy(ib) + vtide(ib)
+            zcorr = - flagv(jb) * SQRT( grav * hvr_e(ii, ij) ) * ( pssh(ii, ijm1) - spgu(ii,ijp1) )
+            zforc = vbtbdy(jb) + vtide(jb)
             va_e(ii,ij) = zforc + zcorr * vmask(ii,ij,1)
          END DO
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyice.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 2185)
@@ -1,34 +1,32 @@
 MODULE bdyice
-   !!=================================================================================
+   !!======================================================================
    !!                       ***  MODULE  bdyice  ***
-   !! Ocean tracers:   Flow Relaxation Scheme of sea-ice fields on each open boundary
-   !!=================================================================================
-#if defined key_bdy && defined key_lim2
-   !!---------------------------------------------------------------------------------
-   !!   'key_bdy'      :                         Unstructured Open Boundary Conditions
-   !!---------------------------------------------------------------------------------
+   !! Unstructured Open Boundary Cond. :  Flow Relaxation Scheme applied sea-ice
+   !!======================================================================
+   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
+   !!----------------------------------------------------------------------
+#if defined   key_bdy   &&   defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_bdy'            and                 Unstructured Open Boundary Conditions
+   !!   'key_lim2'                                                 LIM-2 sea ice model
+   !!----------------------------------------------------------------------
    !!   bdy_ice        : Relaxation of tracers on unstructured open boundaries
-   !!---------------------------------------------------------------------------------
-   !! * Modules used
+   !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and tracers variables
-#if defined key_lim2
    USE ice_2           ! LIM_2 ice variables
-#endif
    USE dom_oce         ! ocean space and time domain variables 
    USE bdy_oce         ! ocean open boundary conditions
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE in_out_manager  ! write to numout file
+   
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC bdy_ice     ! routine called in stp
+   PUBLIC   bdy_ice    ! routine called in sbcmod
 
-   !! * Substitutions
-
-   !!---------------------------------------------------------------------------------
-   !!   OPA 9.0 , LODYC-IPSL  (2003)
-   !!---------------------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -37,59 +35,38 @@
       !!                 ***  SUBROUTINE bdy_ice  ***
       !!                    
-      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the  
-      !!              case of unstructured open boundaries. Currently only tested
-      !!              for LIM2.
+      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 
+      !!              of unstructured open boundaries. Currently only tested for LIM2.
       !! 
-      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in 
-      !!             a three-dimensional baroclinic ocean model with realistic
-      !!             topography. Tellus, 365-382.
-      !!  History :
-      !!   NEMO 3.3  !  2010-09 (D. Storkey) Original code
+      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
+      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
       !!------------------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt
-
-      !! * Local declarations
-      REAL(wp) :: zwgt                       ! boundary weight
-      INTEGER ::   jb, jk, jgrd              ! dummy loop indices
-      INTEGER ::   ii, ij                    ! 2D addresses
+      INTEGER, INTENT( in ) ::   kt   ! model time step index
+      !!
+      INTEGER  ::   jb, jk, jgrd   ! dummy loop indices
+      INTEGER  ::   ii, ij         ! local scalar
+      REAL(wp) ::   zwgt, zwgt1    ! local scalar
       !!------------------------------------------------------------------------------
-
-      jgrd=1 !: Everything is at T-points here
- 
-      IF(ln_bdy_ice_frs) THEN
- 
-        DO jb = 1, nblen(jgrd)
-          DO jk = 1, jpkm1
-            ii = nbi(jb,jgrd)
-            ij = nbj(jb,jgrd)
-            zwgt = nbw(jb,jgrd)
-
-            ! Leads fraction relaxation at the boundary   
-            frld(ii,ij) = ( frld(ii,ij)*(1.-zwgt)  +  frld_bdy(jb)*zwgt ) &
-                                                          * tmask(ii,ij,1)         
-
-            ! Ice depth relaxation at the boundary   
-            hicif(ii,ij) = ( hicif(ii,ij)*(1.-zwgt)  +  hicif_bdy(jb)*zwgt ) &
-                                                          * tmask(ii,ij,1)
-            ! Snow depth relaxation at the boundary   
-            hsnif(ii,ij) = ( hsnif(ii,ij)*(1.-zwgt)  +  hsnif_bdy(jb)*zwgt ) &
-                                                          * tmask(ii,ij,1)
-   
-          END DO
-        END DO 
-
-        CALL lbc_lnk( frld, 'T', 1. )  ! Boundary points should be updated
-        CALL lbc_lnk( hicif, 'T', 1. ) !
-        CALL lbc_lnk( hsnif, 'T', 1. ) !
-      
-      ELSE
-        ! we have called this routine without ln_bdy_ice_frs not set
-        IF(kt .EQ. nit000) THEN
-          WRITE(numout,*) 'E R R O R (possible) called bdy_ice when'
-          WRITE(numout,*) 'ln_bdy_ice_frs is false?'
-        ENDIF
-      ENDIF ! if ln_bdy_ice_frs
-      
+      !
+      jgrd = 1      ! Everything is at T-points here
+      !
+      IF( ln_bdy_ice_frs ) THEN     ! update ice fields by relaxation at the boundary
+         DO jb = 1, nblen(jgrd)
+            DO jk = 1, jpkm1
+               ii    = nbi(jb,jgrd)
+               ij    = nbj(jb,jgrd)
+               zwgt  = nbw(jb,jgrd)
+               zwgt1 = 1.e0 - nbw(jb,jgrd)
+               frld (ii,ij) = ( frld (ii,ij) * zwgt1 + frld_bdy (jb) * zwgt ) * tmask(ii,ij,1)     ! Leads fraction 
+               hicif(ii,ij) = ( hicif(ii,ij) * zwgt1 + hicif_bdy(jb) * zwgt ) * tmask(ii,ij,1)     ! Ice depth 
+               hsnif(ii,ij) = ( hsnif(ii,ij) * zwgt1 + hsnif_bdy(jb) * zwgt ) * tmask(ii,ij,1)     ! Snow depth
+            END DO
+         END DO 
+         CALL lbc_lnk( frld, 'T', 1. )                                         ! lateral boundary conditions
+         CALL lbc_lnk( hicif, 'T', 1. )   ;   CALL lbc_lnk( hsnif, 'T', 1. )
+         !
+      ELSE                          ! we have called this routine without ln_bdy_ice_frs not set
+         IF( kt == nit000 )   CALL ctl_warn( 'E R R O R (possible) called bdy_ice when ln_bdy_ice_frs is false?' )
+      ENDIF
+      !      
    END SUBROUTINE bdy_ice
 #else
@@ -99,4 +76,5 @@
 CONTAINS
    SUBROUTINE bdy_ice( kt )      ! Empty routine
+      WRITE(*,*) 'bdy_ice: You should not have seen this print! error?', kt
    END SUBROUTINE bdy_ice
 #endif
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 2185)
@@ -19,4 +19,5 @@
    USE oce             ! ocean dynamics and tracers variables
    USE dom_oce         ! ocean space and time domain
+   USE obc_par         ! ocean open boundary conditions
    USE bdy_oce         ! unstructured open boundary conditions
    USE bdytides        ! tides at open boundaries initialization (tide_init routine)
@@ -32,9 +33,8 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!---------------------------------------------------------------------------------
-
+   !!----------------------------------------------------------------------
 CONTAINS
    
@@ -50,5 +50,4 @@
       !!
       !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
-      !!
       !!----------------------------------------------------------------------      
       INTEGER ::   ii, ij, ik, igrd, ib, ir   ! dummy loop indices
@@ -56,6 +55,6 @@
       INTEGER ::   ib_len, ibr_max
       INTEGER ::   iw, ie, is, in 
-      INTEGER ::   inum                 ! temporary logical unit
-      INTEGER ::   id_dummy             ! temporary integers
+      INTEGER ::   inum                 ! local logical unit
+      INTEGER ::   id_dummy             ! local integers
       INTEGER ::   igrd_start, igrd_end ! start and end of loops on igrd
       INTEGER, DIMENSION (2)             ::   kdimsz
@@ -68,9 +67,8 @@
       !!
       NAMELIST/nambdy/filbdy_mask, filbdy_data_T, filbdy_data_U, filbdy_data_V,          &
-         &            filbdy_data_bt_T, filbdy_data_bt_U, filbdy_data_bt_V,          &
+         &            filbdy_data_bt_T, filbdy_data_bt_U, filbdy_data_bt_V,              &
          &            ln_bdy_tides, ln_bdy_clim, ln_bdy_vol, ln_bdy_mask,                &
          &            ln_bdy_dyn_fla, ln_bdy_dyn_frs, ln_bdy_tra_frs,ln_bdy_ice_frs,     &
          &            nbdy_dta, nb_rimwidth, volbdy
-
       !!----------------------------------------------------------------------
 
@@ -79,86 +77,71 @@
       IF(lwp) WRITE(numout,*) '~~~~~~~~'
       !
-      IF( jperio /= 0 ) CALL ctl_stop( 'Cyclic or symmetric,',   &
-           ' and unstructured open boundary condition are not compatible' )
-
-#if defined key_obc
-      CALL ctl_stop( 'Straight open boundaries,',   &
-           ' and unstructured open boundaries are not compatible' )
-#endif
-
-      ! Read namelist parameters
+      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
+         &                               ' and unstructured open boundary condition are not compatible' )
+
+      IF( lk_obc      )   CALL ctl_stop( 'Straight open boundaries,',   &
+         &                               ' and unstructured open boundaries are not compatible' )
+
       ! ---------------------------
-      REWIND( numnam )
+      REWIND( numnam )                    ! Read namelist parameters
       READ  ( numnam, nambdy )
 
-      ! control prints
+      !                                   ! control prints
       IF(lwp) WRITE(numout,*) '         nambdy'
 
-      ! Check nbdy_dta value
+      !                                         ! check type of data used (nbdy_dta value)
       IF(lwp) WRITE(numout,*) 'nbdy_dta =', nbdy_dta      
-      IF(lwp) WRITE(numout,*) ' '
-      SELECT CASE( nbdy_dta )
-      CASE( 0 )
-        IF(lwp) WRITE(numout,*) '         initial state used for bdy data'        
-      CASE( 1 )
-        IF(lwp) WRITE(numout,*) '         boundary data taken from file'
-      CASE DEFAULT
-        CALL ctl_stop( 'nbdy_dta must be 0 or 1' )
+      IF(lwp) WRITE(numout,*)
+      SELECT CASE( nbdy_dta )                   ! 
+      CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'        
+      CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
+      CASE DEFAULT   ;   CALL ctl_stop( 'nbdy_dta must be 0 or 1' )
       END SELECT
 
-      IF(lwp) WRITE(numout,*) ' '
+      IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS nb_rimwidth = ', nb_rimwidth
 
-      IF(lwp) WRITE(numout,*) ' '
-      IF(lwp) WRITE(numout,*) '         volbdy = ', volbdy
-
-      IF (ln_bdy_vol) THEN
-        SELECT CASE ( volbdy ) ! Check volbdy value
-        CASE( 1 )
-          IF(lwp) WRITE(numout,*) '         The total volume will be constant'
-        CASE( 0 )
-          IF(lwp) WRITE(numout,*) '         The total volume will vary according'
-          IF(lwp) WRITE(numout,*) '         to the surface E-P flux'
-        CASE DEFAULT
-          CALL ctl_stop( 'volbdy must be 0 or 1' )
-        END SELECT
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '      volbdy = ', volbdy
+
+      IF( ln_bdy_vol ) THEN                     ! check volume conservation (volbdy value)
+         SELECT CASE ( volbdy )
+         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
+         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
+         CASE DEFAULT   ;   CALL ctl_stop( 'volbdy must be 0 or 1' )
+         END SELECT
+         IF(lwp) WRITE(numout,*)
       ELSE
-        IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      IF (ln_bdy_tides) THEN
-        IF(lwp) WRITE(numout,*) ' '
+         IF(lwp) WRITE(numout,*) 'No volume correction with unstructured open boundaries'
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_tides ) THEN
         IF(lwp) WRITE(numout,*) 'Tidal harmonic forcing at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      IF (ln_bdy_dyn_fla) THEN
-        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_dyn_fla ) THEN
         IF(lwp) WRITE(numout,*) 'Flather condition on U, V at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      IF (ln_bdy_dyn_frs) THEN
-        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_dyn_frs ) THEN
         IF(lwp) WRITE(numout,*) 'FRS condition on U and V at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      IF (ln_bdy_tra_frs) THEN
-        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_tra_frs ) THEN
         IF(lwp) WRITE(numout,*) 'FRS condition on T & S fields at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      IF (ln_bdy_ice_frs) THEN
-        IF(lwp) WRITE(numout,*) ' '
+        IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_ice_frs ) THEN
         IF(lwp) WRITE(numout,*) 'FRS condition on ice fields at unstructured open boundaries'
-        IF(lwp) WRITE(numout,*) ' '
-      ENDIF
-
-      ! Read tides namelist 
-      ! ------------------------
-      IF( ln_bdy_tides )   CALL tide_init
+        IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ln_bdy_tides )   CALL tide_init      ! Read tides namelist 
+
 
       ! Read arrays defining unstructured open boundaries
@@ -170,8 +153,8 @@
       !          = 0  elsewhere   
  
-      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
+      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
          zmask(         :                ,:) = 0.e0
          zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0          
-      ELSE IF ( ln_bdy_mask ) THEN
+      ELSE IF( ln_bdy_mask ) THEN
          CALL iom_open( filbdy_mask, inum )
          CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
@@ -181,6 +164,5 @@
       ENDIF
 
-      ! Save mask over local domain      
-      DO ij = 1, nlcj
+      DO ij = 1, nlcj      ! Save mask over local domain      
          DO ii = 1, nlci
             bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
@@ -197,8 +179,6 @@
          END DO
       END DO
-
-      ! Lateral boundary conditions
-      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
-      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
+      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
+
 
       ! Read discrete distance and mapping indices
@@ -210,18 +190,16 @@
       IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
          icount = 0
-         ! Define west boundary (from ii=2 to ii=1+nb_rimwidth):
-         DO ir = 1, nb_rimwidth         
+         DO ir = 1, nb_rimwidth                  ! Define west boundary (from ii=2 to ii=1+nb_rimwidth):
             DO ij = 3, jpjglo-2
-               icount=icount+1
+               icount = icount + 1
                nbidta(icount,:) = ir + 1 + (jpizoom-1)
-               nbjdta(icount,:) = ij + (jpjzoom-1) 
+               nbjdta(icount,:) = ij     + (jpjzoom-1) 
                nbrdta(icount,:) = ir
             END DO
          END DO
-
-         ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth):
-         DO ir=1,nb_rimwidth         
+         !
+         DO ir = 1, nb_rimwidth                  ! Define east boundary (from ii=jpiglo-1 to ii=jpiglo-nb_rimwidth):
             DO ij=3,jpjglo-2
-               icount=icount+1
+               icount = icount + 1
                nbidta(icount,:) = jpiglo-ir + (jpizoom-1)
                nbidta(icount,2) = jpiglo-ir-1 + (jpizoom-1) ! special case for u points
@@ -230,15 +208,13 @@
             END DO
          END DO
-            
+         !       
       ELSE            ! Read indices and distances in unstructured boundary data files 
-
-         IF( ln_bdy_tides ) THEN 
-            ! Read tides input files for preference in case there are
-            ! no bdydata files. 
+         !
+         IF( ln_bdy_tides ) THEN             ! Read tides input files for preference in case there are no bdydata files
             clfile(4) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_T.nc'
             clfile(5) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_U.nc'
             clfile(6) = TRIM(filtide)//TRIM(tide_cpt(1))//'_grid_V.nc'
          ENDIF
-         IF( ln_bdy_dyn_fla .and. .not. ln_bdy_tides ) THEN 
+         IF( ln_bdy_dyn_fla .AND. .NOT. ln_bdy_tides ) THEN 
             clfile(4) = filbdy_data_bt_T
             clfile(5) = filbdy_data_bt_U
@@ -248,36 +224,26 @@
          IF( ln_bdy_tra_frs ) THEN 
             clfile(1) = filbdy_data_T
-            IF( .not. ln_bdy_dyn_frs ) THEN 
-               clfile(2) = filbdy_data_T  ! Dummy read re read T file for sake of 6 files
-               clfile(3) = filbdy_data_T  !
+            IF( .NOT. ln_bdy_dyn_frs ) THEN 
+               clfile(2) = filbdy_data_T     ! Dummy read re read T file for sake of 6 files
+               clfile(3) = filbdy_data_T     !
             ENDIF
          ENDIF          
          IF( ln_bdy_dyn_frs ) THEN 
-            IF( .not. ln_bdy_tra_frs ) THEN 
-               clfile(1) = filbdy_data_U ! Dummy Read 
-            ENDIF
+            IF( .NOT. ln_bdy_tra_frs )   clfile(1) = filbdy_data_U      ! Dummy Read 
             clfile(2) = filbdy_data_U
             clfile(3) = filbdy_data_V 
          ENDIF
 
-         ! how many files are we to read in?
-         IF(ln_bdy_tides .or. ln_bdy_dyn_fla) then
-             igrd_start = 4
+         !                                   ! how many files are we to read in?
+         IF(ln_bdy_tides .OR. ln_bdy_dyn_fla)   igrd_start = 4
+         !
+         IF(ln_bdy_tra_frs    ) THEN   ;   igrd_start = 1
+         ELSEIF(ln_bdy_dyn_frs) THEN   ;   igrd_start = 2
          ENDIF
-
-         IF(ln_bdy_tra_frs) then
-             igrd_start = 1
-         ELSEIF(ln_bdy_dyn_frs) then
-             igrd_start = 2
-         ENDIF
-
-         IF( ln_bdy_tra_frs ) then
-             igrd_end = 1
-         ENDIF
-
-         IF(ln_bdy_dyn_fla .or. ln_bdy_tides) THEN
-             igrd_end = 6
-         ELSEIF( ln_bdy_dyn_frs) THEN
-             igrd_end = 3
+         !
+         IF( ln_bdy_tra_frs   )   igrd_end = 1
+         !
+         IF(ln_bdy_dyn_fla .OR. ln_bdy_tides) THEN   ;   igrd_end = 6
+         ELSEIF( ln_bdy_dyn_frs             ) THEN   ;   igrd_end = 3
          ENDIF
 
@@ -287,8 +253,6 @@
             IF(lwp) WRITE(numout,*) 'kdimsz : ',kdimsz
             ib_len = kdimsz(1)
-            IF( ib_len > jpbdta) CALL ctl_stop(          &
-                'Boundary data array in file too long.', &
-                'File :', TRIM(clfile(igrd)),            &
-                'increase parameter jpbdta.' )
+            IF( ib_len > jpbdta)   CALL ctl_stop(  'Boundary data array in file too long.',                  &
+                &                                  'File :', TRIM(clfile(igrd)),'increase parameter jpbdta.' )
 
             CALL iom_get( inum, jpdom_unknown, 'nbidta', zdta(1:ib_len,:) )
@@ -298,20 +262,18 @@
             CALL iom_get( inum, jpdom_unknown, 'nbjdta', zdta(1:ib_len,:) )
             DO ii = 1,ib_len
-              nbjdta(ii,igrd) = INT( zdta(ii,1) )
-            END DO
-            CALL iom_get ( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )
+               nbjdta(ii,igrd) = INT( zdta(ii,1) )
+            END DO
+            CALL iom_get( inum, jpdom_unknown, 'nbrdta', zdta(1:ib_len,:) )
             DO ii = 1,ib_len
-              nbrdta(ii,igrd) = INT( zdta(ii,1) )
+               nbrdta(ii,igrd) = INT( zdta(ii,1) )
             END DO
             CALL iom_close( inum )
 
-            ! Check that rimwidth in file is big enough for Frs case(barotropic is one):
-            IF( igrd < 4) then
-                ibr_max = MAXVAL( nbrdta(:,igrd) )
-                IF(lwp) WRITE(numout,*)
-                IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
-                IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth
-                IF (ibr_max < nb_rimwidth) CALL ctl_stop( &
-                    'nb_rimwidth is larger than maximum rimwidth in file' )
+            IF( igrd < 4) THEN            ! Check that rimwidth in file is big enough for Frs case(barotropic is one):
+               ibr_max = MAXVAL( nbrdta(:,igrd) )
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
+               IF(lwp) WRITE(numout,*) ' nb_rimwidth from namelist is ', nb_rimwidth
+               IF (ibr_max < nb_rimwidth)   CALL ctl_stop( 'nb_rimwidth is larger than maximum rimwidth in file' )
             ENDIF !Check igrd < 4
             !
@@ -329,19 +291,19 @@
 
       DO igrd = igrd_start, igrd_end
-        icount  = 0
-        icountr = 0
-        nblen(igrd) = 0
-        nblenrim(igrd) = 0
-        nblendta(igrd) = 0
-        DO ir=1, nb_rimwidth
-          DO ib = 1, jpbdta
-          ! check if point is in local domain and equals ir
-            IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.   &
-               & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   &
-               & nbrdta(ib,igrd) == ir  ) THEN
-               !
-               icount = icount  + 1
-               !
-               IF( ir == 1 )   icountr = icountr+1
+         icount  = 0
+         icountr = 0
+         nblen   (igrd) = 0
+         nblenrim(igrd) = 0
+         nblendta(igrd) = 0
+         DO ir=1, nb_rimwidth
+            DO ib = 1, jpbdta
+               ! check if point is in local domain and equals ir
+               IF(  nbidta(ib,igrd) >= iw .AND. nbidta(ib,igrd) <= ie .AND.   &
+                  & nbjdta(ib,igrd) >= is .AND. nbjdta(ib,igrd) <= in .AND.   &
+                  & nbrdta(ib,igrd) == ir  ) THEN
+                  !
+                  icount = icount  + 1
+                  !
+                  IF( ir == 1 )   icountr = icountr+1
                   IF (icount > jpbdim) THEN
                      IF(lwp) WRITE(numout,*) 'bdy_ini: jpbdim too small'
@@ -364,10 +326,7 @@
       DO igrd = igrd_start, igrd_end
          DO ib = 1, nblen(igrd)
-            ! tanh formulation
-            nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )
-            ! quadratic
-!           nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2
-            ! linear
-!           nbw(ib,igrd) =  FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth)
+            nbw(ib,igrd) = 1.- TANH( FLOAT( nbr(ib,igrd) - 1 ) *0.5 )                     ! tanh formulation
+!           nbw(ib,igrd) = (FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth))**2      ! quadratic
+!           nbw(ib,igrd) =  FLOAT(nb_rimwidth+1-nbr(ib,igrd))/FLOAT(nb_rimwidth)          ! linear
          END DO
       END DO 
@@ -420,8 +379,6 @@
 
       ! Lateral boundary conditions
-      CALL lbc_lnk( fmask        , 'F', 1. )
-      CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
-      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )
-      CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
+      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
+      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
 
       IF( ln_bdy_vol .OR. ln_bdy_dyn_fla ) THEN      ! Indices and directions of rim velocity components
@@ -473,5 +430,4 @@
       ! Compute total lateral surface for volume correction:
       ! ----------------------------------------------------
- 
       bdysurftot = 0.e0 
       IF( ln_bdy_vol ) THEN  
@@ -491,5 +447,5 @@
                &                    * tmask_i(nbi(ib,igrd),nbj(ib,igrd)+1)
          END DO
-
+         !
          IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
       END IF   
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytides.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytides.F90	(revision 2185)
@@ -41,7 +41,6 @@
    PUBLIC   tide_update   ! routine called in bdydyn
 
-   LOGICAL, PUBLIC            ::   ln_tide_date            !: =T correct tide phases and amplitude for model start date
-
-   INTEGER, PARAMETER,PUBLIC  ::   jptides_max = 15      !: Max number of tidal contituents
+   LOGICAL, PUBLIC            ::   ln_tide_date          !: =T correct tide phases and amplitude for model start date
+   INTEGER, PUBLIC, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents
    INTEGER, PUBLIC            ::   ntide                 !: Actual number of tidal constituents
 
@@ -49,17 +48,16 @@
    CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used.
 
-   INTEGER , DIMENSION(jptides_max), PUBLIC ::   nindx        !: ???
-   REAL(wp), DIMENSION(jptides_max), PUBLIC ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
+   INTEGER , PUBLIC, DIMENSION(jptides_max) ::   nindx        !: ???
+   REAL(wp), PUBLIC, DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
    
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   !: Tidal constituents : SSH
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     !: Tidal constituents : U
-   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     !: Tidal constituents : V
+   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   ! Tidal constituents : SSH
+   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     ! Tidal constituents : U
+   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     ! Tidal constituents : V
    
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -88,11 +86,11 @@
       READ  ( numnam, nambdy_tide )
       !                                               ! Count number of components specified
-      ntide=jptides_max
-      do itide = 1, jptides_max
-        if ( tide_cpt(itide) == '' ) then
+      ntide = jptides_max
+      DO itide = 1, jptides_max
+        IF( tide_cpt(itide) == '' ) THEN
            ntide = itide-1
            exit
-        endif
-      enddo
+        ENDIF
+      END DO
 
       !                                               ! find constituents in standard list
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytra.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdytra.F90	(revision 2185)
@@ -25,9 +25,8 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -48,28 +47,27 @@
       !!----------------------------------------------------------------------
       !
-      IF(ln_bdy_tra_frs) THEN ! If this is false, then this routine does nothing. 
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'bdy_tra : Flow Relaxation Scheme for tracers'
-         IF(lwp) WRITE(numout,*) '~~~~~~~'
-      ENDIF
+      IF(ln_bdy_tra_frs) THEN       ! If this is false, then this routine does nothing. 
+         !
+         IF( kt == nit000 ) THEN
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'bdy_tra : Flow Relaxation Scheme for tracers'
+            IF(lwp) WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !
+         igrd = 1                       ! Everything is at T-points here
+         DO ib = 1, nblen(igrd)
+            DO ik = 1, jpkm1
+               ii = nbi(ib,igrd)
+               ij = nbj(ib,igrd)
+               zwgt = nbw(ib,igrd)
+               ta(ii,ij,ik) = ( ta(ii,ij,ik) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)         
+               sa(ii,ij,ik) = ( sa(ii,ij,ik) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)
+            END DO
+         END DO 
+         !
+         CALL lbc_lnk( ta, 'T', 1. )   ; CALL lbc_lnk( sa, 'T', 1. )    ! Boundary points should be updated
+         !
+      ENDIF ! ln_bdy_tra_frs
       !
-      igrd = 1                       ! Everything is at T-points here
-      DO ib = 1, nblen(igrd)
-         DO ik = 1, jpkm1
-            ii = nbi(ib,igrd)
-            ij = nbj(ib,igrd)
-            zwgt = nbw(ib,igrd)
-            ta(ii,ij,ik) = ( ta(ii,ij,ik) * (1.-zwgt) + tbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)         
-            sa(ii,ij,ik) = ( sa(ii,ij,ik) * (1.-zwgt) + sbdy(ib,ik) * zwgt ) * tmask(ii,ij,ik)
-        END DO
-      END DO 
-      !
-      CALL lbc_lnk( ta, 'T', 1. )   ! Boundary points should be updated
-      CALL lbc_lnk( sa, 'T', 1. )   !
-      !
-      ENDIF ! ln_bdy_tra_frs
-
    END SUBROUTINE bdy_tra
    
Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyvol.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 2184)
+++ branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyvol.F90	(revision 2185)
@@ -11,5 +11,5 @@
 #if   defined key_bdy   &&   defined key_dynspg_flt
    !!----------------------------------------------------------------------
-   !!   'key_bdy'            and      unstructured open boundary conditions
+   !!   'key_bdy'            AND      unstructured open boundary conditions
    !!   'key_dynspg_flt'                              filtered free surface
    !!----------------------------------------------------------------------
@@ -30,9 +30,8 @@
 #  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008) 
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -73,5 +72,5 @@
       INTEGER  ::   ji, jj, jk, jb, jgrd
       INTEGER  ::   ii, ij
-      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst, zraur
+      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
       !!-----------------------------------------------------------------------------
 
@@ -85,10 +84,9 @@
       ! -----------------------------------------------------------------------
       z_cflxemp = 0.e0
-      zraur = 1.e0 / rau0
-      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) * zraur ) 
-      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )   ! sum over the global domain
+      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * bdytmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
+      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
 
-      ! Barotropic velocity through the unstructured open boundary
-      ! ----------------------------------------------------------
+      ! Transport through the unstructured open boundary
+      ! ------------------------------------------------
       zubtpecor = 0.e0
       jgrd = 2                               ! cumulate u component contribution first 
@@ -112,6 +110,6 @@
       ! The normal velocity correction
       ! ------------------------------
-      IF (volbdy==1) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 
-      ELSE                  ;   zubtpecor =   zubtpecor             / bdysurftot
+      IF( volbdy==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / bdysurftot 
+      ELSE                   ;   zubtpecor =   zubtpecor             / bdysurftot
       END IF
 
@@ -141,5 +139,4 @@
       ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
       ! ------------------------------------------------------
-
       IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
          IF(lwp) WRITE(numout,*)
