C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_cost_gen.F,v 1.9 2015/10/29 03:43:59 gforget Exp $
C $Name:  $

c ----------------------------------------------------------------
c --- ctrl_cost_gen2d
c --- ctrl_cost_gen3d
c ----------------------------------------------------------------

c ----------------------------------------------------------------

#include "CTRL_OPTIONS.h"
#ifdef ALLOW_ECCO
# include "ECCO_OPTIONS.h"
#endif

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: ctrl_cost_gen2d
C     !INTERFACE:
      subroutine CTRL_COST_GEN2D(
     I                       startrec,
     I                       endrec,
     I                       xx_gen_file,
     I                       xx_gen_dummy,
     I                       xx_gen_period,
     I                       xx_gen_weight,
     I                       dodimensionalcost,
     O                       num_gen_anom,
     O                       objf_gen_anom,
#ifdef ECCO_CTRL_DEPRECATED
     I                       xx_gen_wmean,
     O                       num_gen_mean,
     O                       objf_gen_mean,
     O                       objf_gen_smoo,
     I                       xx_gen_remo_intercept,
     I                       xx_gen_remo_slope,
#endif /* ECCO_CTRL_DEPRECATED */
     I                       xx_gen_mask,
     I                       myThid
     &                         )

C     !DESCRIPTION: \bv
C     Generic routine for all 2D control penalty terms
C     \ev

      implicit none

c     == global variables ==

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

#ifdef ALLOW_ECCO
#  include "ecco.h"
#endif
#ifdef ALLOW_CTRL
# include "ctrl.h"
# include "optim.h"
#endif

