C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_generic.F,v 1.39 2016/03/22 22:29:25 gforget Exp $
C $Name:  $

#include "ECCO_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: cost_generic
C     !INTERFACE:
      subroutine COST_GENERIC(
     &     nnzbar, localbarfile, dummy,
     &     nnzobs, localobsfile, localerrfile,
     &     mult_local, nrecloc, nrecobs,
     &     localstartdate, localperiod,
     &     ylocmask, spminloc, spmaxloc, spzeroloc,
     &     preproc, preproc_c, preproc_i, preproc_r,
     &     posproc, posproc_c, posproc_i, posproc_r,
     &     outlev, outname,
     &     objf_local, num_local,
     &     myiter, mytime, mythid )

C     !DESCRIPTION: \bv
C     Generic routine for evaluating time-dependent
c       cost function contribution
C     \ev

C     !USES:
      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
#ifdef ALLOW_SEAICE
# include "SEAICE_COST.h"
#endif

c     == routine arguments ==

      integer myiter
      integer mythid
      integer nnzbar, nnzobs
      integer nrecloc, nrecobs
      integer localstartdate(4)
      integer outlev
      integer preproc_i(NGENPPROC)
      integer posproc_i(NGENPPROC)

      _RL objf_local(nsx,nsy)
      _RL num_local(nsx,nsy)
      _RL dummy
      _RL mult_local
      _RL mytime
      _RL localperiod
      _RL spminloc
      _RL spmaxloc
      _RL spzeroloc
      _RL preproc_r(NGENPPROC)
      _RL posproc_r(NGENPPROC)

      character*(1) ylocmask
      character*(MAX_LEN_FNAM) localbarfile
      character*(MAX_LEN_FNAM) localobsfile
      character*(MAX_LEN_FNAM) localerrfile
      character*(MAX_LEN_FNAM) preproc(NGENPPROC)
      character*(MAX_LEN_FNAM) preproc_c(NGENPPROC)
      character*(MAX_LEN_FNAM) posproc(NGENPPROC)
      character*(MAX_LEN_FNAM) posproc_c(NGENPPROC)
      character*(MAX_LEN_FNAM) outname

#ifdef ALLOW_ECCO

c     == local variables ==

      integer bi,bj,k2
      integer itlo,ithi
      integer jtlo,jthi
      logical domean, doanom

      _RL localdifmean1  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)
      _RL localdifmean2  (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy)

CEOP

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

c--   Initialise local variables.

      do bj = jtlo,jthi
        do bi = itlo,ithi
          objf_local(bi,bj) = 0. _d 0
          num_local(bi,bj) = 0. _d 0
        enddo
      enddo

      call ECCO_ZERO(localdifmean1,Nr,zeroRL,myThid)
      call ECCO_ZERO(localdifmean2,Nr,zeroRL,myThid)

      domean=.FALSE.
      doanom=.FALSE.
      do k2 = 1, NGENPPROC
          if (preproc(k2).EQ.'mean') domean=.TRUE.
          if (preproc(k2).EQ.'anom') doanom=.TRUE.
      enddo

C Extra time loop to compute time-mean fields and costs
      if ( (.NOT. ( localobsfile.EQ.' ' ) ) 
     &   .AND. ( domean .OR. doanom ) ) then
        call COST_GENLOOP(
     &     localdifmean1,localdifmean2,.FALSE.,
     &     nnzbar, localbarfile, dummy,
     &     nnzobs, localobsfile, localerrfile,
     &     mult_local, nrecloc, nrecobs,
     &     localstartdate, localperiod,
     &     ylocmask, spminloc, spmaxloc, spzeroloc,
     &     preproc, preproc_c, preproc_i, preproc_r,
     &     posproc, posproc_c, posproc_i, posproc_r,
     &     outlev, outname,
     &     objf_local, num_local,
     &     myiter, mytime, mythid )
      endif

      call ECCO_ZERO(localdifmean1,Nr,zeroRL,myThid)

      if ((.NOT.(localobsfile.EQ.' ')).AND.(.NOT.domean)) then
        call COST_GENLOOP(
     &     localdifmean2,localdifmean1,.TRUE.,
     &     nnzbar, localbarfile, dummy,
     &     nnzobs, localobsfile, localerrfile,
     &     mult_local, nrecloc, nrecobs,
     &     localstartdate, localperiod,
     &     ylocmask, spminloc, spmaxloc, spzeroloc,
     &     preproc, preproc_c, preproc_i, preproc_r,
     &     posproc, posproc_c, posproc_i, posproc_r,
     &     outlev, outname,
     &     objf_local, num_local,
     &     myiter, mytime, mythid )
      endif

#endif /* ALLOW_ECCO */

      return
      end


