C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/heimbach/admtlm_setup/code_ad_metric/cost_test.F,v 1.1 2005/11/04 19:00:29 heimbach Exp $

#include "CPP_OPTIONS.h"

      subroutine cost_test( myThid )
C     /==========================================================\
C     | subroutine cost_test                                     |
C     | o this routine computes the cost function for the tiles  |
C     |   of this processor                                      |
C     |==========================================================|
C     |                                                          |
C     | Notes                                                    |
C     | =====                                                    |
C     \==========================================================/
      IMPLICIT NONE

C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"

#include "cost.h"

C     == Routine arguments ==
C     myThid - Thread number for this instance of the routine.
      integer bi, bj
      integer myThid

#ifdef ALLOW_COST_TEST
C     == Local variables
      _RL thetaRef
      _RL myMetric
      _RL multAll

      _RL stdDevSST
      _RL stdDevSSS
      _RL stdDevtheta
      _RL stdDevsalt
      _RL stdDevuvel
      _RL stdDevvvel
      _RL stdDevetan

      _RL fcTheta
      _RL fcSalt
      _RL fcUvel
      _RL fcVvel
      _RL fcEtan

      _RL numTheta
      _RL numSalt
      _RL numUvel
      _RL numVvel
      _RL numEtan

      integer i, j, k
      integer ig, jg
      integer itlo,ithi
      integer jtlo,jthi

      jtlo = mybylo(mythid)
      jthi = mybyhi(mythid)
      itlo = mybxlo(mythid)
      ithi = mybxhi(mythid)

      myMetric  = 0. _d 0
      stdDevSST = 0.52 _d 0
      stdDevSSS = 0.134 _d 0
c
      stdDevTheta = 1.
      stdDevSalt  = 1.
      stdDevUvel  = 1.
      stdDevVvel  = 1.
      stdDevEtan  = 1.

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
c
        fcTheta = 0.
        fcSalt = 0.
        fcUvel = 0.
        fcVvel = 0.
        fcEtan = 0.
c
        numTheta = 0.
        numSalt = 0.
        numUvel = 0.
        numVvel = 0.
        numEtan = 0.
c
        DO j=1,sNy
         DO i=1,sNx
          DO k=1,Nr
c
           fcTheta = fcTheta +
     &            hFacC(i,j,k,bi,bj)
     &            *theta(i,j,k,bi,bj)**2/stdDevTheta**2
           numTheta = numTheta + 
     &            maskC(i,j,k,bi,bj)
c
           fcSalt = fcSalt +
     &            hFacC(i,j,k,bi,bj)
     &            *salt(i,j,k,bi,bj)**2/stdDevSalt**2
           numSalt = numSalt + 
     &            hFacC(i,j,k,bi,bj)
c
           fcUvel = fcUvel +
     &            hFacW(i,j,k,bi,bj)
     &            *uvel(i,j,k,bi,bj)**2/stdDevUvel**2
           numUvel = numUvel + 
     &            hFacW(i,j,k,bi,bj)
c
           fcVvel = fcVvel +
     &            hFacS(i,j,k,bi,bj)
     &            *vvel(i,j,k,bi,bj)**2/stdDevVvel**2
           numVvel = numVvel + 
     &            hFacS(i,j,k,bi,bj)
c
           if ( k .EQ. 1 ) then
              fcEtan = fcEtan +
     &             hFacC(i,j,k,bi,bj)
     &             *etan(i,j,bi,bj)**2/stdDevEtan**2
              numEtan = numEtan + 
     &             hFacC(i,j,k,bi,bj)
           endif
c
          END DO
         END DO
        END DO
c
        if ( numTheta .NE. 0. ) 
     &       fcTheta = multTheta*fcTheta/numTheta
        if ( numSalt .NE. 0. ) 
     &       fcSalt = multSalt*fcSalt/numSalt
        if ( numUvel .NE. 0. ) 
     &       fcUvel = multUvel*fcUvel/numUvel
        if ( numVvel .NE. 0. ) 
     &       fcVvel = multVvel*fcVvel/numVvel
        if ( numEtan .NE. 0. ) 
     &       fcEtan = multEtan*fcEtan/numEtan

        multAll = 
     &       multTheta + multSalt + multUvel + multVvel + multEtan

        if ( multAll .NE. 0. ) then
           objf_test(bi,bj) = 1./multAll*
     &       ( fcTheta + fcSalt + fcUvel + fcVvel + fcEtan )
        else
           objf_test(bi,bj) = 0.
        endif
c
       END DO
      END DO

#endif

      END