C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ecco/ecco_cost_init_fixed.F,v 1.39 2014/10/09 21:47:00 gforget Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"
#include "AD_CONFIG.h"

      subroutine ecco_cost_init_fixed( mythid )

c     ==================================================================
c     SUBROUTINE ecco_cost_init_fixed
c     ==================================================================
c
c     o Set contributions to the cost function and the cost function
c       itself to zero. The cost function and the individual contribu-
c       tions are defined in the header file "ecco_cost.h".
c
c     started: Christian Eckert eckert@mit.edu 30-Jun-1999
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 18-Jan-2001
c
c              - move namelist reading to cost_readparms.F
c
c     ==================================================================
c     SUBROUTINE ecco_cost_init_fixed
c     ==================================================================

      implicit none

c     == global variables ==

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

#ifdef ALLOW_CAL
#include "cal.h"
#endif
#ifdef ALLOW_ECCO
c# ifdef ALLOW_OLD_ESTIM_CODES
#  include "ecco_cost.h"
c# else
c#  include "ecco.h"
c# endif
#endif
#ifdef ALLOW_CTRL
# include "optim.h"
#endif

c     == routine arguments ==

      integer mythid

c     == local variables ==

      integer k
      logical exst
c#ifdef ALLOW_OLD_ESTIM_CODES
      _RL     dummy
      integer tempDate1(4)
      integer tempDate2(4)
      integer gwunit
      integer il, ilo,ihi
      character*(max_len_mbuf) msgbuf
      integer irec
      _RL     missingObsFlag
      PARAMETER ( missingObsFlag = 1. _d 23 )
c#endif
#ifdef ALLOW_GENCOST_CONTRIBUTION
      integer bi,bj
      integer i,j,kk,k2
#endif

c     == external functions ==

      integer  cal_IntYears
      external cal_IntYears
      integer  cal_IntMonths
      external cal_IntMonths
      integer  cal_IntDays
      external cal_IntDays
      integer  ifnblnk
      external ifnblnk
      integer  ilnblnk
      external ilnblnk

c     == end of interface ==

#ifdef ALLOW_CTRL
      eccoiter=optimcycle
#else
      eccoiter=0
#endif

#ifdef ALLOW_CAL

c--   The number of monthly and daily averages generated by the
c--   current model integration.
      nyearsrec = cal_IntYears( mythid )
      nmonsrec = cal_IntMonths( mythid )
      ndaysrec = cal_IntDays( mythid )

      _BEGIN_MASTER( myThid )

c#ifdef ALLOW_OLD_ESTIM_CODES

c--     Get the complete dates of the ...
c--     ... TMI data.
        if ( tmidatfile .ne. ' ' )
     &   call cal_FullDate( tmistartdate1,   tmistartdate2,
     &                      tmistartdate,    mythid )
c--     ... SST data.
        if ( sstdatfile .ne. ' ' )
     &   call cal_FullDate( sststartdate1,   sststartdate2,
     &                      sststartdate,    mythid )
c--     ... SSS data.
        if ( sssdatfile .ne. ' ' )
     &   call cal_FullDate( sssstartdate1,   sssstartdate2,
     &                      sssstartdate,    mythid )
c--     ... BP data.
        if ( bpdatfile .ne. ' ' )
     &   call cal_FullDate( bpstartdate1,   bpstartdate2,
     &                      bpstartdate,    mythid )
c--     ... IES data.
        if ( iesdatfile .ne. ' ' )
     &   call cal_FullDate( iesstartdate1,   iesstartdate2,
     &                      iesstartdate,    mythid )
#ifdef ALLOW_SSH_MEAN_COST_CONTRIBUTION
c--     ... mdt data.
        if ( mdtdatfile .ne. ' ' )
     &  call cal_FullDate( mdtstartdate1, mdtstartdate2,
     &                     mdtstartdate,  mythid )
c--     ... mdt data.
        if ( mdtdatfile .ne. ' ' )
     &  call cal_FullDate( mdtenddate1, mdtenddate2,
     &                     mdtenddate,  mythid )
