C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_gencost_transp.F,v 1.4 2015/11/12 13:22:14 atn Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"

      subroutine COST_GENCOST_TRANSP(mythid)

c     ==================================================================
c     SUBROUTINE cost_gencost_transp
c     ==================================================================
c
c     o Evaluate cost function contributions of section transport.
c
c     ==================================================================
c     SUBROUTINE cost_gencost_boxmean
c     ==================================================================

      implicit none

c     == global variables ==

#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif

c     == routine arguments ==
      integer mythid

#ifdef ALLOW_GENCOST_CONTRIBUTION

c     == local variables ==

      integer nnzobs, nnzbar
      parameter (nnzbar = Nr, nnzobs = Nr)
      integer nrecloc, localrec
      integer localstartdate(4)

      _RL myobs     (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL mybar     (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL localdif  (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL difmask   (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL localweight(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      _RL localtmp   (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)

      _RL dummy,facCost,facNum
c- facNum is 1 divided by the number of tiles in SIZE dot h
      parameter  ( facNum = 1./float(nPx) )
      _RL localperiod
      _RL spminloc, spmaxloc, spzeroloc

      _RL tmpMeanTile(nSx,nSy),tmpNumTile(nSx,nSy)
      _RL tmpMeanGlo,tmpNumGlo

      character*(MAX_LEN_FNAM) mybarfile
      character*(MAX_LEN_FNAM) myobsfile

      integer kgen(NGENCOST3D)
      integer bi,bj
      integer i,j,k
      integer itlo,ithi
      integer jtlo,jthi
      integer obsrec,irec,jrec
      integer il,k2
      integer icount,icount_transp
      logical dosumsq, dovarwei, doreadobs

      integer preproc_i(NGENPPROC)
      _RL preproc_r(NGENPPROC)
      character*(MAX_LEN_FNAM) preproc(NGENPPROC)
      character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)


      logical doglobalread
      logical ladinit
      character*(MAX_LEN_MBUF) msgbuf
      character*(128) fname1, fname0

      logical exst

c     == external functions ==

      integer  ilnblnk
      external 

c     == end of interface ==

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

c=============== PART 0: initilization ===================

c-- detect the relevant gencost indices
      do k=1,NGENCOST3D
        kgen(k)=0
      enddo

c-- write a report of how many transport costs
      write(msgbuf,'(A)') 'Inside cost_gencost_transp:'
      CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                    SQUEEZE_RIGHT, myThid )

      icount_transp=0
      do k=1,NGENCOST
        if ( (gencost_name(k)(1:6).EQ.'transp').AND.
     &     (using_gencost(k)) ) then
          icount_transp=icount_transp+1
          kgen(icount_transp)=k
          il=ilnblnk(gencost_barfile(kgen(icount_transp)))
          write(msgbuf,'(A,i4,A,A)') 'Cost ',kgen(icount_transp),
     &    ': ',gencost_barfile(kgen(icount_transp))(1:il)
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )

        endif
      enddo

c-- write a report of how many transport costs

      do icount=1,icount_transp
      if (kgen(icount).NE.0) then

c ========

c-- initialize objf and num:
        DO bj=jtlo,jthi
         DO bi=itlo,ithi
            objf_gencost(bi,bj,kgen(icount))=0. _d 0
            num_gencost(bi,bj,kgen(icount))=0. _d 0
         ENDDO
        ENDDO

c--   Initialise local variables.
        nrecloc=0
        nrecloc=gencost_nrec(kgen(icount))

c-- only enters if there is at least 1 record 
        if(nrecloc.gt.0) then

          facCost=1. _d 0 / float(nrecloc)

          localperiod=0.
          localperiod=gencost_period(kgen(icount))
          dummy=gencost_dummy(kgen(icount))
          spminloc=gencost_spmin(kgen(icount))
          spmaxloc=gencost_spmax(kgen(icount))
          spzeroloc=gencost_spzero(kgen(icount))

c prefer to have preproc match nosumsq but can not seem to get syntax
c to work for comparison of characters so match dosumsq for now.
          dosumsq=.FALSE.
          dovarwei=.FALSE.
          do k2 = 1, NGENPPROC
            preproc(k2)=gencost_preproc(k2,kgen(icount))
            preproc_i(k2)=gencost_preproc_i(k2,kgen(icount))
            preproc_c(k2)=gencost_preproc_c(k2,kgen(icount))
            preproc_r(k2)=gencost_preproc_r(k2,kgen(icount))
            if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE.
            if (preproc(k2)(1:7).EQ.'dosumsq') dosumsq=.TRUE.
          enddo

c-- report of dosumsq flag to make sure it is false
          il=ilnblnk(gencost_name(kgen(icount)))
          write(msgbuf,'(3A,L5,2A)') 
     &    'Cost ',gencost_name(kgen(icount))(1:il),
     &    ' dosumsq: ',dosumsq,' preproc(1): ',preproc(1)(1:7)
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )

c-- set mybarfile
          mybarfile=gencost_barfile(kgen(icount))

c-- set obsfile if defined. Not use for now. send warning
          doreadobs=.FALSE.
          if( .not. gencost_datafile(kgen(icount)).eq.' ') then
c            doreadobs=.TRUE.
            write(msgbuf,'(A)') 
     &      '**WARNING** S/R COST_GENCOST_TRANSP: gencost_datafile '
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
            write(msgbuf,'(A)') 
     &      'are currently ignored. Adjust the S/R to add code.'
            CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )
          endif