C-------------- subroutine COST_GENLOOP( & localdifmeanIn,localdifmeanOut, addVariaCost, & nnzbar, localbarfile, dummy, & nnzobs, localobsfile, localerrfile, & mult_local, nrecloc, nrecobs, & localstartdate, localperiod, & ylocmask, spminloc, spmaxloc, spzeroloc, & preproc, preproc_c, preproc_i, preproc_r, & posproc, posproc_c, posproc_i, posproc_r, & outlev, outname, & objf_local, num_local, & myiter, mytime, mythid ) C !DESCRIPTION: \bv C Generic routine for evaluating time-dependent c cost function contribution C \ev C !USES: 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 #ifdef ALLOW_SEAICE # include "SEAICE_COST.h" #endif c == routine arguments == integer myiter integer mythid integer nnzbar, nnzobs integer nrecloc, nrecobs integer localstartdate(4) integer outlev integer preproc_i(NGENPPROC) integer posproc_i(NGENPPROC) _RL objf_local(nsx,nsy) _RL num_local(nsx,nsy) _RL dummy _RL mult_local _RL mytime _RL localperiod _RL spminloc _RL spmaxloc _RL spzeroloc _RL preproc_r(NGENPPROC) _RL posproc_r(NGENPPROC) character*(1) ylocmask character*(MAX_LEN_FNAM) localbarfile character*(MAX_LEN_FNAM) localobsfile character*(MAX_LEN_FNAM) localerrfile character*(MAX_LEN_FNAM) preproc(NGENPPROC) character*(MAX_LEN_FNAM) preproc_c(NGENPPROC) character*(MAX_LEN_FNAM) posproc(NGENPPROC) character*(MAX_LEN_FNAM) posproc_c(NGENPPROC) character*(MAX_LEN_FNAM) outname logical addVariaCost _RL localdifmeanIn (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) _RL localdifmeanOut (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) #ifdef ALLOW_ECCO c == local variables == integer bi,bj integer itlo,ithi integer jtlo,jthi integer irec, jrec integer il, k2 integer localrec, obsrec integer nrecloop, nrecclim, k2smooth logical domean, doanom, dovarwei, doclim, dosmooth, dosumsq _RL localmask (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) _RL localbar (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 localobs (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 localdifmsk (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) _RL localdifsum (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) _RL localdifnum (1-olx:snx+olx,1-oly:sny+oly,Nr,nsx,nsy) _RL fac character*(128) fname1, fname2, fname3 character*200 msgbuf logical exst c == external functions == integer ilnblnk external CEOP call ECCO_ZERO(localbar,Nr,zeroRL,myThid) call ECCO_ZERO(localweight,Nr,zeroRL,myThid) call ECCO_ZERO(localtmp,Nr,zeroRL,myThid) call ECCO_ZERO(localmask,Nr,zeroRL,myThid) call ECCO_ZERO(localobs,Nr,zeroRL,myThid) call ECCO_ZERO(localdif,Nr,zeroRL,myThid) call ECCO_ZERO(difmask,Nr,zeroRL,myThid) call ECCO_ZERO(localdifmsk,Nr,zeroRL,myThid) call ECCO_ZERO(localdifsum,Nr,zeroRL,myThid) call ECCO_ZERO(localdifnum,Nr,zeroRL,myThid) dosumsq=.TRUE. domean=.FALSE. doanom=.FALSE. dovarwei=.FALSE. dosmooth=.FALSE. k2smooth=1 doclim=.FALSE. nrecclim=nrecloc fac=oneRL do k2 = 1, NGENPPROC if (preproc(k2).EQ.'mean') domean=.TRUE. if (preproc(k2).EQ.'anom') doanom=.TRUE. if (preproc(k2).EQ.'variaweight') dovarwei=.TRUE. if (preproc(k2).EQ.'nosumsq') dosumsq=.FALSE. if (posproc(k2).EQ.'smooth') then dosmooth=.TRUE. k2smooth=k2 endif if (preproc(k2).EQ.'clim') then doclim=.TRUE. nrecclim=preproc_i(k2) endif if (preproc(k2).EQ.'factor') then fac=preproc_r(k2) endif enddo c-- Assign mask if ( ylocmask .EQ. 'C' .OR. ylocmask .EQ. 'c' ) then call ECCO_CPRSRL(maskC,nr,localmask,nr,myThid) elseif ( ylocmask .EQ. 'S' .OR. ylocmask .EQ. 's' ) then call ECCO_CPRSRL(maskS,nr,localmask,nr,myThid) elseif ( ylocmask .EQ. 'W' .OR. ylocmask .EQ. 'w' ) then call ECCO_CPRSRL(maskW,nr,localmask,nr,myThid) else STOP 'cost_generic: wrong ylocmask' endif c-- set nrecloop to nrecloc nrecloop=nrecloc c-- reset nrecloop, if needed, according to preproc if ( doclim ) nrecloop=MIN(nrecloop,nrecclim) c-- loop over obsfile records do irec = 1, nrecloop c-- load weights exst=.FALSE. jrec=1 if( dovarwei ) jrec = irec call COST_GENCAL(localbarfile, localerrfile, & jrec, localstartdate, localperiod, fname1, & fname3, localrec, obsrec, exst, mythid ) call ECCO_ZERO(localweight,nnzobs,zeroRL,myThid) if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) & call ECCO_READWEI(fname3,localweight,localrec,nnzobs,mythid) c-- determine records and file names exst=.FALSE. call COST_GENCAL(localbarfile, localobsfile, & irec, localstartdate, localperiod, fname1, & fname2, localrec, obsrec, exst, mythid ) c-- load model average and observed average call ECCO_ZERO(localbar,nnzbar,zeroRL,myThid) call COST_GENREAD( fname1, localbar, localtmp, irec, nnzbar, & nrecloc, preproc, preproc_c, preproc_i, preproc_r, & dummy, mythid ) call ECCO_MULT(localbar,nnzbar,fac,myThid) call ECCO_ZERO(localobs,nnzobs,spzeroloc,myThid) if ( (localrec .GT. 0).AND.(obsrec .GT. 0).AND.(exst) ) & call MDSREADFIELD( fname2, cost_iprec, cost_yftype, nnzobs, & localobs, localrec, mythid ) c-- Compute masked model-data difference call ECCO_DIFFMSK( localbar, nnzbar, localobs, nnzobs, & localmask, spminloc, spmaxloc, spzeroloc, & localdif, difmask, myThid ) if ( doanom ) call ECCO_SUBTRACT( localdifmeanIn, & nnzobs, localdif, nnzobs, myThid ) if ( domean.OR.doanom ) & call ECCO_ADDMASK(localdif,difmask, nnzobs,localdifsum, & localdifnum, nnzobs,myThid) if (addVariaCost) then #ifdef ALLOW_SMOOTH if ( useSMOOTH.AND.dosmooth.AND. & (nnzbar.EQ.1).AND.(nnzobs.EQ.1) ) & call SMOOTH_HETERO2D(localdif,maskc, & posproc_c(k2smooth),posproc_i(k2smooth),mythid) #endif c-- Compute normalized model-obs cost function call ECCO_ADDCOST( I localdif, localweight, difmask, nnzobs, dosumsq, U objf_local, num_local, I myThid & ) c-- output model-data difference to disk if ( outlev.GT.0 ) then il=ilnblnk(outname) write(fname3(1:128),'(2a)') 'misfit_', outname(1:il) if ( nnzobs.EQ.1 ) CALL & WRITE_REC_XY_RL( fname3, localdif,irec, eccoiter, mythid ) if ( nnzobs.EQ.nr ) CALL & WRITE_REC_XYZ_RL( fname3, localdif,irec, eccoiter, mythid ) endif endif enddo c-- End of loop over obsfile records. call ECCO_ZERO(localdifmeanOut,Nr,zeroRL,myThid) call ECCO_CP(localdifsum,nnzobs,localdifmeanOut,nnzobs,myThid) call ECCO_DIVFIELD(localdifmeanOut,nnzobs,localdifnum,myThid) call ECCO_CP(localdifnum,nnzobs,localdifmsk,nnzobs,myThid) call ECCO_DIVFIELD(localdifmsk,nnzobs,localdifnum,myThid) if ( domean ) then c-- Compute normalized model-obs cost function call ECCO_ADDCOST( I localdifmeanOut, localweight, localdifmsk, nnzobs, dosumsq, U objf_local, num_local, myThid) c-- output model-data difference to disk if ( outlev.GT.0 ) then il=ilnblnk(outname) write(fname3(1:128),'(2a)') 'misfit_', outname(1:il) if ( nnzobs.EQ.1 ) CALL & WRITE_REC_XY_RL(fname3,localdifmeanOut,1,eccoiter,mythid) if ( nnzobs.EQ.nr ) CALL & WRITE_REC_XYZ_RL(fname3,localdifmeanOut,1,eccoiter,mythid) endif endif if ( outlev.GT.1 ) then il=ilnblnk(outname) write(fname3(1:128),'(2a)') 'weight_', outname(1:il) if ( nnzobs.EQ.1 ) CALL & WRITE_REC_XY_RL( fname3, localweight,irec, eccoiter, mythid ) if ( nnzobs.EQ.nr ) CALL & WRITE_REC_XYZ_RL( fname3, localweight,irec, eccoiter, mythid ) endif #endif /* ALLOW_ECCO */ return end