Index: branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyini.F90
===================================================================
--- branches/devukmo2010/NEMO/OPA_SRC/BDY/bdyini.F90	(revision 2128)
+++ 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   
