Index: branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90
===================================================================
--- branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90	(revision 2796)
+++ branches/2011/UKMO_MERCATOR_obc_bdy_merge/NEMOGCM/NEMO/OPA_SRC/OBC/obcvol.F90	(revision 2797)
@@ -1,288 +1,167 @@
 MODULE obcvol
-   !!=================================================================================
+   !!======================================================================
    !!                       ***  MODULE  obcvol  ***
-   !! Ocean dynamic :  Volume constraint when OBC and Free surface are used
-   !!=================================================================================
-#if   defined key_obc   &&   ! defined key_vvl
-   !!---------------------------------------------------------------------------------
-   !!   'key_obc'               and   NOT                 open boundary conditions
-   !!   'key_vvl'                                         constant volume free surface
-   !!---------------------------------------------------------------------------------
-   !! * Modules used
-   USE oce             ! ocean dynamics and tracers 
-   USE dom_oce         ! ocean space and time domain 
-   USE sbc_oce         ! ocean surface boundary conditions
-   USE phycst          ! physical constants
-   USE obc_oce         ! ocean open boundary conditions
-   USE lib_mpp         ! for mppsum
-   USE in_out_manager  ! I/O manager
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Accessibility
-   PUBLIC obc_vol        ! routine called by dynspg_flt
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "obc_vectopt_loop_substitute.h90"
-   !!---------------------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!---------------------------------------------------------------------------------
-
+   !! Ocean dynamic :  Volume constraint when unstructured boundary 
+   !!                  and Free surface are used
+   !!======================================================================
+   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
+   !!             -   !  2006-01  (J. Chanut) Bug correction
+   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
+   !!----------------------------------------------------------------------
+#if   defined key_obc   &&   defined key_dynspg_flt
+!!$   !!----------------------------------------------------------------------
+!!$   !!   'key_obc'            AND      unstructured open boundary conditions
+!!$   !!   'key_dynspg_flt'                              filtered free surface
+!!$   !!----------------------------------------------------------------------
+!!$   USE oce             ! ocean dynamics and tracers 
+!!$   USE dom_oce         ! ocean space and time domain 
+!!$   USE phycst          ! physical constants
+!!$   USE obc_oce         ! ocean open boundary conditions
+!!$   USE lib_mpp         ! for mppsum
+!!$   USE in_out_manager  ! I/O manager
+!!$   USE sbc_oce         ! ocean surface boundary conditions
+!!$
+!!$   IMPLICIT NONE
+!!$   PRIVATE
+!!$
+!!$   PUBLIC obc_vol        ! routine called by dynspg_flt.h90
+!!$
+!!$   !! * Substitutions
+!!$#  include "domzgr_substitute.h90"
+!!$   !!----------------------------------------------------------------------
+!!$   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+!!$   !! $Id$ 
+!!$   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+!!$   !!----------------------------------------------------------------------
+!!$CONTAINS
+!!$
+!!$   SUBROUTINE obc_vol( kt )
+!!$      !!----------------------------------------------------------------------
+!!$      !!                      ***  ROUTINE obcvol  ***
+!!$      !!
+!!$      !! ** Purpose :   This routine is called in dynspg_flt to control 
+!!$      !!      the volume of the system. A correction velocity is calculated
+!!$      !!      to correct the total transport through the unstructured OBC. 
+!!$      !!      The total depth used is constant (H0) to be consistent with the 
+!!$      !!      linear free surface coded in OPA 8.2
+!!$      !!
+!!$      !! ** Method  :   The correction velocity (zubtpecor here) is defined calculating
+!!$      !!      the total transport through all open boundaries (trans_obc) minus
+!!$      !!      the cumulate E-P flux (z_cflxemp) divided by the total lateral 
+!!$      !!      surface (obcsurftot) of the unstructured boundary. 
+!!$      !!         zubtpecor = [trans_obc - z_cflxemp ]*(1./obcsurftot)
+!!$      !!      with z_cflxemp => sum of (Evaporation minus Precipitation)
+!!$      !!                       over all the domain in m3/s at each time step.
+!!$      !!      z_cflxemp < 0 when precipitation dominate
+!!$      !!      z_cflxemp > 0 when evaporation dominate
+!!$      !!
+!!$      !!      There are 2 options (user's desiderata): 
+!!$      !!         1/ The volume changes according to E-P, this is the default
+!!$      !!            option. In this case the cumulate E-P flux are setting to
+!!$      !!            zero (z_cflxemp=0) to calculate the correction velocity. So
+!!$      !!            it will only balance the flux through open boundaries.
+!!$      !!            (set nn_volctl to 0 in tne namelist for this option)
+!!$      !!         2/ The volume is constant even with E-P flux. In this case
+!!$      !!            the correction velocity must balance both the flux 
+!!$      !!            through open boundaries and the ones through the free
+!!$      !!            surface. 
+!!$      !!            (set nn_volctl to 1 in tne namelist for this option)
+!!$      !!----------------------------------------------------------------------
+!!$      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+!!$      !!
+!!$      INTEGER  ::   ji, jj, jk, jb, jgrd
+!!$      INTEGER  ::   ii, ij
+!!$      REAL(wp) ::   zubtpecor, z_cflxemp, ztranst
+!!$      !!-----------------------------------------------------------------------------
+!!$
+!!$      IF( ln_vol ) THEN
+!!$
+!!$      IF( kt == nit000 ) THEN 
+!!$         IF(lwp) WRITE(numout,*)
+!!$         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along unstructured OBC'
+!!$         IF(lwp) WRITE(numout,*)'~~~~~~~'
+!!$      END IF 
+!!$
+!!$      ! Calculate the cumulate surface Flux z_cflxemp (m3/s) over all the domain
+!!$      ! -----------------------------------------------------------------------
+!!$      z_cflxemp = SUM ( ( emp(:,:)-rnf(:,:) ) * obctmask(:,:) * e1t(:,:) * e2t(:,:) ) / rau0
+!!$      IF( lk_mpp )   CALL mpp_sum( z_cflxemp )     ! sum over the global domain
+!!$
+!!$      ! Transport through the unstructured open boundary
+!!$      ! ------------------------------------------------
+!!$      zubtpecor = 0.e0
+!!$      jgrd = 2                               ! cumulate u component contribution first 
+!!$      DO jb = 1, nblenrim(jgrd)
+!!$         DO jk = 1, jpkm1
+!!$            ii = nbi(jb,jgrd)
+!!$            ij = nbj(jb,jgrd)
+!!$            zubtpecor = zubtpecor + flagu(jb) * ua(ii,ij, jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+!!$         END DO
+!!$      END DO
+!!$      jgrd = 3                               ! then add v component contribution
+!!$      DO jb = 1, nblenrim(jgrd)
+!!$         DO jk = 1, jpkm1
+!!$            ii = nbi(jb,jgrd)
+!!$            ij = nbj(jb,jgrd)
+!!$            zubtpecor = zubtpecor + flagv(jb) * va(ii,ij, jk) * e1v(ii,ij) * fse3v(ii,ij,jk) 
+!!$         END DO
+!!$      END DO
+!!$      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
+!!$
+!!$      ! The normal velocity correction
+!!$      ! ------------------------------
+!!$      IF( nn_volctl==1 ) THEN   ;   zubtpecor = ( zubtpecor - z_cflxemp) / obcsurftot 
+!!$      ELSE                   ;   zubtpecor =   zubtpecor             / obcsurftot
+!!$      END IF
+!!$
+!!$      ! Correction of the total velocity on the unstructured boundary to respect the mass flux conservation
+!!$      ! -------------------------------------------------------------
+!!$      ztranst = 0.e0
+!!$      jgrd = 2                               ! correct u component
+!!$      DO jb = 1, nblenrim(jgrd)
+!!$         DO jk = 1, jpkm1
+!!$            ii = nbi(jb,jgrd)
+!!$            ij = nbj(jb,jgrd)
+!!$            ua(ii,ij,jk) = ua(ii,ij,jk) - flagu(jb) * zubtpecor * umask(ii,ij,jk)
+!!$            ztranst = ztranst + flagu(jb) * ua(ii,ij,jk) * e2u(ii,ij) * fse3u(ii,ij,jk)
+!!$         END DO
+!!$      END DO
+!!$      jgrd = 3                              ! correct v component
+!!$      DO jb = 1, nblenrim(jgrd)
+!!$         DO jk = 1, jpkm1
+!!$            ii = nbi(jb,jgrd)
+!!$            ij = nbj(jb,jgrd)
+!!$            va(ii,ij,jk) = va(ii,ij,jk) -flagv(jb) * zubtpecor * vmask(ii,ij,jk)
+!!$            ztranst = ztranst + flagv(jb) * va(ii,ij,jk) * e1v(ii,ij) * fse3v(ii,ij,jk)
+!!$         END DO
+!!$      END DO
+!!$      IF( lk_mpp )   CALL mpp_sum( ztranst )   ! sum over the global domain
+!!$ 
+!!$      ! Check the cumulated transport through unstructured OBC once barotropic velocities corrected
+!!$      ! ------------------------------------------------------
+!!$      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
+!!$         IF(lwp) WRITE(numout,*)
+!!$         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
+!!$         IF(lwp) WRITE(numout,*)'~~~~~~~ '
+!!$         IF(lwp) WRITE(numout,*)'          cumulate flux EMP             =', z_cflxemp  , ' (m3/s)'
+!!$         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC  =', obcsurftot, '(m2)'
+!!$         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor =', zubtpecor , '(m/s)'
+!!$         IF(lwp) WRITE(numout,*)'          cumulated transport ztranst   =', ztranst   , '(m3/s)'
+!!$      END IF 
+!!$      !
+!!$      END IF ! ln_vol
+!!$
+!!$   END SUBROUTINE obc_vol
+!!$
+!!$#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                   NO Unstruct Open Boundary Conditions
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE obc_vol ( kt )
-      !!------------------------------------------------------------------------------
-      !!                      ***  ROUTINE obcvol  ***
-      !!
-      !! ** Purpose : 
-      !!      This routine is called in dynspg_flt to control 
-      !!      the volume of the system. A correction velocity is calculated
-      !!      to correct the total transport through the OBC. 
-      !!      The total depth used is constant (H0) to be consistent with the 
-      !!      linear free surface coded in OPA 8.2
-      !!
-      !! ** Method :  
-      !!      The correction velocity (zubtpecor here) is defined calculating
-      !!      the total transport through all open boundaries (trans_obc) minus
-      !!      the cumulate E-P flux (zCflxemp) divided by the total lateral 
-      !!      surface (obcsurftot) of these OBC. 
-      !!
-      !!      zubtpecor = [trans_obc - zCflxemp ]*(1./obcsurftot)
-      !!
-      !!      with zCflxemp => sum of (Evaporation minus Precipitation)
-      !!                       over all the domain in m3/s at each time step.
-      !!
-      !!      zCflxemp < 0 when precipitation dominate
-      !!      zCflxemp > 0 when evaporation dominate
-      !!
-      !!      There are 2 options (user's desiderata): 
-      !!
-      !!         1/ The volume changes according to E-P, this is the default
-      !!            option. In this case the cumulate E-P flux are setting to
-      !!            zero (zCflxemp=0) to calculate the correction velocity. So
-      !!            it will only balance the flux through open boundaries.
-      !!            (set volemp to 0 in tne namelist for this option)
-      !!
-      !!         2/ The volume is constant even with E-P flux. In this case
-      !!            the correction velocity must balance both the flux 
-      !!            through open boundaries and the ones through the free
-      !!            surface. 
-      !!            (set volemp to 1 in tne namelist for this option)
-      !!
-      !! History :
-      !!   8.5  !  02-10 (C. Talandier, A-M. Treguier) Original code
-      !!----------------------------------------------------------------------------
-      !! * Arguments
-      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
-
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk
-      REAL(wp) ::   zubtpecor
-      REAL(wp) ::   zCflxemp
-      REAL(wp) ::   ztransw, ztranse, ztransn, ztranss, ztranst
-      !!-----------------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN 
-         IF(lwp) WRITE(numout,*)'        '
-         IF(lwp) WRITE(numout,*)'obc_vol : Correction of velocities along OBC'
-         IF(lwp) WRITE(numout,*)'~~~~~~~'
-         IF(lwp) WRITE(numout,*)'        '
-      END IF 
-
-      ! 1. Calculate the cumulate surface Flux zCflxemp (m3/s) over all the domain.
-      ! ---------------------------------------------------------------------------
-
-      zCflxemp = SUM ( ( emp(:,:)-rnf(:,:) )*obctmsk(:,:)* e1t(:,:) * e2t(:,:)  / rau0 ) 
-
-      IF( lk_mpp )   CALL mpp_sum( zCflxemp )   ! sum over the global domain
-
-      ! 2. Barotropic velocity for each open boundary
-      ! ---------------------------------------------
-
-      zubtpecor = 0.e0
-
-      ! ... East open boundary
-      IF( lp_obc_east ) THEN                      ! ... Total transport through the East OBC
-         DO ji = fs_nie0, fs_nie1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  zubtpecor = zubtpecor - ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * &
-             &     uemsk(jj,jk)*MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
-               END DO
-            END DO
-         END DO
-      END IF 
-
-      ! ... West open boundary
-      IF( lp_obc_west ) THEN                      ! ... Total transport through the West OBC
-         DO ji = fs_niw0, fs_niw1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  zubtpecor = zubtpecor + ua(ji,jj,jk)*e2u(ji,jj)*fse3u(ji,jj,jk) * &
-             &    uwmsk(jj,jk) *MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
-               END DO
-            END DO
-         END DO
-       ENDIF
-
-      ! ... North open boundary
-      IF( lp_obc_north ) THEN                     ! ... Total transport through the North OBC
-         DO jj = fs_njn0, fs_njn1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zubtpecor = zubtpecor - va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * &
-             &    vnmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
-               END DO
-            END DO
-         END DO
-       ENDIF
-
-      ! ... South open boundary
-      IF( lp_obc_south ) THEN                     ! ... Total transport through the South OBC
-         DO jj = fs_njs0, fs_njs1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO ji = 1, jpi
-                  zubtpecor = zubtpecor + va(ji,jj,jk)*e1v(ji,jj)*fse3v(ji,jj,jk) * &
-             &    vsmsk(ji,jk) * MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
-               END DO
-            END DO
-         END DO
-       ENDIF
-
-      IF( lk_mpp )   CALL mpp_sum( zubtpecor )   ! sum over the global domain
-
-
-      ! 3. The normal velocity correction
-      ! ---------------------------------
-      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-         IF(lwp) WRITE(numout,*)'        '
-         IF(lwp) WRITE(numout,*)'obc_vol : time step :', kt
-         IF(lwp) WRITE(numout,*)'~~~~~~~ '
-         IF(lwp) WRITE(numout,*)'          cumulate flux EMP :', zCflxemp,' (m3/s)'
-         IF(lwp) WRITE(numout,*)'          lateral transport :',zubtpecor,'(m3/s)'
-         IF(lwp) WRITE(numout,*)'          net inflow        :',zubtpecor-zCflxemp,'(m3/s)'
-      ENDIF
-
-      zubtpecor = (zubtpecor - zCflxemp*volemp)*(1./obcsurftot)
-
-      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-         IF(lwp) WRITE(numout,*)'          total lateral surface of OBC :',obcsurftot,'(m2)'
-         IF(lwp) WRITE(numout,*)'          correction velocity zubtpecor :',zubtpecor,'(m/s)'
-         IF(lwp) WRITE(numout,*)'        '
-      END IF 
-
-      ! 4. Correction of the total velocity on each open 
-      !    boundary to respect the mass flux conservation
-      ! -------------------------------------------------
-
-      ztranse = 0.e0   ; ztransw = 0.e0 ; ztransn = 0.e0 ; ztranss = 0.e0
-      ztranst = 0.e0  ! total
-
-      IF( lp_obc_west ) THEN
-         ! ... correction of the west velocity
-         DO ji = fs_niw0, fs_niw1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  ua(ji,jj,jk) = ua(ji,jj,jk) - zubtpecor*uwmsk(jj,jk)
-                  ztransw= ztransw + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uwmsk(jj,jk) * &
-             &    MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
-               END DO
-            END DO
-         END DO
-
-         IF( lk_mpp )   CALL mpp_sum( ztransw )   ! sum over the global domain
-
-         IF( lwp .AND. MOD( kt, nwrite ) == 0)  WRITE(numout,*)'          West OB transport ztransw :', ztransw,'(m3/s)'
-      END IF 
-
-      IF( lp_obc_east ) THEN
-
-         ! ... correction of the east velocity
-         DO ji = fs_nie0, fs_nie1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO jj = 1, jpj
-                  ua(ji,jj,jk) = ua(ji,jj,jk) + zubtpecor*uemsk(jj,jk)
-                  ztranse= ztranse + ua(ji,jj,jk)*fse3u(ji,jj,jk)*e2u(ji,jj)*uemsk(jj,jk) * &
-            &     MAX(obctmsk(ji,jj),obctmsk(ji+1,jj) )
-               END DO
-            END DO
-         END DO
-
-         IF( lk_mpp )   CALL mpp_sum( ztranse )   ! sum over the global domain
-
-         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-            IF(lwp) WRITE(numout,*)'          East OB transport ztranse :', ztranse,'(m3/s)'
-         END IF 
-
-      END IF 
-
-      IF( lp_obc_north ) THEN
-
-         ! ... correction of the north velocity
-         DO jj = fs_njn0, fs_njn1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO ji =  1, jpi
-                  va(ji,jj,jk) = va(ji,jj,jk) + zubtpecor*vnmsk(ji,jk)
-                  ztransn= ztransn + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vnmsk(ji,jk) * &
-           &      MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztransn )   ! sum over the global domain
-
-         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-            IF(lwp) WRITE(numout,*)'          North OB transport ztransn :', ztransn,'(m3/s)'
-         END IF 
-
-      END IF 
-
-      IF( lp_obc_south ) THEN
-
-         ! ... correction of the south velocity
-         DO jj = fs_njs0, fs_njs1 ! Vector opt.
-            DO jk = 1, jpkm1
-               DO ji =  1, jpi
-                  va(ji,jj,jk) = va(ji,jj,jk) - zubtpecor*vsmsk(ji,jk)
-                  ztranss= ztranss + va(ji,jj,jk)*fse3v(ji,jj,jk)*e1v(ji,jj)*vsmsk(ji,jk) * &
-            &     MAX(obctmsk(ji,jj),obctmsk(ji,jj+1) )
-               END DO
-            END DO
-         END DO
-         IF( lk_mpp )   CALL mpp_sum( ztranss )   ! sum over the global domain
-
-         IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-            IF(lwp) WRITE(numout,*)'          South OB transport ztranss :', ztranss,'(m3/s)'
-         END IF 
-
-      END IF 
-
-      ! 5. Check the cumulate transport through OBC
-      !    once barotropic velocities corrected
-      ! -------------------------------------------
-
-
-      IF( lwp .AND. MOD( kt, nwrite ) == 0) THEN
-         ztranst = ztransw - ztranse + ztranss - ztransn
-         IF(lwp) WRITE(numout,*)'        '
-         IF(lwp) WRITE(numout,*)'          Cumulate transport ztranst =', ztranst,'(m3/s)'
-         IF(lwp) WRITE(numout,*)'          Balance  =', ztranst - zCflxemp ,'(m3/s)'
-         IF(lwp) WRITE(numout,*)'        '
-      END IF 
-
-   END SUBROUTINE obc_vol
-
-#else
-   !!---------------------------------------------------------------------------------
-   !!  Default option :                                                   Empty module
-   !!---------------------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE obc_vol        ! Empty routine
+   SUBROUTINE obc_vol( kt )        ! Empty routine
+      WRITE(*,*) 'obc_vol: You should not have seen this print! error?', kt
    END SUBROUTINE obc_vol
 #endif
 
-   !!=================================================================================
+   !!======================================================================
 END MODULE obcvol