c     == routine arguments ==

      integer startrec
      integer endrec
      character*(MAX_LEN_FNAM) xx_gen_file
      _RL xx_gen_dummy
      _RL xx_gen_period
      _RL xx_gen_weight(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
      logical dodimensionalcost
      _RL num_gen_anom(nsx,nsy)
      _RL objf_gen_anom(nsx,nsy)
#ifdef ECCO_CTRL_DEPRECATED
      _RL xx_gen_wmean
      _RL num_gen_mean(nsx,nsy)
      _RL num_gen_smoo(nsx,nsy)
      _RL objf_gen_mean(nsx,nsy)
      _RL objf_gen_smoo(nsx,nsy)
      _RL xx_gen_remo_intercept
      _RL xx_gen_remo_slope
#endif /* ECCO_CTRL_DEPRECATED */
      _RS xx_gen_mask(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
      integer myThid

#ifdef ALLOW_CTRL

c     == local variables ==

      integer bi,bj
      integer i,j,kk
      integer itlo,ithi
      integer jtlo,jthi
      integer jmin,jmax
      integer imin,imax
      integer nrec
      integer irec
      integer ilfld

      _RL fctile
      _RL fctilem
      _RL tmpx
      _RL lengthscale

#ifdef ECCO_CTRL_DEPRECATED
      _RL xx_mean(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
#endif /* ECCO_CTRL_DEPRECATED */

      logical doglobalread
      logical ladinit

      character*(80) fnamefld

c     == external functions ==

      integer  ilnblnk
      external 

CEOP

      jtlo = mybylo(myThid)
      jthi = mybyhi(myThid)
      itlo = mybxlo(myThid)
      ithi = mybxhi(myThid)
      jmin = 1
      jmax = sny
      imin = 1
      imax = snx

      lengthscale = 1. _d 0

c--   Read state record from global file.
      doglobalread = .false.
      ladinit      = .false.

c     Number of records to be used.
      nrec = endrec-startrec+1

      if (optimcycle .ge. 0) then
        ilfld=ilnblnk( xx_gen_file )
        write(fnamefld(1:80),'(2a,i10.10)')
     &       xx_gen_file(1:ilfld),'.',optimcycle
      endif

c--   >>> Loop 1 to compute mean forcing:
      do bj = jtlo,jthi
        do bi = itlo,ithi
          num_gen_anom(bi,bj)  = 0. _d 0
          objf_gen_anom(bi,bj) = 0. _d 0
#ifdef ECCO_CTRL_DEPRECATED
          do j = jmin,jmax
            do i = imin,imax
              xx_mean(i,j,bi,bj)   = 0. _d 0
            enddo
          enddo
          num_gen_mean(bi,bj)  = 0. _d 0
          num_gen_smoo(bi,bj)  = 0. _d 0
          objf_gen_mean(bi,bj) = 0. _d 0
          objf_gen_smoo(bi,bj) = 0. _d 0
#endif /* ECCO_CTRL_DEPRECATED */
        enddo
      enddo

#ifndef ECCO_CTRL_DEPRECATED
c--   >>> Loop over records.
      do irec = 1,nrec

#ifdef ALLOW_AUTODIFF
        call ACTIVE_READ_XY(
     &        fnamefld, tmpfld2d, irec, doglobalread,
     &        ladinit, optimcycle, myThid, xx_gen_dummy )
#else
        CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif

c--     Loop over this thread tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

c--         Determine the weights to be used.
            kk = 1
            fctile = 0. _d 0
            do j = jmin,jmax
              do i = imin,imax
                if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then

                tmpx = tmpfld2d(i,j,bi,bj)
                IF ( dodimensionalcost ) THEN
                  fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
                ELSE
                  fctile = fctile + tmpx*tmpx
                ENDIF
                if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
     &            num_gen_anom(bi,bj) = num_gen_anom(bi,bj)
     &            + 1. _d 0
                endif

              enddo
            enddo

            objf_gen_anom(bi,bj) = objf_gen_anom(bi,bj) + fctile

          enddo
        enddo

c--   End of loop over records.
      enddo

#else /* ECCO_CTRL_DEPRECATED */

      IF ( dodimensionalcost ) THEN
      do irec = 1,nrec

#ifdef ALLOW_AUTODIFF
        call ACTIVE_READ_XY(
     &        fnamefld, tmpfld2d, irec, doglobalread,
     &        ladinit, optimcycle, myThid, xx_gen_dummy )
#else
        CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif

c--     Loop over this thread tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi
            do j = jmin,jmax
              do i = imin,imax
                xx_mean(i,j,bi,bj) = xx_mean(i,j,bi,bj)
     &                + tmpfld2d(i,j,bi,bj)
     &                - ( xx_gen_remo_intercept +
     &                    xx_gen_remo_slope*(irec-1)*xx_gen_period )
              enddo
            enddo
          enddo
        enddo

      enddo

      if ( xx_gen_wmean .NE. 0. ) then
       do bj = jtlo,jthi
        do bi = itlo,ithi
c--     Determine the weights to be used.
        kk = 1
        fctilem = 0. _d 0
        do j = jmin,jmax
          do i = imin,imax
            xx_mean(i,j,bi,bj)
     &            = xx_mean(i,j,bi,bj)/float(nrec)
            tmpx = xx_mean(i,j,bi,bj)/xx_gen_wmean
            if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
#ifdef ALLOW_ECCO
              if ( ABS(R_low(i,j,bi,bj)) .LT. 100.  _d 0)
     &              tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100.  _d 0
              fctilem = fctilem + cosphi(i,j,bi,bj)*tmpx*tmpx
              if ( cosphi(i,j,bi,bj) .ne. 0. _d 0)
     &             num_gen_mean(bi,bj) = num_gen_mean(bi,bj) + 1. _d 0
#else
              fctilem = fctilem + tmpx*tmpx
                   num_gen_mean(bi,bj) = num_gen_mean(bi,bj) + 1. _d 0
#endif
            endif
          enddo
        enddo
        objf_gen_mean(bi,bj) = objf_gen_mean(bi,bj) + fctilem
        enddo
       enddo
      endif
      ENDIF !IF ( dodimensionalcost ) THEN

c--   >>> Loop 2 over records.
      do irec = 1,nrec

#ifdef ALLOW_AUTODIFF
        call ACTIVE_READ_XY(
     &        fnamefld, tmpfld2d, irec, doglobalread,
     &        ladinit, optimcycle, myThid, xx_gen_dummy )
#else
        CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif

c--     Loop over this thread tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

c--         Determine the weights to be used.
            kk = 1
            fctile = 0. _d 0
            do j = jmin,jmax
              do i = imin,imax
                if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then

                IF ( dodimensionalcost ) THEN
                  tmpx = tmpfld2d(i,j,bi,bj)
     &                   - xx_mean(i,j,bi,bj)
     &                   - ( xx_gen_remo_intercept +
     &                       xx_gen_remo_slope*(irec-1)*xx_gen_period )
#ifdef ALLOW_ECCO
                  if ( ABS(R_low(i,j,bi,bj)) .LT. 100. _d 0 )
     &              tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100. _d 0
                  fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
     &                   * cosphi(i,j,bi,bj)
#else
                  fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
#endif
                ELSE !IF ( dodimensionalcost ) THEN
                  tmpx = tmpfld2d(i,j,bi,bj)
                  fctile = fctile + tmpx*tmpx
                ENDIF !IF ( dodimensionalcost ) THEN
#ifdef ALLOW_ECCO
                  if ( xx_gen_weight(i,j,bi,bj)
     &                *cosphi(i,j,bi,bj) .ne. 0. _d 0 )
#else
                  if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
#endif
     &                 num_gen_anom(bi,bj) = num_gen_anom(bi,bj)
     &                 + 1. _d 0

                endif

              enddo
            enddo

            objf_gen_anom(bi,bj) = objf_gen_anom(bi,bj) + fctile

          enddo
        enddo

c--   End of loop over records.
      enddo

      IF ( dodimensionalcost ) THEN
#ifdef ALLOW_SMOOTH_BC_COST_CONTRIBUTION

c--   >>> Loop 2 over records.
      do irec = 1,nrec

#ifdef ALLOW_AUTODIFF
        call ACTIVE_READ_XY(
     &        fnamefld, tmpfld2d, irec, doglobalread,
     &        ladinit, optimcycle, myThid, xx_gen_dummy )
#else
        CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif

        _EXCH_XY_RL(tmpfld2d, myThid)

c--     Loop over this thread tiles.
        do bj = jtlo,jthi
          do bi = itlo,ithi

c--         Determine the weights to be used.
            kk = 1
            fctile = 0. _d 0
            do j = jmin,jmax
              do i = imin,imax
                if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
                  tmpx =
     &                 ( tmpfld2d(i+2,j,bi,bj)-tmpfld2d(i+1,j,bi,bj) )
     &                   *maskW(i+1,j,kk,bi,bj)*maskW(i+2,j,kk,bi,bj)
     &               + ( tmpfld2d(i+1,j,bi,bj)-tmpfld2d(i,j,bi,bj) )
     &                   *maskW(i+1,j,kk,bi,bj)
     &               + ( tmpfld2d(i,j+2,bi,bj)-tmpfld2d(i,j+1,bi,bj) )
     &                   *maskS(i,j+1,kk,bi,bj)*maskS(i,j+2,kk,bi,bj)
     &               + ( tmpfld2d(i,j+1,bi,bj)-tmpfld2d(i,j,bi,bj) )
     &                   *maskS(i,j+1,kk,bi,bj)
#ifdef ALLOW_ECCO
                  if ( ABS(R_low(i,j,bi,bj)) .LT. 100. _d 0 )
     &              tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100. _d 0
                  fctile = fctile
     &               + xx_gen_weight(i,j,bi,bj)*cosphi(i,j,bi,bj)
#else
                  fctile = fctile
     &               + xx_gen_weight(i,j,bi,bj)
#endif
     *                 *0.0161 _d 0*lengthscale/4.0 _d 0
     &                 *tmpx*tmpx
#ifdef ALLOW_ECCO
                  if ( xx_gen_weight(i,j,bi,bj)*cosphi(i,j,bi,bj)
     &                 .ne. 0.  _d 0 )
#else
                  if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
#endif
     &                 num_gen_smoo(bi,bj) = num_gen_smoo(bi,bj)
     &                 + 1. _d 0
                endif
              enddo
            enddo

            objf_gen_smoo(bi,bj) = objf_gen_smoo(bi,bj) + fctile

          enddo
        enddo

c--   End of loop over records.
      enddo

#endif /* ALLOW_SMOOTH_BC_COST_CONTRIBUTION */
      ENDIF !IF ( dodimensionalcost ) THEN

#endif /* ECCO_CTRL_DEPRECATED */

#endif /* ALLOW_CTRL */

      return
      end


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: ctrl_cost_gen3d C !INTERFACE: subroutine CTRL_COST_GEN3D( I xx_gen_file, I xx_gen_dummy, I xx_gen_weight, I dodimensionalcost, O num_gen, O objf_gen, I xx_gen_mask, I myThid & ) C !DESCRIPTION: \bv C Generic routine for all 3D control penalty terms C \ev implicit none c == global variables == #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" #ifdef ALLOW_ECCO # include "ecco.h" #endif #ifdef ALLOW_CTRL # include "ctrl.h" # include "optim.h" #endif c == routine arguments == character*(MAX_LEN_FNAM) xx_gen_file _RL xx_gen_dummy _RL xx_gen_weight(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy) logical dodimensionalcost _RL num_gen(nsx,nsy) _RL objf_gen(nsx,nsy) _RS xx_gen_mask(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy) INTEGER myThid #ifdef ALLOW_CTRL c == local variables == integer bi,bj integer i,j,k integer itlo,ithi integer jtlo,jthi integer jmin,jmax integer imin,imax integer irec integer ilfld _RL tmpx logical doglobalread logical ladinit character*(80) fnamefld c == external functions == integer ilnblnk external CEOP jtlo = mybylo(myThid) jthi = mybyhi(myThid) itlo = mybxlo(myThid) ithi = mybxhi(myThid) jmin = 1 jmax = sny imin = 1 imax = snx c-- Read state record from global file. doglobalread = .false. ladinit = .false. if (optimcycle .ge. 0) then ilfld = ilnblnk( xx_gen_file ) write(fnamefld(1:80),'(2a,i10.10)') & xx_gen_file(1:ilfld),'.',optimcycle endif c-- >>> Loop 1 to compute mean forcing: do bj = jtlo,jthi do bi = itlo,ithi num_gen(bi,bj) = 0. _d 0 objf_gen(bi,bj) = 0. _d 0 enddo enddo irec = 1 #ifdef ALLOW_AUTODIFF call ACTIVE_READ_XYZ( fnamefld, tmpfld3d, irec, doglobalread, & ladinit, optimcycle, myThid & , xx_gen_dummy ) #else CALL READ_REC_XYZ_RL( fnamefld, tmpfld3d, iRec, 1, myThid ) #endif c-- Loop over this thread tiles. do bj = jtlo,jthi do bi = itlo,ithi num_gen(bi,bj) = 0. _d 0 objf_gen(bi,bj) = 0. _d 0 do k = 1,nr do j = jmin,jmax do i = imin,imax if (xx_gen_mask(i,j,k,bi,bj) .ne. 0. _d 0) then tmpx = tmpfld3d(i,j,k,bi,bj) IF ( dodimensionalcost ) THEN objf_gen(bi,bj) = objf_gen(bi,bj) & + xx_gen_weight(i,j,k,bi,bj) & *tmpx*tmpx ELSE objf_gen(bi,bj) = objf_gen(bi,bj) + tmpx*tmpx ENDIF if ( xx_gen_weight(i,j,k,bi,bj) .ne. 0. _d 0 ) & num_gen(bi,bj) = num_gen(bi,bj) + 1. _d 0 endif enddo enddo enddo enddo enddo #endif return end