c model mask[W,S]: already included in transp calc in ecco_phys

c--------- PART 0.1 read weights --------------------

c-- read mask in errfile, already stored under msktrVol[W,S]
c-- for now, assume non-time-variable mask

c=============== PART 1: main loop ===================
          do irec = 1,nrecloc

c--------- PART 1.1 read barfiles ------------------

c-- set all bars to zeros:
            call ECCO_ZERO(mybar,Nr,zeroRL,myThid)

c gencost_errfile and fname0 are dummy, get fname1 from mybarfile
            exst=.FALSE.
            call ECCO_ZERO(localtmp,Nr,zeroRL,myThid)
            call COST_GENCAL(mybarfile,gencost_errfile(kgen(icount)),
     &       irec,localstartdate,localperiod,fname1,
     &       fname0,localrec,obsrec,exst,mythid)
            call COST_GENREAD(fname1,mybar,localtmp,irec,nnzbar,
     &       nrecloc,preproc,preproc_c,preproc_i,preproc_r,
     &       dummy,mythid)


c--------- PART 1.2 read data --------------------

c-- ignore for now, but use doreadobs flag if needed
c-- be careful of recomputation when put inside if-end block
c            if(doreadobs) then
            call ECCO_ZERO(myobs,Nr,zeroRL,myThid)
c            endif


c--------- PART 1.3 Cost calculation -------------


c-- keep total at each irec to print out for time-series
            DO bj = jtlo,jthi
              DO bi = itlo,ithi
                tmpMeanTile(bi,bj) = 0. _d 0
                tmpNumTile(bi,bj) = 0. _d 0
              ENDDO
            ENDDO

c compute obs minus bar (localdif) and mask (difmask) 
c note localtmp is set to 1.
            call ECCO_ZERO(localtmp,nnzobs,oneRL,mythid)
            call ECCO_ZERO(localdif,nnzobs,zeroRL,mythid)
            call ECCO_ZERO(difmask,nnzobs,zeroRL,mythid)

c take care to set sp[min,max,zero]loc carefully to not
c filter out signal.  Can consider skip diffmsk step
            call ECCO_DIFFMSK(
     I       mybar, nnzbar, myobs, nnzobs, localtmp,
     I       spminloc, spmaxloc, spzeroloc,
     O       localdif, difmask,
     I       myThid )

            call ECCO_ZERO(localtmp,nnzobs,oneRL,mythid)
            call ECCO_ADDCOST(
     I       localdif,localtmp,difmask,nnzobs,dosumsq,
     O       tmpMeanTile,tmpNumTile,
     I       mythid)

cc either use ecco_diffmsk and ecco_addcost from above or simplify
cc to call below. For now keep syntax consistent with gencost
c            call ecco_zero(localtmp,nnzobs,oneRL,mythid)
c            call ecco_addcost(
c     I       mybar,localtmp,localtmp,nnzobs,dosumsq,
c     O       tmpMeanTile,tmpNumTile,
c     I       mythid)

c global sums for display of time series
c note tmpNumGlo is the constant total wet points in the msktrVol[W,S]
            tmpMeanGlo = 0. _d 0
            tmpNumGlo = 0. _d 0
            il=ilnblnk(gencost_barfile(kgen(icount)))
            CALL GLOBAL_SUM_TILE_RL( tmpMeanTile, tmpMeanGlo, myThid )
            CALL GLOBAL_SUM_TILE_RL( tmpNumTile, tmpNumGlo, myThid )
            WRITE(msgBuf,'(2A,I3,A,1PE21.14,1PE21.14)')
     &        'globalsum transp ',gencost_barfile(kgen(icount))(1:il),
     &        irec,' ',tmpMeanGlo,tmpNumGlo
            CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                          SQUEEZE_RIGHT, myThid )

c sum that is actually be used in cost function
            DO bj = jtlo,jthi
              DO bi = itlo,ithi
                objf_gencost(bi,bj,kgen(icount))=
     &             objf_gencost(bi,bj,kgen(icount))+tmpMeanTile(bi,bj)
              ENDDO
            ENDDO

          enddo !irec

c-- last step: 
c-- divide by number of record to get mean transport:
c-- make num_gencost equals number of months/days used
          do bj = jtlo,jthi
            do bi = itlo,ithi
              objf_gencost(bi,bj,kgen(icount))=
     &           objf_gencost(bi,bj,kgen(icount))*facCost
              num_gencost(bi,bj,kgen(icount))=nrecloc*facNum
            enddo
          enddo

        endif !if (nrecloc.gt.0)
      endif !if (kgen.NE.0)
      enddo !icount_transp

#endif /* ALLOW_GENCOST_CONTRIBUTION */

      RETURN
      end