C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/cost/Attic/cost_ssh_mean.F,v 1.1.2.1 2002/02/05 20:23:57 heimbach Exp $

#include "COST_CPPOPTIONS.h"


      subroutine cost_ssh_mean(
     I                     psmean, offset
     O                   , costmean
     I                   , mythid
     &                   )

c     ==================================================================
c     SUBROUTINE cost_ssh_mean
c     ==================================================================
c
c     o Evaluate cost function contribution of sea surface height.
c       using of geoid error covariances requires regular model grid
c     o TODO: interpolate model grid to regular grid
c             check mask
c
c     started: Detlef Stammer, Ralf Giering Jul-1996
c
c     changed: Ralf Giering Ralf.Giering@FastOpt.de 12-Jun-2001
c
c              - totally rewrite for parallel processing
c
c     ==================================================================
c     SUBROUTINE cost_ssh_mean
c     ==================================================================

      implicit none

c     == global variables ==

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

#include "EESUPPORT.h"

#include "cost.h"

#ifdef ALLOW_EGM96_ERROR_COV
# include "cost_sph.h"
#endif

c     == routine arguments ==

      _RL     psmean ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL     offset
      _RL     costmean
      integer mythid

c     == local variables ==

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

#ifdef ALLOW_EGM96_ERROR_COV
      _RL cphi
      _RL xmean
      _RL shc( ncshc )
      _RL misfit  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy )
      _RL misfitgl( 1-olx:snx+olx, 1-oly:sny+oly, nsx,nsy, npx,npy )

      _RL pwork     ( (lshcmax+1)*(lshcmax+2)/2 )
      _RL diffearth ( lonshc, latshc )
      integer joffset
      integer ipx, ipy
      integer iglobe, jglobe
      integer iproc
      integer mpirc
      integer mpicrd(2)
#endif

      _RL diff
      Real*8 sumc

c     == end of interface ==

c--   Initialise local variables.
      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)

      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

      costmean   = 0.

#ifdef ALLOW_EGM96_ERROR_COV

c--   compute misfit on local domain
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = jmin,jmax
               if ( tpmeanmask(i,j,bi,bj) .gt. 0. ) then
                 misfit(i,j,bi,bj) =
     &           tpmean(i,j,bi,bj) - offset - psmean(i,j,bi,bj)
              else
                 misfit(i,j,bi,bj) = 0.
              endif
            enddo
          enddo
        enddo
      enddo

c--   communicate to get misfit on full 2d model surface grid
      call exch_all_2d_rl( misfit, misfitgl, mythid )

c--   set meridional index offset,
c--   it is non zero if the grid does no reach the poles
      joffset = (latshc - ny)/2

      _BEGIN_MASTER( mythid )
      print *
      print *, ' cost_SSH: grid starts at point ', joffset
     &       , ' of full earth.'
      _END_MASTER( mythid )

c--   preset field
c--   necessary if the grid does no reach the poles
      do j = 1, latshc
        do i = 1, lonshc
          diffearth(i,j) = 0.
        enddo
      enddo

c--   -interpolate- from model grid to regular grid
c--   so far we assume that the model grid is already regular
      _BEGIN_MASTER( mythid )

      do iproc = 1, nprocs

c--     get coordinates of processor (iporc-1)
        call MPI_Cart_coords(
     I          MPI_COMM_MODEL, iproc-1, 2, mpicrd
     O        , mpirc
     &        )

        ipx = 1 + mpicrd(1)
        ipy = 1 + mpicrd(2)

        do bj = jtlo,jthi
          do bi = itlo,ithi

            do j = jmin,jmax
              do i = jmin,jmax
                jglobe = joffset+ j + sNy*(bj-1) + sNy*nSy*(ipy-1)
                iglobe =          i + sNx*(bi-1) + sNx*nSx*(ipx-1)

                diffearth(iglobe,jglobe) = misfitgl(i,j,bi,bj,ipx,ipy)
              enddo
            enddo

          enddo
        enddo

      enddo

c--   Project regular grid misfit onto sperical harmonics
      call shc4grid(
     I               lshcmax
     O             , shc
     I             , latshc, lonshc
     I             , diffearth
     W             , pwork
     &             )

c--   Remove the C(0,0) component, i.e. the global mean.
      shc(1) = 0.

C--   Compute the cost function for the mean SSH.
      call cost_geoid(
     O                 costmean
     I               , shc
     I               , mythid
     &               )

      _END_MASTER( mythid )

      _BARRIER

#else

c--   Compute cost function for SSH by using the diagonal
c--   of the error covariance only.
c--
c--   Note: wp is assumed to include latitude dependence
c--         of error due to meridian convergence;
c--         --> no weighting by cosphi.

      sumc = 0.
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = jmin,jmax
              diff = psmean(i,j,bi,bj) - tpmean(i,j,bi,bj) - offset
              sumc = sumc + diff*diff
     &                    * wp(i,j,bi,bj)*tpmeanmask(i,j,bi,bj)
ce            --> check masks!
            enddo
          enddo
        enddo
      enddo

c--   Do the global summation.
      _GLOBAL_SUM_R8( sumc, mythid )
      costmean = sumc
#endif

      end