#endif /* ALLOW_SSH_MEAN_COST_CONTRIBUTION */
c--     ... T/P data.
        if ( topexfile .ne. ' ' )
     &  call cal_FullDate( topexstartdate1, topexstartdate2,
     &                     topexstartdate,  mythid )
c--     ... ERS data.
        if ( ersfile .ne. ' ' )
     &   call cal_FullDate( ersstartdate1,   ersstartdate2,
     &                      ersstartdate,    mythid )
c--     ... GFO data.
        if ( gfofile .ne. ' ' )
     &   call cal_FullDate( gfostartdate1,   gfostartdate2,
     &                      gfostartdate,    mythid )
c--     ... SCAT data.
        if ( scatxdatfile .ne. ' ' )
     &   call cal_FullDate( scatstartdate1,   scatstartdate2,
     &                      scatxstartdate,    mythid )
        if ( scatydatfile .ne. ' ' )
     &   call cal_FullDate( scatstartdate1,   scatstartdate2,
     &                      scatystartdate,    mythid )
c--     ... ARGO data.
        if ( argotfile .ne. ' ' )
     &  call cal_FullDate( argotstartdate1,   argotstartdate2,
     &                     argotstartdate,    mythid )
         if ( argosfile .ne. ' ' )
     &  call cal_FullDate( argosstartdate1,   argotstartdate2,
     &                     argosstartdate,    mythid )
c#endif /* ALLOW_OLD_ESTIM_CODES */

#ifdef ALLOW_GENCOST_CONTRIBUTION
      do k = 1, NGENCOST

c--    skip averaging when several cost terms use the 
c      same barfile or when barfile is undefined
       gencost_barskip(k)=.FALSE.
       if ( gencost_barfile(k).EQ.' ' ) 
     &      gencost_barskip(k)=.TRUE.
       do k2 = 1,k-1
         if ( gencost_barfile(k2).EQ.gencost_barfile(k) )
     &      gencost_barskip(k)=.TRUE.
       enddo

c--    set time averaging parameters
       if ( (using_gencost(k)).AND.( (gencost_flag(k).GE.1).OR.
     &         (gencost_avgperiod(k).NE.'     ') ) ) then
         if ( gencost_avgperiod(k) .EQ. 'day' .OR.
     &        gencost_avgperiod(k) .EQ. 'DAY' ) then
            gencost_nrec(k)   = ndaysrec
            gencost_period(k) = 86400.
         else if ( gencost_avgperiod(k) .EQ. 'month' .OR.
     &        gencost_avgperiod(k) .EQ. 'MONTH' ) then
            gencost_nrec(k)   =nmonsrec
            gencost_period(k) = 0.
         else if ( gencost_avgperiod(k) .EQ. 'year' .OR.
     &        gencost_avgperiod(k) .EQ. 'YEAR' ) then
           STOP
     &      'ecco_cost_init_fixed: yearly data not yet implemented'
         else
           STOP
     &      'ecco_cost_init_fixed: gencost_avgperiod wrongly specified'
         endif
         if ( gencost_nrecperiod(k).EQ.0)
     &     gencost_nrecperiod(k)=gencost_nrec(k)
       endif

c--    set observation start/enddate
       if (gencost_startdate1(k).GT.0) then
         call cal_FullDate(
     &     gencost_startdate1(k), gencost_startdate2(k),
     &     gencost_startdate(1,k), mythid )
       else
         call cal_CopyDate(modelStartDate,
     &     gencost_startdate(1,k),mythid)
         gencost_startdate1(k)=startdate_1
         gencost_startdate2(k)=startdate_2
       endif

       if (gencost_enddate1(k).GT.0) then
         call cal_FullDate(
     &     gencost_enddate1(k), gencost_enddate2(k),
     &     gencost_enddate(1,k), mythid )
       else
         call cal_CopyDate(modelEndDate,
     &     gencost_enddate(1,k),mythid)
       endif

c--   set weights
      if ( gencost_errfile(k) .NE. ' ' ) then

      kk=gencost_pointer3d(k)
      il = ilnblnk(gencost_errfile(k))
      inquire( file=gencost_errfile(k)(1:il), exist=exst )
      if (.NOT.exst) using_gencost(k)=.FALSE.

      if ( (.NOT. gencost_timevaryweight(k)).AND.
     &     (.NOT.gencost_is3d(k)).AND.(exst) ) then

           call mdsreadfield( gencost_errfile(k), cost_iprec,
     &     cost_yftype, 1, gencost_weight(:,:,1,1,k), 1, mythid )

         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j = 1-Oly,sNy+Oly
            DO i = 1-Olx,sNx+Olx
