C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_sshv4.F,v 1.23 2014/06/09 17:55:22 gforget Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"


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

c     ==================================================================
c     SUBROUTINE cost_gencost_sshv4
c     ==================================================================
c
c     o Evaluate cost function contributions of sea surface height.
c
c        started: Gael Forget, Oct-2009
c
c        working assumption for the time mean dynamic topography (MDT) constraint:
c        the various SLA data sets (tp, ers, gfo) have been consistenty cross-referenced,
c        as done in the RADS data sets. We do not need to know the reference dynamic
c        topography (or SSH/Geoid). Simply it is assumed that the biases
c        between instruments have been taken care of. This is only a matter
c        for the MDT constraint, not for SLA constraints (see below).
c
cgf 1) there are a few hardcoded numbers that will eventually be put in common
cgf     blocks/namelists
cgf 2) there are a several refinements that should be considered, such as
cgf     modulating weights with respect to numbers of samples
c
c     ==================================================================
c     SUBROUTINE cost_gencost_sshv4
c     ==================================================================

      implicit none

c     == global variables ==

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

#include "ecco_cost.h"
#include "CTRL_SIZE.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"
#include "DYNVARS.h"
#ifdef ALLOW_PROFILES
#include "profiles.h"
#endif
#include "cal.h"
#ifdef ALLOW_SMOOTH
#include "SMOOTH.h"
#endif

c     == routine arguments ==

      integer myiter
      _RL     mytime
      integer mythid

#ifdef ALLOW_SSH_COST_CONTRIBUTION
#ifdef ALLOW_GENCOST_CONTRIBUTION

c     == functions ==
      LOGICAL  MASTER_CPU_THREAD
      EXTERNAL 

c     == local variables ==

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

      logical doglobalread
      logical ladinit

c mapping to gencost
      integer igen_mdt, igen_lsc, igen_gmsl
      integer igen_tp, igen_ers, igen_gfo

      _RL offset
      _RL offset_sum
      _RL slaoffset
      _RL slaoffset_sum
      _RL slaoffset_weight

      _RL diagnosfld ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )

