C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/SOSE/code_ad/cost_ssh_new.F,v 1.1 2010/04/23 19:55:11 mmazloff Exp $
C $Name:  $

#include "COST_CPPOPTIONS.h"


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

c     ==================================================================
c     SUBROUTINE cost_ssh
c     ==================================================================
c
c     o Evaluate cost function contribution of sea surface height.
c       using of geoid error covariances requires regular model grid
c
c     started: Detlef Stammer, Ralf Giering Jul-1996
c
c     changed: Christian Eckert eckert@mit.edu 25-Feb-2000
c
c              - Restructured the code in order to create a package
c                for the MITgcmUV.
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
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.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_SSH_COST_CONTRIBUTION
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
      integer gwunit

      logical doglobalread
      logical ladinit

      _RL offset
      _RL costmean
      _RL offset_sum
      _RL psmean    ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL psmeantp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL psmeaners ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL psmeangfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL sumtp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL sumers ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL sumgfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )

      _RL wwwtp1  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL wwwers1 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL wwwgfo1 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL wwwtp2  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL wwwers2 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL wwwgfo2 ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL junk,diff
CMM(
      _RL obsmeantp  ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL obsmeaners ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL obsmeangfo ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL tpmean_loc ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL mean_slaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMssh_all(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMssh_all_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMssh_all_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL diagnosfld(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL sumc, num_hmean_mm
      _RL factor, spval
      integer IDOT, k
      character*(80) fname4test
CMM)
      character*(80) fname
      character*(MAX_LEN_MBUF) msgbuf

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--   Initialise local variables.
      costmean   = 0. _d 0

c--   First, read tiled data.
      doglobalread = .false.
      ladinit      = .false.

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

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

      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
              psmean(i,j,bi,bj)    = 0. _d 0
              psmeantp(i,j,bi,bj)  = 0. _d 0
              psmeaners(i,j,bi,bj) = 0. _d 0
              psmeangfo(i,j,bi,bj) = 0. _d 0
CMM(  initialize
              obsmeantp(i,j,bi,bj) = 0. _d 0
              obsmeaners(i,j,bi,bj) = 0. _d 0
              obsmeangfo(i,j,bi,bj) = 0. _d 0 
              mean_slaobs(i,j,bi,bj)  = 0. _d 0
              mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0 
              mean_psMssh_all(i,j,bi,bj) = 0. _d 0
              mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
              mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
              diagnosfld(i,j,bi,bj)  = 0. _d 0
              tpmean_loc(i,j,bi,bj) = 0. _d 0
CMM)
              sumtp(i,j,bi,bj)  = 0. _d 0
              sumers(i,j,bi,bj) = 0. _d 0
              sumgfo(i,j,bi,bj) = 0. _d 0
              wwwtp1(i,j,bi,bj)  = 0. _d 0
              wwwers1(i,j,bi,bj) = 0. _d 0
              wwwgfo1(i,j,bi,bj) = 0. _d 0
              wwwtp2(i,j,bi,bj)  = 0. _d 0
              wwwers2(i,j,bi,bj) = 0. _d 0
              wwwgfo2(i,j,bi,bj) = 0. _d 0
            enddo
          enddo
        enddo
      enddo

c--   Read mean field and generate mask
      call cost_ReadTopexMean( mythid )

CMM(  HAD TO MOVE OFFSET CALC TO GET IT INTO RECORD LOOP
c--   Compute and remove offset for current tile and sum over all
c--   tiles of this instance.
      offset     = 0. _d 0
      offset_sum = 0. _d 0
CMM)
c--   Loop over records for the first time.
      do irec = 1, ndaysrec

c--     Compute the mean over all psbar records.
        call active_read_xy( fname, psbar, irec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_psbar_mean_dummy )

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
cph these if would be more efficient, but need recomputations
c        if ( tpTimeMask(irec) .NE. 0. )
        call cost_readtopex( irec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
cph these if would be more efficient, but need recomputations
c        if ( ersTimeMask(irec) .NE. 0. )
        call cost_readers( irec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
cph these if would be more efficient, but need recomputations
c        if ( gfoTimeMask(irec) .NE. 0. )
        call cost_readgfo( irec, mythid )
#endif

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
                psmean(i,j,bi,bj) = psmean(i,j,bi,bj) +
     &                psbar(i,j,bi,bj) / float(ndaysrec)
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
                if ( tpmask(i,j,bi,bj)*wtp(i,j,bi,bj)
     &               *tpTimeMask(irec) .NE. 0. ) then
                   psmeantp(i,j,bi,bj) = psmeantp(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)
                   obsmeantp(i,j,bi,bj) = obsmeantp(i,j,bi,bj) +
     &                 tpobs(i,j,bi,bj)
                   sumtp(i,j,bi,bj) = sumtp(i,j,bi,bj) + 1. _d 0
                endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
                if ( ersmask(i,j,bi,bj)*wers(i,j,bi,bj)
     &               *ersTimeMask(irec) .NE. 0. ) then
                   psmeaners(i,j,bi,bj) = psmeaners(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)
                   obsmeaners(i,j,bi,bj) = obsmeaners(i,j,bi,bj) +
     &                 ersobs(i,j,bi,bj)
                   sumers(i,j,bi,bj) = sumers(i,j,bi,bj) + 1. _d 0
                endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
                if ( gfomask(i,j,bi,bj)*wgfo(i,j,bi,bj)
     &               *gfoTimeMask(irec) .NE. 0. ) then
                   psmeangfo(i,j,bi,bj) = psmeangfo(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)
                   obsmeangfo(i,j,bi,bj) = obsmeangfo(i,j,bi,bj) +
     &                 gfoobs(i,j,bi,bj)
                   sumgfo(i,j,bi,bj) = sumgfo(i,j,bi,bj) + 1. _d 0
                endif
#endif
              enddo
            enddo
          enddo
        enddo

c--   END loop over records for the first time.
      enddo

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
CMM(  average up all SLA
               mean_slaobs(i,j,bi,bj) = obsmeantp(i,j,bi,bj) +
     &               obsmeaners(i,j,bi,bj) + obsmeangfo(i,j,bi,bj) 
               mean_slaobs_NUM(i,j,bi,bj) = sumtp(i,j,bi,bj) +
     &               sumers(i,j,bi,bj) + sumgfo(i,j,bi,bj)
               if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then
                  mean_slaobs(i,j,bi,bj) = mean_slaobs(i,j,bi,bj) /
     &                 mean_slaobs_NUM(i,j,bi,bj)
               else
                  mean_slaobs(i,j,bi,bj) = 0. _d 0
               endif
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
               if ( sumtp(i,j,bi,bj) .NE. 0. ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 psmeantp(i,j,bi,bj) - obsmeantp(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 sumtp(i,j,bi,bj)
                  psmeantp(i,j,bi,bj) = psmeantp(i,j,bi,bj) /
     &                 sumtp(i,j,bi,bj)
                  obsmeantp(i,j,bi,bj) = obsmeantp(i,j,bi,bj) /
     &                 sumtp(i,j,bi,bj)
                  wwwtp1(i,j,bi,bj) = 1. _d 0
               endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
               if ( sumers(i,j,bi,bj) .NE. 0. ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 psmeaners(i,j,bi,bj) - obsmeaners(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 sumers(i,j,bi,bj)
                  psmeaners(i,j,bi,bj) = psmeaners(i,j,bi,bj) /
     &                 sumers(i,j,bi,bj)
                  obsmeaners(i,j,bi,bj) = obsmeaners(i,j,bi,bj) /
     &                 sumers(i,j,bi,bj)
                  wwwers1(i,j,bi,bj) = 1. _d 0
               endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
               if ( sumgfo(i,j,bi,bj) .NE. 0. ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 psmeangfo(i,j,bi,bj)-obsmeangfo(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 sumgfo(i,j,bi,bj)
                  psmeangfo(i,j,bi,bj) = psmeangfo(i,j,bi,bj) /
     &                 sumgfo(i,j,bi,bj)
                  obsmeangfo(i,j,bi,bj) = obsmeangfo(i,j,bi,bj) /
     &                 sumgfo(i,j,bi,bj)
                  wwwgfo1(i,j,bi,bj) = 1. _d 0
               endif
#endif
               if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( tpmeanmask(i,j,bi,bj) .NE. 0. ) ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) /
     &                 mean_psMssh_all_NUM(i,j,bi,bj)-tpmean(i,j,bi,bj)
     &                 -mean_slaobs(i,j,bi,bj)
                  mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                  offset=offset+mean_psMssh_all(i,j,bi,bj)*
     &                 mean_psMssh_all_NUM(i,j,bi,bj)
                  offset_sum=offset_sum+mean_psMssh_all_NUM(i,j,bi,bj)
               else
                  mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                  mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
               endif
CMM(  OFFSET IS < DOT + obs - model > 
              enddo
            enddo
          enddo
        enddo
CMM( output SLA mean for referencing
      write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
      call mdswritefield(fname4test,32,.false.,'RL',
     & 1,mean_slaobs,1,1,mythid)
CMM)

c--   Do a global summation.
      _GLOBAL_SUM_RL( offset     , mythid )
      _GLOBAL_SUM_RL( offset_sum , mythid )

      if (offset_sum .eq. 0.0) then
        _BEGIN_MASTER( mythid )
        write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'
        call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
        stop   '  ... stopped in cost_ssh.'
      else
        _BEGIN_MASTER( mythid )
        write(msgbuf,'(a,d22.15)')
     &        ' cost_ssh: offset_sum = ',offset_sum
        call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        write(msgbuf,'(a,d22.15)')
     &          ' cost_ssh: offset_tot = ',offset
        call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
      endif

      offset = offset / offset_sum

#ifdef ALLOW_PROFILES
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = 1,sny
            do i = 1,snx
              prof_etan_mean(i,j,bi,bj)=offset+psmean(i,j,bi,bj)
            enddo
          enddo
        enddo
      enddo
      _EXCH_XY_RL( prof_etan_mean, mythid )
#endif

#ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
#ifndef ALLOW_SSH_TOT
c--   ==========
c--      Mean
c--   ==========
c--   compute mean ssh difference cost contribution
c--      Mean
CMM   mean_psMssh_all  = <model - obs> - (DOT + <obs>)
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
       if ( psmean(i,j,bi,bj) .ne. 0. ) then
          diagnosfld(i,j,bi,bj) = (mean_psMssh_all(i,j,bi,bj) -  offset)
     &     *mean_psMssh_all_MSK(i,j,bi,bj)
         diff = mean_psMssh_all(i,j,bi,bj) -  offset
         sumc = sumc + diff*diff
     &          * wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj)
         if ( wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) .ne. 0. )
     &             num_hmean_mm = num_hmean_mm + 1. _d 0
CMM         diagnosfld(i,j,bi,bj) = diff*diff 
CMM     &      * wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj)
         endif
            enddo
          enddo
        enddo
      enddo

      CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean1', ' ', diagnosfld,
     &                           optimcycle, mythid )

      objf_hmean = objf_hmean + sumc
      num_hmean = num_hmean + num_hmean_mm
CMM)
CMM(((((HERE ADD OTHER DOT FOR DOT PROJECT
CMM     THIS IS ALL HARDCODED!
CMM     DO THIS for 4 DOTS -- first read the field
CMM first get orig cost
c--   Do a global summation.
       _GLOBAL_SUM_RL( sumc         , mythid )
       _GLOBAL_SUM_RL( num_hmean_mm , mythid )
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT tmp ref: offset ', offset
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT tmp ref: offset N ', offset_sum
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT tmp ref: cost integral ', sumc
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT tmp ref: cost N ', num_hmean_mm
CMM  now to 
       factor =  0.01
       spval  = -9990.
       write(standardmessageunit,'(A)')
     & 'CMM ssh_cost - reading DOT in order: EGM08 AVISO GRACE MAXI'
       do IDOT = 1,4
CMM reinitialize due to recomps
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
                tpmean_loc(i,j,bi,bj) = 0. _d 0
              enddo
            enddo
          enddo
        enddo
        sumc = 0.
        num_hmean_mm = 0.
        if     (IDOT.eq.1) then
        call mdsreadfield( 'DOT_egm08.bin', cost_iprec,cost_yftype
     &                  ,1,tpmean_loc, 1, mythid )
        elseif (IDOT.eq.2) then
        call mdsreadfield( 'DOT_aviso.bin', cost_iprec,cost_yftype,
     &                   1,tpmean_loc, 1, mythid )
        elseif (IDOT.eq.3) then
        call mdsreadfield( 'DOT_grace2.bin', cost_iprec,cost_yftype,
     &                   1,tpmean_loc, 1, mythid )
        elseif (IDOT.eq.4) then
        call mdsreadfield( 'DOT_maxi.bin', cost_iprec,cost_yftype,
     &                   1,tpmean_loc, 1, mythid )
        endif
CMM  NOW MAKE MASK
        do bj = jtlo,jthi
        do bi = itlo,ithi
          k = 1
          do j = jmin,jmax
            do i = imin,imax
              mean_psMssh_all_MSK(i,j,bi,bj) =  tpmeanmask(i,j,bi,bj)
              if (tpmean_loc(i,j,bi,bj) .lt. spval) then
                mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
              endif
              if (tpmean_loc(i,j,bi,bj) .eq. 0. _d 0 ) then
                mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
              endif
              tpmean_loc(i,j,bi,bj) = tpmean_loc(i,j,bi,bj)*
     &                            mean_psMssh_all_MSK(i,j,bi,bj)*
     &                             factor
            enddo
          enddo
        enddo
        enddo
CMM NOW GET OFFSET
c--   Compute and remove offset for current tile and sum over all
c--   tiles of this instance.
      offset     = 0. _d 0
      offset_sum = 0. _d 0
c--   Sum over this thread tiles.
        do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = 1,sny
            do i = 1,snx
              offset     = offset +
     &                mean_psMssh_all_MSK(i,j,bi,bj)*cosphi(i,j,bi,bj)*
     &                     (tpmean_loc(i,j,bi,bj) - psmean(i,j,bi,bj))
              offset_sum = offset_sum +
     &                mean_psMssh_all_MSK(i,j,bi,bj)*cosphi(i,j,bi,bj)
            enddo
          enddo
        enddo
        enddo

c--   Do a global summation.
        _GLOBAL_SUM_RL( offset     , mythid )
        _GLOBAL_SUM_RL( offset_sum , mythid )

        if (offset_sum .eq. 0.0) then
        _BEGIN_MASTER( mythid )
          write(msgbuf,'(a)') ' cost_ssh: offset_sum = zero!'
          call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
          stop   '  ... stopped in cost_ssh.'
        else
        _BEGIN_MASTER( mythid )
          write(msgbuf,'(a,d22.15)')
     &          ' cost_ssh: offset = ',offset
          call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
          write(msgbuf,'(a,d22.15)')
     &          ' cost_ssh: offset_sum = ',offset_sum
          call print_message( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)

        _END_MASTER( mythid )
        endif
        offset = offset / offset_sum
        do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
          diagnosfld(i,j,bi,bj)=(psmean(i,j,bi,bj)-tpmean_loc(i,j,bi,bj)
     &              + offset)*mean_psMssh_all_MSK(i,j,bi,bj)
              diff = psmean(i,j,bi,bj) - tpmean_loc(i,j,bi,bj) + offset
              sumc = sumc + diff*diff
     &           *wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj)
             if ( wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj) .ne. 0. )
     &             num_hmean_mm = num_hmean_mm + 1. _d 0
CMM              diagnosfld(i,j,bi,bj) = diff*diff
CMM     &           *wp(i,j,bi,bj)*mean_psMssh_all_MSK(i,j,bi,bj)
            enddo
          enddo
        enddo
        enddo

        if     (IDOT.eq.1) then
          CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_A', ' ', diagnosfld,
     &                           optimcycle, mythid )
        elseif (IDOT.eq.2) then
          CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_B', ' ', diagnosfld,
     &                           optimcycle, mythid )
        elseif (IDOT.eq.3) then
          CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_C', ' ', diagnosfld,
     &                           optimcycle, mythid )
        elseif (IDOT.eq.4) then
          CALL WRITE_FLD_XY_RL( 'DiagnosSSHmean_mm_D', ' ', diagnosfld,
     &                           optimcycle, mythid )
        endif

CMM ADD TO COST FUNCTION
      objf_hmean = objf_hmean + sumc
      num_hmean = num_hmean + num_hmean_mm

c--   Do a global summation.
      _GLOBAL_SUM_RL( sumc     , mythid )
      _GLOBAL_SUM_RL( num_hmean_mm , mythid )

       write(standardmessageunit,'(A,I6)')
     & 'CMM DOT sta ref number: ', IDOT
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT sta ref: offset ', offset
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT sta ref: offset N ', offset_sum  
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT sta ref: cost integral ', sumc
       write(standardmessageunit,'(A,D22.15)')
     & 'CMM DOT sta ref: cost N ', num_hmean_mm

      enddo
CMM end of loop over DOT products


#endif /* ALLOW_SSH_TOT */
#endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */

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

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

        call active_read_xy( fname, psbar, irec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_psbar_mean_dummy )

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
        call cost_readtopex( irec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
        call cost_readers( irec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
        call cost_readgfo( irec, mythid )
#endif

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

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
            do j = jmin,jmax
              do i = imin,imax
c--             The array psobs contains SSH anomalies.
                wwwtp2(i,j,bi,bj) = wwwtp1(i,j,bi,bj)
     &                             *wtp(i,j,bi,bj)
     &                             *tpmask(i,j,bi,bj)
     &                             *cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
CMM                junk = ( psbar(i,j,bi,bj)
CMM     &                   - psmeantp(i,j,bi,bj) - tpobs(i,j,bi,bj) )
                junk =  ( psbar(i,j,bi,bj) - psmeantp(i,j,bi,bj)  )
     &                - ( tpobs(i,j,bi,bj) - obsmeantp(i,j,bi,bj) )
                objf_tp(bi,bj) = objf_tp(bi,bj)
     &              + junk*junk*wwwtp2(i,j,bi,bj)
                if ( wwwtp2(i,j,bi,bj) .ne. 0. )
     &               num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0
#else
                if (tpmeanmask(i,j,bi,bj)*tpmask(i,j,bi,bj)
     &                  *wp(i,j,bi,bj)*wwwtp2(i,j,bi,bj) .ne.0.) then
                junk       = ( psbar(i,j,bi,bj) -
     &                 (tpobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) )
                objf_tp(bi,bj) = objf_tp(bi,bj)
     &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwtp2(i,j,bi,bj) )
                num_tp(bi,bj) = num_tp(bi,bj) + 1. _d 0
                endif
#endif /* ALLOW_SSH_TOT */
              enddo
            enddo
#endif

#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
            do j = jmin,jmax
              do i = imin,imax
c--             The array ersobs contains SSH anomalies.
                wwwers2(i,j,bi,bj) = wwwers1(i,j,bi,bj)
     &                              *wers(i,j,bi,bj)
     &                              *ersmask(i,j,bi,bj)
     &                              *cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
CMM                junk = ( psbar(i,j,bi,bj) 
CMM     &                   - psmeaners(i,j,bi,bj) - ersobs(i,j,bi,bj) )
                junk =  ( psbar(i,j,bi,bj) - psmeaners(i,j,bi,bj)  )
     &                - ( ersobs(i,j,bi,bj) - obsmeaners(i,j,bi,bj) )
                objf_ers(bi,bj) = objf_ers(bi,bj)
     &              + junk*junk*wwwers2(i,j,bi,bj)
                if ( wwwers2(i,j,bi,bj) .ne. 0. )
     &               num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0
#else
                if (tpmeanmask(i,j,bi,bj)*ersmask(i,j,bi,bj)
     &                  *wp(i,j,bi,bj)*wwwers2(i,j,bi,bj) .ne.0.) then
                junk       = ( psbar(i,j,bi,bj) -
     &                 (ersobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) )
                objf_ers(bi,bj) = objf_ers(bi,bj)
     &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwers2(i,j,bi,bj) )
                num_ers(bi,bj) = num_ers(bi,bj) + 1. _d 0
                endif
#endif /* ALLOW_SSH_TOT */
              enddo
            enddo
#endif

#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
            do j = jmin,jmax
              do i = imin,imax
c--             The array gfoobs contains SSH anomalies.
                wwwgfo2(i,j,bi,bj) = wwwgfo1(i,j,bi,bj)
     &                        *wgfo(i,j,bi,bj)
     &                        *gfomask(i,j,bi,bj)
     &                        *cosphi(i,j,bi,bj)
#ifndef ALLOW_SSH_TOT
CMM                junk = ( psbar(i,j,bi,bj)
CMM     &                   - psmeangfo(i,j,bi,bj) - gfoobs(i,j,bi,bj) )
                junk =  ( psbar(i,j,bi,bj) - psmeangfo(i,j,bi,bj)  )
     &                - ( gfoobs(i,j,bi,bj) - obsmeangfo(i,j,bi,bj) )
                objf_gfo(bi,bj) = objf_gfo(bi,bj)
     &              + junk*junk*wwwgfo2(i,j,bi,bj)
                if ( wwwgfo2(i,j,bi,bj) .ne. 0. )
     &               num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0
#else
                if (tpmeanmask(i,j,bi,bj)*gfomask(i,j,bi,bj)
     &                  *wp(i,j,bi,bj)*wwwgfo2(i,j,bi,bj) .ne.0.) then
                junk       = ( psbar(i,j,bi,bj) -
     &                 (gfoobs(i,j,bi,bj)+tpmean(i,j,bi,bj)-offset) )
                objf_gfo(bi,bj) = objf_gfo(bi,bj)
     &              +junk*junk/( 1/wp(i,j,bi,bj)+1/wwwgfo2(i,j,bi,bj) )
                num_gfo(bi,bj) = num_gfo(bi,bj) + 1. _d 0
                endif
#endif /* ALLOW_SSH_TOT */
              enddo
            enddo
#endif

          enddo
        enddo

      enddo
c--   End of second loop over records.

      do bj = jtlo,jthi
        do bi = itlo,ithi
          objf_h(bi,bj) = objf_h(bi,bj) +
     &        objf_tp(bi,bj) + objf_ers(bi,bj) + objf_gfo(bi,bj)
          num_h(bi,bj) = num_h(bi,bj) +
     &        num_tp(bi,bj) + num_ers(bi,bj) + num_gfo(bi,bj)
        enddo
      enddo


#endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */

      end