C $Header: /u/gcmpack/MITgcm/pkg/ctrl/adctrl_bound.F,v 1.11 2012/08/12 19:34:22 jmc Exp $
C $Name:  $

#include "CTRL_OPTIONS.h"

C--  File ctrl_bound.F:
C--   Contents
C--   o ADCTRL_BOUND_3D
C--   o ADCTRL_BOUND_2D

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C     !ROUTINE: ADCTRL_BOUND_3D
C     !INTERFACE:
      SUBROUTINE ADCTRL_BOUND_3D(
     U                fieldCur, adjFieldCur,
     I                maskFld3d, boundsVec, myThid )

C     !DESCRIPTION: \bv
C     *==========================================================*
C     | started: Gael Forget gforget@mit.edu 20-Aug-2007
C     |
C     | o in forward mode: impose bounds on ctrl vector values
C     | o in adjoint mode: do nothing ... or emulate local minimum
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"

C     !INPUT/OUTPUT PARAMETERS:
      _RL fieldCur   (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL adjFieldCur(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RS maskFld3d  (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL boundsVec(5)
      INTEGER myThid

#ifdef ALLOW_ADCTRLBOUND
C     !LOCAL VARIABLES:
      INTEGER bi,bj,i,j,k
      _RL x0,x0p5,l0p5,x1,x2,x2p5,l2p5,x3
      _RL tmpCur,xCur,adxCur
CEOP

      x0=boundsVec(1)
      x1=boundsVec(2)
        x0p5=(x0+x1)/2.0
        l0p5=(x1-x0)/2.0
      x2=boundsVec(3)
      x3=boundsVec(4)
        x2p5=(x2+x3)/2.0
        l2p5=(x3-x2)/2.0

C  x0 ctrl_bound and adctrl_bound   act on xx/adxx
C  x0=x3        => ctrl_bound and adctrl_bound   DO nothing
C  otherwise    => error

      IF ( x0.LT.x3 ) THEN
       IF ( (x0.LT.x1).AND.(x1.LT.x2).AND.(x2.LT.x3) ) THEN

        DO bj=myByLo(myThid), myByHi(myThid)
         DO bi=myBxLo(myThid), myBxHi(myThid)

          DO k = 1,Nr
           DO j = 1,sNy
            DO i = 1,sNx
              IF (maskFld3d(i,j,k,bi,bj).NE.0.) THEN
                xCur=fieldCur(i,j,k,bi,bj)
                adxCur=adjFieldCur(i,j,k,bi,bj)
                IF ( (xCur.GT.x2).AND.(adxCur.LT.0.) ) THEN
                  tmpCur=1.0
                  adjFieldCur(i,j,k,bi,bj)=abs(adxCur)*
     &            min((xCur-x2p5)/l2p5,tmpCur)
                ENDIF
                IF ( (xCur.LT.x1).AND.(adxCur.GT.0.) ) THEN
                  tmpCur=-1.0
                  adjFieldCur(i,j,k,bi,bj)=abs(adxCur)*
     &            max((xCur-x0p5)/l0p5,tmpCur)
                ENDIF
              ENDIF
            ENDDO
           ENDDO
          ENDDO

         ENDDO
        ENDDO

       ELSE
          PRINT*,"boundsVec is not self-consistent"
          STOP 'ABNORMAL END: S/R ADCTRL_BOUND_3D'
       ENDIF
      ENDIF

#endif /* ALLOW_ADCTRLBOUND */

      RETURN
      END


C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| CBOP C !ROUTINE: ADCTRL_BOUND_2D C !INTERFACE: SUBROUTINE ADCTRL_BOUND_2D( U fieldCur, adjFieldCur, I maskFld3d, boundsVec, myThid ) C !DESCRIPTION: \bv C *==========================================================* C | started: Gael Forget gforget@mit.edu 20-Aug-2007 C | C | o in forward mode: impose bounds on ctrl vector values C | o in adjoint mode: DO nothing ... or emulate local minimum C *==========================================================* C \ev C !USES: IMPLICIT NONE #include "SIZE.h" #include "EEPARAMS.h" C !INPUT/OUTPUT PARAMETERS: _RL fieldCur (1-Olx:sNx+Olx,1-Oly:sNy+Oly, nSx,nSy) _RL adjFieldCur(1-Olx:sNx+Olx,1-Oly:sNy+Oly, nSx,nSy) _RS maskFld3d (1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy) _RL boundsVec(5) INTEGER myThid #ifdef ALLOW_ADCTRLBOUND C !LOCAL VARIABLES: INTEGER bi,bj,i,j _RL x0,x0p5,l0p5,x1,x2,x2p5,l2p5,x3 _RL tmpCur,xCur,adxCur CEOP x0=boundsVec(1) x1=boundsVec(2) x0p5=(x0+x1)/2.0 l0p5=(x1-x0)/2.0 x2=boundsVec(3) x3=boundsVec(4) x2p5=(x2+x3)/2.0 l2p5=(x3-x2)/2.0 C x0 ctrl_bound and adctrl_bound act on xx/adxx C x0=x3 => ctrl_bound and adctrl_bound DO nothing C otherwise => error IF ( x0.LT.x3 ) THEN IF ( (x0.LT.x1).AND.(x1.LT.x2).AND.(x2.LT.x3) ) THEN DO bj=myByLo(myThid), myByHi(myThid) DO bi=myBxLo(myThid), myBxHi(myThid) DO j = 1,sNy DO i = 1,sNx IF (maskFld3d(i,j,1,bi,bj).NE.0.) THEN xCur=fieldCur(i,j,bi,bj) adxCur=adjFieldCur(i,j,bi,bj) IF ( (xCur.GT.x2).AND.(adxCur.LT.0.) ) THEN tmpCur=1.0 adjFieldCur(i,j,bi,bj)=abs(adxCur)* & min((xCur-x2p5)/l2p5,tmpCur) ENDIF IF ( (xCur.LT.x1).AND.(adxCur.GT.0.) ) THEN tmpCur=-1.0 adjFieldCur(i,j,bi,bj)=abs(adxCur)* & max((xCur-x0p5)/l0p5,tmpCur) ENDIF ENDIF ENDDO ENDDO ENDDO ENDDO ELSE PRINT*,"boundsVec is not self-consistent" STOP 'ABNORMAL END: S/R ADCTRL_BOUND_2D' ENDIF ENDIF #endif /* ALLOW_ADCTRLBOUND */ RETURN END