c--           initialize weight to 0 if needed
              if( gencost_timevaryweight(k)) then
                gencost_weight(i,j,bi,bj,k) = 0. _d 0
              endif
c--           Test for missing values.
              if (gencost_weight(i,j,bi,bj,k) .lt. -9900.) then
                gencost_weight(i,j,bi,bj,k) = 0. _d 0
              endif
c--           Convert to weight
              if (gencost_weight(i,j,bi,bj,k) .ne. 0.) then
                 gencost_weight(i,j,bi,bj,k) =
     &                1./gencost_weight(i,j,bi,bj,k)/
     &                   gencost_weight(i,j,bi,bj,k)
              endif
            enddo
           enddo
          enddo
         enddo

#ifdef ALLOW_GENCOST3D
      elseif ( (.NOT. gencost_timevaryweight(k)).AND.
     &     (gencost_is3d(k)).AND.(kk.LE.NGENCOST3D).AND.(exst) ) then

           call mdsreadfield( gencost_errfile(k), cost_iprec,
     &     cost_yftype, nr, gencost_wei3d(:,:,1,1,1,kk), 1, mythid )

         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j = 1-Oly,sNy+Oly
            DO i = 1-Olx,sNx+Olx
             DO k2 = 1,nr
c--           initialize weight to 0 if needed
              if( gencost_timevaryweight(k)) then
                gencost_wei3d(i,j,k2,bi,bj,kk) = 0. _d 0
              endif
c--           Test for missing values.
              if (gencost_wei3d(i,j,k2,bi,bj,kk) .lt. -9900.) then
                gencost_wei3d(i,j,k2,bi,bj,kk) = 0. _d 0
              endif
c--           Convert to weight
              if (gencost_wei3d(i,j,k2,bi,bj,kk) .ne. 0.) then
                 gencost_wei3d(i,j,k2,bi,bj,kk) =
     &                1./gencost_wei3d(i,j,k2,bi,bj,kk)/
     &                   gencost_wei3d(i,j,k2,bi,bj,kk)
              endif
             enddo
            enddo
           enddo
          enddo
         enddo
