C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_bp.F,v 1.3 2010/02/15 23:10:36 gforget Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"


      subroutine COST_BP(
     I                     myiter,
     I                     mytime,
     I                     mythid
     &                   )

c     ==================================================================
c     SUBROUTINE cost_bp
c     ==================================================================
c
c     o Evaluate cost function contribution of bottom pressure anoamlies 
c       => GRACE data
c
c     started: Gael Forget Oct-2009
c
c     ==================================================================
c     SUBROUTINE cost_bp
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"

#include "ecco_cost.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"
#include "DYNVARS.h"
#ifdef ALLOW_PROFILES
#include "profiles.h"
#endif

c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid

#ifdef ALLOW_BP_COST_CONTRIBUTION
#ifdef ALLOW_SMOOTH

c     == local variables ==

      integer bi,bj
      integer i,j
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer irec
      integer ilps

      logical doglobalread
      logical ladinit

      _RL bpdifmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL bpdifanom ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL bpdatmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL bpdatanom ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL bpcount ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL junk,junkweight

      character*(80) fname
      character*(80) fname4test
      character*(MAX_LEN_MBUF) msgbuf
      _RL fac

c     == external functions ==

      integer  ilnblnk
      external 

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--   Initialise local variables.
cgf convert phibot from m2/s2 to cm
      fac = 1. _d 2 / 9.81 _d 0

      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
              bpdifmean(i,j,bi,bj) = 0. _d 0
              bpdifanom(i,j,bi,bj) = 0. _d 0
              bpdatmean(i,j,bi,bj) = 0. _d 0
              bpdatanom(i,j,bi,bj) = 0. _d 0
              bpcount(i,j,bi,bj) = 0. _d 0
            enddo
          enddo
        enddo
      enddo

      doglobalread = .false.
      ladinit      = .false.

      write(fname(1:80),'(80a)') ' '
      ilps=ilnblnk( bpbarfile )
      write(fname(1:80),'(2a,i10.10)')
     &     bpbarfile(1:ilps),'.',optimcycle

c--   ============
c--   Mean values.
c--   ============

      do irec = 1, nmonsrec

c--     Compute the mean over all bpdat records.
        call ACTIVE_READ_XY( fname, bpbar, irec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_bpbar_mean_dummy )

        call COST_BP_READ( irec, mythid )

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
      if (bpmask(i,j,bi,bj).NE.0.) then
                bpdifmean(i,j,bi,bj) = bpdifmean(i,j,bi,bj) +
     &              ( fac*maskc(i,j,1,bi,bj)*bpbar(i,j,bi,bj) 
     &              - bpdat(i,j,bi,bj) )
                bpdatmean(i,j,bi,bj) = bpdatmean(i,j,bi,bj) +
     &                bpdat(i,j,bi,bj) 
                bpcount(i,j,bi,bj) = bpcount(i,j,bi,bj) +1.
      endif
              enddo
            enddo
          enddo
        enddo

      enddo

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
      if (bpcount(i,j,bi,bj).GT.0.) then
        bpdifmean(i,j,bi,bj) = bpdifmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
        bpdatmean(i,j,bi,bj) = bpdatmean(i,j,bi,bj)/bpcount(i,j,bi,bj)
      endif
              enddo
            enddo
          enddo
        enddo

c--   ==========
c--   Anomalies.
c--   ==========

c--   Loop over records for the second time.
      do irec = 1, nmonsrec

        call ACTIVE_READ_XY( fname, bpbar, irec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_bpbar_mean_dummy )

        call COST_BP_READ( irec, mythid )

c--    Compute field of anomalies
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
                bpdifanom(i,j,bi,bj) = 
     &                      ( fac*maskc(i,j,1,bi,bj)*bpbar(i,j,bi,bj) 
     &                      - bpdat(i,j,bi,bj) - bpdifmean(i,j,bi,bj)
     &                      )*bpmask(i,j,bi,bj)
                bpdatanom(i,j,bi,bj) =
     &                     ( bpdat(i,j,bi,bj) - bpdatmean(i,j,bi,bj)
     &                     )*bpmask(i,j,bi,bj)
              enddo
            enddo
          enddo
        enddo

      _EXCH_XY_RL ( bpdifanom, myThid )
      _EXCH_XY_RL ( bpdatanom, myThid )

c--    Smooth field of anomalies
      write(fname4test(1:80),'(1a)') 'bpdifanom_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,bpdifanom,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'bpdatanom_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,bpdatanom,irec,1,mythid)

      call SMOOTH_BASIC2D(bpdifanom,bpmask,750000. _d 0,5000,mythid)
      call SMOOTH_BASIC2D(bpdatanom,bpmask,750000. _d 0,5000,mythid)

      write(fname4test(1:80),'(1a)') 'bpdifanom_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,bpdifanom,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'bpdatanom_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,bpdatanom,irec,1,mythid)

c--    Compute cost function
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
c--             The array bpdat contains SSH anomalies.
                junkweight = wbp(i,j,bi,bj) *cosphi(i,j,bi,bj)
     &                      *bpmask(i,j,bi,bj)*maskc(i,j,1,bi,bj)
                junk       = bpdifanom(i,j,bi,bj)
                objf_bp(bi,bj) = objf_bp(bi,bj)
     &              + junk*junk*junkweight
                if ( junkweight .ne. 0. )
     &               num_bp(bi,bj) = num_bp(bi,bj) + 1. _d 0

              enddo
            enddo
          enddo
        enddo

      enddo

#endif /* ifdef ALLOW_SMOOTH */
#endif /* ifdef ALLOW_BP_COST_CONTRIBUTION */

      end