c for PART 1: re-reference MDT (mdt) to the inferred SLA reference field
      _RL psmean    ( 1-olx:snx+olx, 1-oly:sny+oly, nsx, nsy )
      _RL mean_slaobs_mdt(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_slaobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

c for PART 2: compute time mean differences over the model period
      _RL mean_slaobs_model(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 mean_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMslaobs_MSK(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL mean_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMtpobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL mean_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMersobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL mean_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL mean_psMgfoobs_NUM(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

c for PART 4/5: compute smooth/raw anomalies
      _RL anom_psMslaobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_psMslaobs_NUM (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_psMslaobs_MSK (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_psMslaobs_SUB (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_slaobs (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_slaobs_SUB  (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL anom_psMtpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_tpobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL anom_psMersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_ersobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      _RL anom_psMgfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      _RL anom_gfoobs(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)

      integer mdt_y0,mdt_y1,year,day
      integer num_var, num0

      _RL junk,junkweight

      integer ndaysave
      _RL ndaysaveRL

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

      LOGICAL doReference

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-- detect the relevant gencost indices
      igen_mdt=0
      igen_tp =0
      igen_ers=0
      igen_gfo=0
      igen_lsc=0
      igen_gmsl=0
      do k=1,NGENCOST
        if (gencost_name(k).EQ.'sshv4-mdt') igen_mdt=k
        if (gencost_name(k).EQ.'sshv4-tp') igen_tp=k
        if (gencost_name(k).EQ.'sshv4-ers') igen_ers=k
        if (gencost_name(k).EQ.'sshv4-gfo') igen_gfo=k
        if (gencost_name(k).EQ.'sshv4-lsc') igen_lsc=k
        if (gencost_name(k).EQ.'sshv4-gmsl') igen_gmsl=k
      enddo

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


cgf =======================================================
cgf PART 1:
cgf        x Get the MDT (mdt) ... compute the sample mean
cgf        (mean_slaobs_mdt) of the SLA data (i.e. RADS for tp, ers, and gfo
cgf        together) over the time interval of the MDT ... subtract
cgf        mean_slaobs_mdt from mdt.
cgf        x At this point, mdt is the inferred SLA reference field.
cgf        x From there on, mdt+sla will be directly comparable to
cgf        the model SSH (psbar).
cgf =======================================================

c--   Read mean field and mask
      call COST_READTOPEXMEAN( mythid )

c--   Compute mean_slaobs_mdt: sample mean SLA over the time period of mdt.

c pavlis and ecco/rio
c      mdt_y0=1993
c      mdt_y1=2004
c maximenko
c      mdt_y0=1992
c      mdt_y1=2002
c rio
c      mdt_y0=1993
c      mdt_y1=1999
c generalized:
       mdt_y0=mdtstartdate(1)/10000
       mdt_y1=mdtenddate(1)/10000

       write(msgbuf,'(a,i8,a,i8)')
     &  'mdt[start,end]date(1): ',mdtstartdate(1),
     &  ',',mdtenddate(1)
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)

        write(msgbuf,'(a,i4,a,i4)')
     &  '  TP MDT yrrange:          ',
     &                        mdt_y0,',',
     &                        mdt_y1
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                      SQUEEZE_RIGHT , mythid)

c minimum number of observations to consider re-referenced MDT
      num0=100

      doReference=.FALSE.
      if ((modelstartdate(1).GT.1992*10000).AND.
     &    (modelstartdate(1).LT.2011*10000).AND.
     &    (ndaysrec.GE.365))  doReference=.TRUE.

        write(msgbuf,'(a,l)') ' sshv4:re-reference MDT=',doReference
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
              mean_slaobs_mdt(i,j,bi,bj)  = 0. _d 0
              mean_slaobs_NUM(i,j,bi,bj)  = 0. _d 0
          enddo
         enddo
        enddo
       enddo

      do year=mdt_y0,mdt_y1
       do day=1,366
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      call COST_SLA_READ_YD( topexfile, topexstartdate,
     &                tpobs, tpmask,
     &                year, day, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      call COST_SLA_READ_YD( ersfile, ersstartdate,
     &                ersobs, ersmask,
     &                year, day, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      call COST_SLA_READ_YD( gfofile, gfostartdate,
     &                gfoobs, gfomask,
     &                year, day, mythid )
#endif
       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      if ( tpmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
     &    gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
          mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
     &  tpobs(i,j,bi,bj)
          mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
      endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      if ( ersmask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
     &    gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
          mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
     &  ersobs(i,j,bi,bj)
          mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
      endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      if ( gfomask(i,j,bi,bj)*mdtmask(i,j,bi,bj)*
     &    gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
          mean_slaobs_mdt(i,j,bi,bj)= mean_slaobs_mdt(i,j,bi,bj)+
     &  gfoobs(i,j,bi,bj)
          mean_slaobs_NUM(i,j,bi,bj)= mean_slaobs_NUM(i,j,bi,bj)+1. _d 0
      endif
#endif
          enddo
         enddo
        enddo
       enddo

       enddo !do day=1,366
      enddo !do year=mdt_y0,mdt_y1

      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
               if ( ( mean_slaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
     &              ( mdtmask(i,j,bi,bj) .NE. 0. ) ) then
                  mean_slaobs_mdt(i,j,bi,bj) =
     &                 mean_slaobs_mdt(i,j,bi,bj) /
     &                 mean_slaobs_NUM(i,j,bi,bj)
               else
                  mean_slaobs_mdt(i,j,bi,bj) = 0. _d 0
               endif
            enddo
          enddo
        enddo
      enddo


c--   smooth mean_slaobs_mdt:

      write(fname4test(1:80),'(1a)') 'sla2mdt_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,mean_slaobs_mdt,1,1,mythid)

#ifdef ALLOW_SMOOTH
      if ( useSMOOTH )
     &  call SMOOTH_HETERO2D(mean_slaobs_mdt,maskc,
     &     gencost_scalefile(igen_mdt),
     &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
#endif

      write(fname4test(1:80),'(1a)') 'sla2mdt_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,mean_slaobs_mdt,1,1,mythid)

c--   re-reference mdt:
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
               if ( ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
     &              ( doReference ) ) then
                  mdt(i,j,bi,bj) = mdt(i,j,bi,bj)
     &                 -mean_slaobs_mdt(i,j,bi,bj)
               endif
            enddo
          enddo
        enddo
      enddo


cgf =======================================================
cgf PART 2: compute sample means of psbar-slaobs over the
cgf          period that is covered by the model (i.e. psbar).
cgf          x for all SLA data sets together: mean_psMssh_all, mean_psMssh_all_MSK,
cgf             and offset will be used in PART 3 (MDT cost term).
cgf          x for each SLA data individually. mean_psMtpobs, mean_psMtpobs_MS, etc.
cgf             will be used in PART 4&5 (SLA cost terms).
cgf =======================================================

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

              psmean(i,j,bi,bj)    = 0. _d 0
              mean_psMslaobs(i,j,bi,bj)  = 0. _d 0
              mean_psMtpobs(i,j,bi,bj)  = 0. _d 0
              mean_psMersobs(i,j,bi,bj) = 0. _d 0
              mean_psMgfoobs(i,j,bi,bj) = 0. _d 0
              mean_psMssh_all(i,j,bi,bj) = 0. _d 0
              mean_slaobs_model(i,j,bi,bj)  = 0. _d 0

              mean_psMtpobs_NUM(i,j,bi,bj)  = 0. _d 0
              mean_psMersobs_NUM(i,j,bi,bj) = 0. _d 0
              mean_psMgfoobs_NUM(i,j,bi,bj) = 0. _d 0
              mean_psMssh_all_NUM(i,j,bi,bj) = 0. _d 0
              mean_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0

            enddo
          enddo
        enddo
      enddo
      offset     = 0. _d 0
      offset_sum = 0. _d 0


      do irec = 1, ndaysrec

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

#ifndef ALLOW_PSBAR_MEAN
        CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
     &        'psbar', myTime, myThid )
#endif

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      call COST_SLA_READ( topexfile, topexstartdate, topexperiod,
     &                topexintercept, topexslope,
     &                tpobs, tpmask,
     &                irec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      call COST_SLA_READ( ersfile, ersstartdate, ersperiod,
     &                ersintercept, ersslope,
     &                ersobs, ersmask,
     &                irec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      call COST_SLA_READ( gfofile, gfostartdate, gfoperiod,
     &                gfointercept, gfoslope,
     &                gfoobs, gfomask,
     &                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)*
     &             gencost_weight(i,j,bi,bj,igen_tp) .NE. 0. ) then
                   mean_slaobs_model(i,j,bi,bj)=
     &                 mean_slaobs_model(i,j,bi,bj)+tpobs(i,j,bi,bj)
                   mean_psMtpobs(i,j,bi,bj) =
     &                 mean_psMtpobs(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
                   mean_psMtpobs_NUM(i,j,bi,bj) =
     &                 mean_psMtpobs_NUM(i,j,bi,bj) + 1. _d 0
                endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
                if ( ersmask(i,j,bi,bj)*
     &             gencost_weight(i,j,bi,bj,igen_ers) .NE. 0. ) then
                   mean_slaobs_model(i,j,bi,bj)=
     &                 mean_slaobs_model(i,j,bi,bj)+ersobs(i,j,bi,bj)
                   mean_psMersobs(i,j,bi,bj) =
     &                 mean_psMersobs(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
                   mean_psMersobs_NUM(i,j,bi,bj) =
     &                 mean_psMersobs_NUM(i,j,bi,bj) + 1. _d 0
                endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
                if ( gfomask(i,j,bi,bj)*
     &             gencost_weight(i,j,bi,bj,igen_gfo) .NE. 0. ) then
                   mean_slaobs_model(i,j,bi,bj)=
     &                 mean_slaobs_model(i,j,bi,bj)+gfoobs(i,j,bi,bj)
                   mean_psMgfoobs(i,j,bi,bj) =
     &                 mean_psMgfoobs(i,j,bi,bj) +
     &                 psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
                   mean_psMgfoobs_NUM(i,j,bi,bj) =
     &                 mean_psMgfoobs_NUM(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
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
               if ( ( mean_psMtpobs_NUM(i,j,bi,bj) .NE. 0. )
     &              ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 mean_psMtpobs(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 mean_psMtpobs_NUM(i,j,bi,bj)
                  mean_psMtpobs(i,j,bi,bj) =
     &                 mean_psMtpobs(i,j,bi,bj) /
     &                 mean_psMtpobs_NUM(i,j,bi,bj)
               endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
               if ( ( mean_psMersobs_NUM(i,j,bi,bj) .NE. 0. )
     &              ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 mean_psMersobs(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 mean_psMersobs_NUM(i,j,bi,bj)
                  mean_psMersobs(i,j,bi,bj) =
     &                 mean_psMersobs(i,j,bi,bj) /
     &                 mean_psMersobs_NUM(i,j,bi,bj)
               endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
               if ( ( mean_psMgfoobs_NUM(i,j,bi,bj) .NE. 0. )
     &              ) then
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) +
     &                 mean_psMgfoobs(i,j,bi,bj)
                  mean_psMssh_all_NUM(i,j,bi,bj) =
     &                 mean_psMssh_all_NUM(i,j,bi,bj) +
     &                 mean_psMgfoobs_NUM(i,j,bi,bj)
                  mean_psMgfoobs(i,j,bi,bj) =
     &                 mean_psMgfoobs(i,j,bi,bj) /
     &                 mean_psMgfoobs_NUM(i,j,bi,bj)
               endif
#endif

               if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. 0 ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. )
     &              ) then
                  mean_psMslaobs(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) /
     &                 mean_psMssh_all_NUM(i,j,bi,bj)
                  mean_psMslaobs_MSK(i,j,bi,bj) = 1. _d 0
               else
                  mean_psMslaobs(i,j,bi,bj) = 0. _d 0
                  mean_psMslaobs_MSK(i,j,bi,bj) = 0. _d 0
               endif

               if ( ( mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0 ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
     &              ( mdtmask(i,j,bi,bj) .NE. 0. ).AND.
     &              ( doReference ) ) then
                  mean_slaobs_model(i,j,bi,bj) =
     &                 mean_slaobs_model(i,j,bi,bj) /
     &                 mean_psMssh_all_NUM(i,j,bi,bj)
                  mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) /
     &                 mean_psMssh_all_NUM(i,j,bi,bj)-mdt(i,j,bi,bj)
                  mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                  offset=offset+RA(i,j,bi,bj)*mean_psMssh_all(i,j,bi,bj)
                  offset_sum=offset_sum+RA(i,j,bi,bj)
               elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ).AND.
     &              ( .NOT.doReference ) ) then
                  mean_slaobs_model(i,j,bi,bj) = 0.d0
                  mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                  mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0
                  offset=offset+RA(i,j,bi,bj)
     &                  *( psmean(i,j,bi,bj) -mdt(i,j,bi,bj) )
                  offset_sum=offset_sum+RA(i,j,bi,bj)
               else
                  mean_slaobs_model(i,j,bi,bj) = 0.d0
                  mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                  mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
               endif
              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)') ' sshv4: offset_sum = zero!'
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
c        stop   '  ... stopped in cost_ssh.'
      else
        _BEGIN_MASTER( mythid )
        write(msgbuf,'(a,2d22.15)') ' sshv4:offset,sum=',
     &      offset,offset_sum
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
      endif

c--   Compute (average) offset
      offset = offset / offset_sum

c--   subtract offset from mean_psMssh_all and psmean
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax

               if ( (mean_psMssh_all_NUM(i,j,bi,bj) .GT. num0) .AND.
     &              ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND. 
     &              ( doReference ) ) then
c use the re-referencing approach
                   mean_psMssh_all(i,j,bi,bj) =
     &                 mean_psMssh_all(i,j,bi,bj) - offset
                   mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0

               elseif ( ( mdtmask(i,j,bi,bj) .NE. 0. ) .AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) .AND.
     &              ( .NOT.doReference ) ) then
c use the simpler approach
                   mean_psMssh_all(i,j,bi,bj) =
     &             psmean(i,j,bi,bj) -mdt(i,j,bi,bj) - offset
                   mean_psMssh_all_MSK(i,j,bi,bj) = 1. _d 0

               else
                   mean_psMssh_all(i,j,bi,bj) = 0. _d 0
                   mean_psMssh_all_MSK(i,j,bi,bj) = 0. _d 0
               endif

               if ( maskc(i,j,1,bi,bj) .NE. 0. )
     &            psmean(i,j,bi,bj)=psmean(i,j,bi,bj)-offset

            enddo
          enddo
        enddo
       enddo

c--    smooth mean_psMssh_all
      write(fname4test(1:80),'(1a)') 'mdtdiff_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     &    1,mean_psMssh_all,1,1,mythid)

#ifdef ALLOW_SMOOTH
      if ( useSMOOTH )
     &  call SMOOTH_HETERO2D(mean_psMssh_all,maskc,
     &     gencost_scalefile(igen_mdt),
     &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
#endif

      write(fname4test(1:80),'(1a)') 'mdtdiff_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     &    1,mean_psMssh_all,1,1,mythid)

      write(fname4test(1:80),'(1a)') 'sla2model_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     &    1,mean_slaobs_model,1,1,mythid)

#ifdef ALLOW_SMOOTH
      if ( useSMOOTH )
     &  call SMOOTH_HETERO2D(mean_slaobs_model,maskc,
     &     gencost_scalefile(igen_mdt),
     &     gencost_smooth2Ddiffnbt(igen_mdt),mythid)
#endif

      write(fname4test(1:80),'(1a)') 'sla2model_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     &    1,mean_slaobs_model,1,1,mythid)

cgf at this point:
cgf     1) mean_psMssh_all is the sample mean ,
cgf             to which a smoothing filter has been applied.
cgf     2) mean_psMtpobs is the (unsmoothed) sample mean .
cgf             And similarly for ers and gfo, each treated separately.

#ifdef ALLOW_PROFILES
      if ( usePROFILES ) then
      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = 1,sny
            do i = 1,snx
cgf this has not been tested recently
              prof_etan_mean(i,j,bi,bj)=psmean(i,j,bi,bj)
            enddo
          enddo
        enddo
      enddo
      _EXCH_XY_RL( prof_etan_mean, mythid )
      endif
#endif


cgf =======================================================
cgf PART 3: compute MDT cost term
cgf =======================================================


#ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION

      do bj = jtlo,jthi
        do bi = itlo,ithi
          do j = jmin,jmax
            do i = imin,imax
       if (mean_psMssh_all_MSK(i,j,bi,bj).NE.0. _d 0) then
         junk = mean_psMssh_all(i,j,bi,bj)
         junkweight = gencost_weight(i,j,bi,bj,igen_mdt)
     &              * mdtmask(i,j,bi,bj)
       else
         junk = 0. _d 0
         junkweight = 0. _d 0
       endif
       objf_gencost(bi,bj,igen_mdt) = objf_gencost(bi,bj,igen_mdt)
     &      + junk*junk*junkweight
       if ( junkweight .ne. 0. ) num_gencost(bi,bj,igen_mdt) =
     &      num_gencost(bi,bj,igen_mdt) + 1. _d 0
       diagnosfld(i,j,bi,bj) = junk*junk*junkweight
            enddo
          enddo
        enddo
      enddo

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

#endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */



cgf =======================================================
cgf PART 4: compute smooth SLA cost term
cgf =======================================================


      ndaysave=35
      ndaysaveRL=ndaysave

      do irec = 1, ndaysrec-ndaysave+1

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
              anom_psMslaobs(i,j,bi,bj)  = 0. _d 0
              anom_psMslaobs_MSK(i,j,bi,bj)  = 0. _d 0
              anom_psMslaobs_SUB(i,j,bi,bj)  = 0. _d 0
              anom_slaobs(i,j,bi,bj)  = 0. _d 0
              anom_slaobs_SUB(i,j,bi,bj)  = 0. _d 0
              anom_psMslaobs_NUM(i,j,bi,bj)  = 0. _d 0
          enddo
         enddo
        enddo
       enddo

c PART 4.1: compute running sample average over ndaysave
c ------------------------------------------------------

      do jrec=1,ndaysave

        krec=irec+jrec-1

        call ACTIVE_READ_XY( fname, psbar, krec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_psbar_mean_dummy )

#ifndef ALLOW_PSBAR_MEAN
        CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
     &        'psbar', myTime, myThid )
#endif

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      call COST_SLA_READ( topexfile, topexstartdate, topexperiod,
     &                topexintercept, topexslope,
     &                tpobs, tpmask,
     &                krec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      call COST_SLA_READ( ersfile, ersstartdate, ersperiod,
     &                ersintercept, ersslope,
     &                ersobs, ersmask,
     &                krec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      call COST_SLA_READ( gfofile, gfostartdate, gfoperiod,
     &                gfointercept, gfoslope,
     &                gfoobs, gfomask,
     &                krec, mythid )
#endif

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      if ( tpmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
     &  .NE.0. ) then
           anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
     &        psbar(i,j,bi,bj)-tpobs(i,j,bi,bj)
     &        -mean_psMslaobs(i,j,bi,bj)
           anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
     &        tpobs(i,j,bi,bj)
           anom_psMslaobs_NUM(i,j,bi,bj)=
     &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
           if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) ) 
     &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
      endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      if ( ersmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
     &  .NE.0. ) then
           anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
     &        psbar(i,j,bi,bj)-ersobs(i,j,bi,bj)
     &        -mean_psMslaobs(i,j,bi,bj)
           anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
     &        ersobs(i,j,bi,bj)
           anom_psMslaobs_NUM(i,j,bi,bj)=
     &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
           if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
     &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
      endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      if ( gfomask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj)
     &  .NE.0. ) then
           anom_psMslaobs(i,j,bi,bj)= anom_psMslaobs(i,j,bi,bj)+
     &        psbar(i,j,bi,bj)-gfoobs(i,j,bi,bj)
     &        -mean_psMslaobs(i,j,bi,bj)
           anom_slaobs(i,j,bi,bj)= anom_slaobs(i,j,bi,bj)+
     &        gfoobs(i,j,bi,bj)
           anom_psMslaobs_NUM(i,j,bi,bj)=
     &        anom_psMslaobs_NUM(i,j,bi,bj)+1. _d 0
           if ( (2*jrec.EQ.ndaysave).OR.(2*jrec+1.EQ.ndaysave) )
     &        anom_psMslaobs_MSK(i,j,bi,bj)  = 1. _d 0
      endif
#endif
          enddo
         enddo
        enddo
       enddo

      enddo !do jrec=1,ndaysave

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
               if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) 
     &            ) then
                  anom_psMslaobs(i,j,bi,bj) =
     &                 anom_psMslaobs(i,j,bi,bj) /
     &                 anom_psMslaobs_NUM(i,j,bi,bj)
                  anom_slaobs(i,j,bi,bj) =
     &                 anom_slaobs(i,j,bi,bj) /
     &                 anom_psMslaobs_NUM(i,j,bi,bj)
               else
                  anom_psMslaobs(i,j,bi,bj) = 0. _d 0
                  anom_slaobs(i,j,bi,bj) = 0. _d 0
               endif
              enddo
            enddo
          enddo
        enddo

c PART 4.11: separate time variable offset
c ----------------------------------------

      slaoffset     = 0. _d 0
      slaoffset_sum = 0. _d 0
      slaoffset_weight = 0. _d 0

        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
               if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) 
     &            ) then
                  slaoffset=slaoffset+RA(i,j,bi,bj)*
     &                      anom_psMslaobs(i,j,bi,bj)
                  slaoffset_sum=slaoffset_sum+RA(i,j,bi,bj)
c Of interest is the total weight for an average of N indep sample (not the area weighted average of weight)
c Since slaoffset is itself area weighted, however, the total weight is only approx the simple weight sum :
                  slaoffset_weight=slaoffset_weight+
     &                      gencost_weight(i,j,bi,bj,igen_lsc)
               endif
              enddo
            enddo
          enddo
        enddo

      _GLOBAL_SUM_RL( slaoffset     , mythid )
      _GLOBAL_SUM_RL( slaoffset_weight , mythid )
      _GLOBAL_SUM_RL( slaoffset_sum , mythid )

      if (slaoffset_sum .eq. 0.0) then
        _BEGIN_MASTER( mythid )
        write(msgbuf,'(a)') ' sshv4: slaoffset_sum = zero!'
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
c        stop   '  ... stopped in cost_ssh.'
      else
        _BEGIN_MASTER( mythid )
        write(msgbuf,'(a,3d22.15)') ' sshv4:slaoffset,sum,weight=',
     &         slaoffset,slaoffset_sum,slaoffset_weight
        call PRINT_MESSAGE( msgbuf, standardmessageunit,
     &                          SQUEEZE_RIGHT , mythid)
        _END_MASTER( mythid )
      endif

c--   Compute (average) slaoffset
      slaoffset = slaoffset / slaoffset_sum

c--   Subtract slaoffset from anom_psMslaobs
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
               if ( ( anom_psMslaobs_NUM(i,j,bi,bj) .NE. 0. ).AND.
     &              ( maskc(i,j,1,bi,bj) .NE. 0. ) 
     &            ) then
                  anom_psMslaobs(i,j,bi,bj) =
     &                 anom_psMslaobs(i,j,bi,bj) - slaoffset
                  anom_slaobs(i,j,bi,bj) =
     &                 anom_slaobs(i,j,bi,bj) - slaoffset
               else
                  anom_psMslaobs(i,j,bi,bj) = 0. _d 0
                  anom_slaobs(i,j,bi,bj) = 0. _d 0
               endif
              enddo
            enddo
          enddo
        enddo

c PART 4.2: smooth anom_psMslaobs in space
c ----------------------------------------

#ifdef ALLOW_GENCOST_SSHV4_OUTPUT
      write(fname4test(1:80),'(1a)') 'sladiff_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMslaobs,irec,1,mythid)

      write(fname4test(1:80),'(1a)') 'slaobs_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_slaobs,irec,1,mythid)
#endif

#ifdef ALLOW_SMOOTH
      if ( useSMOOTH )
     &  call SMOOTH_HETERO2D(anom_psMslaobs,maskc,
     &     gencost_scalefile(igen_lsc),
     &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
#endif

#ifdef ALLOW_GENCOST_SSHV4_OUTPUT
#ifdef ALLOW_SMOOTH
      if ( useSMOOTH )
     &  call SMOOTH_HETERO2D(anom_slaobs,maskc,
     &     gencost_scalefile(igen_lsc),
     &     gencost_smooth2Ddiffnbt(igen_lsc),mythid)
#endif

      write(fname4test(1:80),'(1a)') 'sladiff_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMslaobs,irec,1,mythid)

      write(fname4test(1:80),'(1a)') 'slaobs_smooth'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_slaobs,irec,1,mythid)
#endif

c PART 4.3: compute cost function term
c ------------------------------------

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
# if (defined (ALLOW_SSH_GFOANOM_COST_CONTRIBUTION)  
      defined (ALLOW_SSH_TPANOM_COST_CONTRIBUTION)  
      defined (ALLOW_SSH_ERSANOM_COST_CONTRIBUTION))
          if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
     &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) ) then
            junk = anom_psMslaobs(i,j,bi,bj)
            anom_slaobs_SUB(i,j,bi,bj) = anom_slaobs(i,j,bi,bj)
            anom_psMslaobs_SUB(i,j,bi,bj) = anom_psMslaobs(i,j,bi,bj)
          else
            junk = 0. _d 0
            anom_slaobs_SUB(i,j,bi,bj)= 0. _d 0
            anom_psMslaobs_SUB(i,j,bi,bj)= 0. _d 0
          endif
          objf_gencost(bi,bj,igen_lsc) = objf_gencost(bi,bj,igen_lsc)
     &        + junk*junk*gencost_weight(i,j,bi,bj,igen_lsc)
          if ( (gencost_weight(i,j,bi,bj,igen_lsc).GT.0.).AND.
     &         (anom_psMslaobs_MSK(i,j,bi,bj).GT.0.) )
     &         num_gencost(bi,bj,igen_lsc) =
     &         num_gencost(bi,bj,igen_lsc) + 1. _d 0
#endif
           enddo
         enddo
        enddo
       enddo

#ifdef ALLOW_GENCOST_SSHV4_OUTPUT
      write(fname4test(1:80),'(1a)') 'sladiff_subsample'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMslaobs_SUB,irec,1,mythid)

      write(fname4test(1:80),'(1a)') 'slaobs_subsample'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_slaobs_SUB,irec,1,mythid)
#endif

c PART 4.4: compute cost function term for global mean sea level
c --------------------------------------------------------------

      IF ( ( MASTER_CPU_THREAD(myThid).AND.(igen_gmsl.NE.0) ).AND.
     &     ( slaoffset_sum.GT.0.25 _d 0 * globalArea ) ) THEN
          junk=slaoffset
          junkweight=slaoffset_weight
          objf_gencost(1,1,igen_gmsl) = objf_gencost(1,1,igen_gmsl)
     &        + junk*junk*junkweight
          num_gencost(1,1,igen_gmsl) = num_gencost(1,1,igen_gmsl)
     &        + 1. _d 0
c      print*,'igen_gmsl',igen_gmsl,
      ENDIF

cgf =======================================================
cgf PART 5: compute raw SLA cost term
cgf =======================================================

c time associated with this ndaysrec mean
        krec = irec + (ndaysave/2)

        call ACTIVE_READ_XY( fname, psbar, krec, doglobalread,
     &                       ladinit, optimcycle, mythid,
     &                       xx_psbar_mean_dummy )

#ifndef ALLOW_PSBAR_MEAN
        CALL REMOVE_MEAN_RL( 1, psbar, maskInC, maskInC, rA, drF,
     &        'psbar', myTime, myThid )
#endif

#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
        call COST_READTOPEX( krec, mythid )
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
        call COST_READERS( krec, mythid )
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
        call COST_READGFO( krec, mythid )
#endif

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
              anom_psMtpobs(i,j,bi,bj)  = 0. _d 0
              anom_psMersobs(i,j,bi,bj) = 0. _d 0
              anom_psMgfoobs(i,j,bi,bj) = 0. _d 0
              anom_tpobs(i,j,bi,bj)  = 0. _d 0
              anom_ersobs(i,j,bi,bj) = 0. _d 0
              anom_gfoobs(i,j,bi,bj) = 0. _d 0
          enddo
         enddo
        enddo
       enddo

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      if ( tpmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
     & then
         anom_psMtpobs(i,j,bi,bj)=
     &       psbar(i,j,bi,bj) - tpobs(i,j,bi,bj)
     &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
     &       - anom_psMslaobs(i,j,bi,bj)
         anom_tpobs(i,j,bi,bj)=tpobs(i,j,bi,bj) - slaoffset
     &       - anom_slaobs(i,j,bi,bj)
      endif
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      if ( ersmask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
     & then
         anom_psMersobs(i,j,bi,bj)=
     &       psbar(i,j,bi,bj) - ersobs(i,j,bi,bj)
     &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
     &       - anom_psMslaobs(i,j,bi,bj)
         anom_ersobs(i,j,bi,bj)=ersobs(i,j,bi,bj) - slaoffset
     &       - anom_slaobs(i,j,bi,bj)
      endif
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      if ( gfomask(i,j,bi,bj)*mean_psMslaobs_MSK(i,j,bi,bj).NE.0. )
     & then
         anom_psMgfoobs(i,j,bi,bj)=
     &       psbar(i,j,bi,bj) - gfoobs(i,j,bi,bj)
     &       - mean_psMslaobs(i,j,bi,bj) - slaoffset
     &       - anom_psMslaobs(i,j,bi,bj)
         anom_gfoobs(i,j,bi,bj)=gfoobs(i,j,bi,bj) - slaoffset
     &       - anom_slaobs(i,j,bi,bj)
      endif
#endif
          enddo
         enddo
        enddo
       enddo

#ifdef ALLOW_GENCOST_SSHV4_OUTPUT
      write(fname4test(1:80),'(1a)') 'sladiff_tp_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMtpobs,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'sladiff_ers_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMersobs,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'sladiff_gfo_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_psMgfoobs,irec,1,mythid)

      write(fname4test(1:80),'(1a)') 'slaobs_tp_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_tpobs,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'slaobs_ers_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_ersobs,irec,1,mythid)
      write(fname4test(1:80),'(1a)') 'slaobs_gfo_raw'
      call MDSWRITEFIELD(fname4test,32,.false.,'RL',
     & 1,anom_gfoobs,irec,1,mythid)
#endif

       do bj = jtlo,jthi
        do bi = itlo,ithi
         do j = jmin,jmax
          do i = imin,imax
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
c--             The array psobs contains SSH anomalies.
                junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
     &                      *gencost_weight(i,j,bi,bj,igen_tp)
     &                      *tpmask(i,j,bi,bj)
     &                      *cosphi(i,j,bi,bj)
                junk = anom_psMtpobs(i,j,bi,bj)
                objf_gencost(bi,bj,igen_tp) =
     &              objf_gencost(bi,bj,igen_tp)+junk*junk*junkweight
                if ( junkweight .ne. 0. )
     &              num_gencost(bi,bj,igen_tp) =
     &              num_gencost(bi,bj,igen_tp) + 1. _d 0
#endif
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
c--             The array ersobs contains SSH anomalies.
                junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
     &                      *gencost_weight(i,j,bi,bj,igen_ers)
     &                      *ersmask(i,j,bi,bj)
     &                      *cosphi(i,j,bi,bj)
                junk = anom_psMersobs(i,j,bi,bj)
                objf_gencost(bi,bj,igen_ers) =
     &               objf_gencost(bi,bj,igen_ers)+junk*junk*junkweight
                if ( junkweight .ne. 0. )
     &              num_gencost(bi,bj,igen_ers) =
     &              num_gencost(bi,bj,igen_ers) + 1. _d 0
#endif
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
c--             The array gfoobs contains SSH anomalies.
                junkweight = mean_psMslaobs_MSK(i,j,bi,bj)
     &                      *gencost_weight(i,j,bi,bj,igen_gfo)
     &                      *gfomask(i,j,bi,bj)
     &                      *cosphi(i,j,bi,bj)
                junk = anom_psMgfoobs(i,j,bi,bj)
                objf_gencost(bi,bj,igen_gfo) =
     &              objf_gencost(bi,bj,igen_gfo)+junk*junk*junkweight
                if ( junkweight .ne. 0. )
     &              num_gencost(bi,bj,igen_gfo) =
     &              num_gencost(bi,bj,igen_gfo) + 1. _d 0
#endif
           enddo
         enddo
        enddo
       enddo

      enddo


#endif /* ifdef ALLOW_GENCOST_CONTRIBUTION */
#endif /* ifdef ALLOW_SSH_COST_CONTRIBUTION */

      end