#endif /* ALLOW_GENCOST3D */

      endif !if ( (.NOT. gencost_timevaryweight(k)).AND.
      endif !if ( gencost_errfile(k) .NE. ' ' ) then
      enddo !do k = 1, NGENCOST
#endif /* ALLOW_GENCOST_CONTRIBUTION */

      _END_MASTER( mythid )

#endif /* ALLOW_CAL */

      call ecco_check( myThid )

c--   Get the weights that are to be used for the individual cost
c--   function contributions.
      call ecco_cost_weights( mythid )

c--   Initialise adjoint of monthly mean files calculated
c--   in cost_averagesfields (and their ad...).
cph(
cph   The following init. shoud not be applied if in the middle
cph   of a divided adjoint run
cph)
#ifndef ALLOW_TANGENTLINEAR_RUN
cph!!! and I think it needs to be seen by TAF
cph!!! for repeated TLM runs
cph!!!
      inquire( file='costfinal', exist=exst )
      if ( .NOT. exst) then
         call ecco_cost_init_barfiles( mythid )
      endif
#endif

c#ifdef ALLOW_OLD_ESTIM_CODES

#ifdef ALLOW_TRANSPORT_COST_CONTRIBUTION
      do irec = 1, ndaysrec
       wtransp(irec)   = 0. _d 0
       transpobs(irec) = 0. _d 0
      enddo

      if ( costTranspDataFile .NE. ' ' ) then
      _BEGIN_MASTER(myThid)
      ilo = ifnblnk(costTranspDataFile)
      ihi = ilnblnk(costTranspDataFile)
      CALL OPEN_COPY_DATA_FILE(
     I                          costTranspDataFile(ilo:ihi),
     I                          'ECCO_COST_INIT_FIXED',
     O                          gwunit,
     I                          myThid )
      do irec = 1, ndaysrec
c-- read daily transport time series
c-- 1st: transport in m/s
c-- 2nd: date in YYYYMMDD
c-- 3rd: uncertainty in m/s
         read(gwunit,*) transpobs(irec), dummy, wtransp(irec)
c-- convert std.dev. to weight
         if ( wtransp(irec) .NE. 0. )
     &        wtransp(irec) =1.0/(wtransp(irec)*wtransp(irec))
c-- set weight to zero for missing values
         if ( transpobs(irec) .EQ. missingObsFlag )
     &        wtransp(irec) = 0. _d 0
      enddo
      _END_MASTER(myThid)
      _BARRIER
      endif
#endif /* ALLOW_TRANSPORT_COST_CONTRIBUTION */

#ifdef ALLOW_NEW_SSH_COST

c--   Read flags for picking SSH time averages
      do irec = 1, ndaysrec
       tpTimeMask(irec)  = 1. _d 0
       ersTimeMask(irec) = 1. _d 0
       gfoTimeMask(irec) = 1. _d 0
      enddo
c
      _BEGIN_MASTER(myThid)
c
#ifdef ALLOW_SSH_TPANOM_COST_CONTRIBUTION
      if ( tpTimeMaskFile .NE. ' ' ) then
       ilo = ifnblnk(tpTimeMaskFile)
       ihi = ilnblnk(tpTimeMaskFile)
       CALL OPEN_COPY_DATA_FILE(
     I                          tpTimeMaskFile(ilo:ihi),
     I                          'cost_ssh tp',
     O                          gwunit,
     I                          myThid )
       do irec = 1, ndaysrec
        read(gwunit,*) tpTimeMask(irec)
       enddo
      endif
#endif
c
#ifdef ALLOW_SSH_ERSANOM_COST_CONTRIBUTION
      if ( ersTimeMaskFile .NE. ' ' ) then
       ilo = ifnblnk(ersTimeMaskFile)
       ihi = ilnblnk(ersTimeMaskFile)
       CALL OPEN_COPY_DATA_FILE(
     I                          ersTimeMaskFile(ilo:ihi),
     I                          'cost_ssh ers',
     O                          gwunit,
     I                          myThid )
       do irec = 1, ndaysrec
        read(gwunit,*) ersTimeMask(irec)
       enddo
      endif
#endif
c
#ifdef ALLOW_SSH_GFOANOM_COST_CONTRIBUTION
      if ( gfoTimeMaskFile .NE. ' ' ) then
       ilo = ifnblnk(gfoTimeMaskFile)
       ihi = ilnblnk(gfoTimeMaskFile)
       CALL OPEN_COPY_DATA_FILE(
     I                          gfoTimeMaskFile(ilo:ihi),
     I                          'cost_ssh gfo',
     O                          gwunit,
     I                          myThid )
       do irec = 1, ndaysrec
        read(gwunit,*) gfoTimeMask(irec)
       enddo
      endif
#endif
c
       do irec = 1, ndaysrec
        if (
     &   ( tpTimeMask(irec).NE.0. .AND. tpTimeMask(irec).NE.1. ) .OR.
     &   ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) .OR.
     &   ( ersTimeMask(irec).NE.0. .AND. ersTimeMask(irec).NE.1. ) )
     &  then
           WRITE(msgBuf,'(2A,I10)')
     &        'ecco_cost_init_fixed: (SSH)TimeMask not 0. or 1. ',
     &        'for irec (=day) ', irec
           CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
     &          SQUEEZE_RIGHT , myThid )
           CALL PRINT_ERROR( msgBuf , myThid )
           STOP 'ABNORMAL END: S/R ECCO_COST_INIT_FIXED'
        endif
       enddo
c
      _END_MASTER(myThid)
      _BARRIER
#endif /* ALLOW_NEW_SSH_COST */

c#endif /* ALLOW_OLD_ESTIM_CODES */

c--   Summarize the cost function setup.
      _BEGIN_MASTER( mythid )
      call ecco_summary( mythid )
      call ecco_cost_summary( mythid )
      _END_MASTER( mythid )

      _BARRIER

      end