#include "COST_CPPOPTIONS.h"
#ifdef ALLOW_OBCS
# include "OBCS_OPTIONS.h"
#endif

      subroutine cost_obcss(
     I                       myiter,
     I                       mytime,
     I                       startrec,
     I                       endrec,
     I                       mythid
     &                     )

c     ==================================================================
c     SUBROUTINE cost_obcss
c     ==================================================================
c
c     o cost function contribution obc
c
c     o G. Gebbie, gebbie@mit.edu, 18-Mar-2003
c     ==================================================================
c     SUBROUTINE cost_obcss
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DYNVARS.h"
#ifdef ALLOW_OBCS
# include "OBCS.h"
#endif

#include "cal.h"
#include "ecco_cost.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"

c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid
cgg(  
      integer startrec
      integer endrec
cgg)

c     == local variables ==

      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec
      integer il
      integer iobcs
      integer jp1
      integer nrec
      integer ilfld
      integer igg

      _RL fctile
      _RL fcthread
      _RL dummy
      _RL gg
      _RL tmpx
      _RL tmpfield (1-olx:snx+olx,nr,nsx,nsy)
      _RL barofield (1-olx:snx+olx,nsx,nsy)
      _RL maskxz   (1-olx:snx+olx,nr,nsx,nsy)
      _RL vtop

      character*(80) fnamefld

      logical doglobalread
      logical ladinit

#ifdef ECCO_VERBOSE
      character*(MAX_LEN_MBUF) msgbuf
#endif

c     == external functions ==

      integer  ilnblnk
      external ilnblnk

c     == end of interface ==

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

c--   Read tiled data.
      doglobalread = .false.
      ladinit      = .false.

c     Number of records to be used.
      nrec = endrec-startrec+1

#ifdef ALLOW_OBCSS_COST_CONTRIBUTION

      jp1 = 1
      fcthread = 0. _d 0

#ifdef ECCO_VERBOSE
      _BEGIN_MASTER( mythid )
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a,i9.8)')
     &  ' cost_obcss: number of records to process: ',nrec
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      write(msgbuf,'(a)') ' '
      call print_message( msgbuf, standardmessageunit,
     &                    SQUEEZE_RIGHT , mythid)
      _END_MASTER( mythid )
#endif

      if (optimcycle .ge. 0) then
        ilfld=ilnblnk( xx_obcss_file )
        write(fnamefld(1:80),'(2a,i10.10)') 
     &       xx_obcss_file(1:ilfld), '.', optimcycle
      endif

c--   Loop over records.
      do irec = 1,nrec

         do bj = jtlo,jthi
           do bi = itlo,ithi
             do i = imin,imax
               barofield(i,bi,bj) = 0. _d 0
             enddo
           enddo
        enddo

        call active_read_xz( fnamefld, tmpfield, irec, doglobalread,
     &                       ladinit, optimcycle, mythid
     &        , xx_obcss_dummy )

cgg    Need to solve for iobcs would have been.
          gg    = (irec-1)/nobcs
          igg   = int(gg)
          iobcs = irec - igg*nobcs
cgg          print*,'S IREC, IOBCS',irec, iobcs

          call active_read_xz( 'maskobcss', maskxz, 
     &                         iobcs,
     &                         doglobalread, ladinit, 0,
     &                         mythid, dummy )

cgg     Need to change xx field to baroclinic vel.
cgg     Level Nr contains barotropic vel, remove it.
cgg     The deepest wet level velocity should be 
cgg     computed in order to ensure zero integrated flux.
        if (iobcs .eq. 3 .or. iobcs .eq. 4) then
c--     Loop over this thread's tiles.
          do bj = jtlo,jthi
            do bi = itlo,ithi
              do i = imin,imax

cgg         The barotropic velocity is stored in level 1.
cgg         Penalize for deviation in it.
                barofield(i,bi,bj) = tmpfield(i,1,bi,bj)
     &                             * maskxz(i,1,bi,bj)

                vtop = 0.d0
                tmpfield(i,1,bi,bj) = 0.d0

                do k = 1,Nr
cgg    If cells are not full, this should be modified with hFac.
cgg    
cgg    The xx field (tmpfldxz) does not contain the velocity at the 
cgg    lowest wet level. This velocity is not independent; it must 
cgg    exactly balance the volume flux, since we are dealing with
cgg    the baroclinic velocity structure.. 
                  vtop = tmpfield(i,k,bi,bj)*
     &               maskxz(i,k,bi,bj) * delZ(k) + vtop
                enddo
                tmpfield(i,1,bi,bj) = -vtop / delZ(1)
              enddo
            enddo
          enddo

        endif

c--     Loop over this thread's tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

c--         Determine the weights to be used.
            fctile = 0. _d 0

            if (iobcs .eq. 3 .or. iobcs .eq. 4) then
              do i = imin,imax
                j = OB_Js(I,bi,bj)
                tmpx = barofield(i,bi,bj)
                  fctile = fctile
     &               + wbaro*cosphi(i,j,bi,bj)
     &                 *tmpx*tmpx
     &                 *maskxz(i,1,bi,bj)
cgg                print*,'S bt fctile',fctile
              enddo
            endif

            do k = 1, Nr
              do i = imin,imax
                 j = OB_Js(I,bi,bj)
cgg                 if (maskS(i,j+jp1,k,bi,bj) .ne. 0.) then
                  tmpx = tmpfield(i,k,bi,bj)
                  fctile = fctile
     &                 + wobcssLev(i,k,bi,bj,iobcs)*cosphi(i,j,bi,bj)
     &                 *tmpx*tmpx
     &                 *maskxz(i,k,bi,bj)
cgg                endif
cgg                print*,'S fctile',fctile
              enddo
            enddo

            objf_obcss(bi,bj) = objf_obcss(bi,bj) + fctile
            fcthread         = fcthread + fctile
          enddo
        enddo

#ifdef ECCO_VERBOSE
c--     Print cost function for all tiles.
        _GLOBAL_SUM_R8( fcthread , myThid )
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,i8.8)')
     &    ' cost_obcss: irec = ',irec
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,a,d22.15)')
     &    ' global cost function value',
     &    ' (obcss) = ',fcthread
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a)') ' '
        call print_message( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)
#endif

      enddo
c--   End of loop over records.

#endif

      return
      end