C $Header: /u/gcmpack/MITgcm/pkg/gmredi/gmredi_slope_limit.F,v 1.37 2015/12/29 20:30:31 gforget Exp $
C $Name: $
#include "GMREDI_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
CBOP
C !ROUTINE: GMREDI_SLOPE_LIMIT
C !INTERFACE:
SUBROUTINE GMREDI_SLOPE_LIMIT(
O SlopeX, SlopeY,
O SlopeSqr, taperFct,
U hTransLay, baseSlope, recipLambda,
U dSigmaDr,
I dSigmaDx, dSigmaDy,
I Lrho, hMixLay, depthZ, kLow,
I k, bi, bj, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE GMREDI_SLOPE_LIMIT
C | o Calculate slopes for use in GM/Redi tensor
C *==========================================================*
C | On entry:
C | dSigmaDr contains the d/dz Sigma
C | dSigmaDx/Dy contains X/Y gradients of sigma
C | depthZ contains the depth (< 0 !) [m]
C | Lrho
C | hMixLay mixed layer depth (> 0)
C U hTransLay transition layer depth (> 0)
C U baseSlope, recipLambda,
C | On exit:
C | dSigmaDr contains the effective dSig/dz
C | SlopeX/Y contains X/Y slopes
C | SlopeSqr contains Sx^2+Sy^2
C | taperFct contains tapering funct. value ;
C | = 1 when using no tapering
C U hTransLay transition layer depth (> 0)
C U baseSlope, recipLambda
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ==
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "GMREDI.h"
#include "PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */
C == Routine arguments ==
C !INPUT PARAMETERS:
C dSigmaDx :: zonal gradient of density
C dSigmaDy :: meridional gradient of density
C Lrho ::
C hMixLay :: mixed layer thickness (> 0) [m]
C depthZ :: model discretized depth (< 0) [m]
C kLow :: bottom level index 2-D array
C k :: level index
C bi, bj :: tile indices
C myTime :: time in simulation
C myIter :: iteration number in simulation
C myThid :: My Thread Id. number
C !OUTPUT PARAMETERS:
C SlopeX :: isopycnal slope in X direction
C SlopeY :: isopycnal slope in Y direction
C SlopeSqr :: square of isopycnal slope
C taperFct :: tapering function
C hTransLay :: depth of the base of the transition layer (> 0) [m]
C baseSlope :: slope at the the base of the transition layer
C recipLambda :: Slope vertical gradient at Trans. Layer Base (=recip.Lambda)
C dSigmaDr :: vertical gradient of density
_RL SlopeX (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SlopeY (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SlopeSqr (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL taperFct (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hTransLay (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL baseSlope (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL recipLambda(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDr (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dSigmaDy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL Lrho (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL hMixLay (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS depthZ(*)
INTEGER kLow (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
INTEGER k, bi,bj
_RL myTime
INTEGER myIter
INTEGER myThid
CEOP
#ifdef ALLOW_GMREDI
C !LOCAL VARIABLES:
C == Local variables ==
#ifdef GMREDI_WITH_STABLE_ADJOINT
_RL slopeSqTmp,slopeTmp,slopeMax
#endif
_RL dSigmMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL SlopeMod(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL dRdSigmaLtd(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RL f1,Smod,f2,Rnondim
_RL maxSlopeSqr
_RL fpi
PARAMETER( fpi = PI )
INTEGER i,j
_RL dTransLay, rLambMin, DoverLamb
_RL taperFctLoc, taperFctHat
_RL minTransLay
C- need to put this one in GM namelist :
_RL GM_bigSlope
GM_bigSlope = 1. _d +02
#ifdef ALLOW_AUTODIFF_TAMC
C TAF thinks for some reason that depthZ is an active variable.
C While this does not make the adjoint code wrong, the resulting
C code inhibits vectorization in some cases so we tell TAF here
C that depthZ is actually a passive grid variable that needs no adjoint
CADJ PASSIVE depthz
act1 = bi - myBxLo(myThid)
max1 = myBxHi(myThid) - myBxLo(myThid) + 1
act2 = bj - myByLo(myThid)
max2 = myByHi(myThid) - myByLo(myThid) + 1
act3 = myThid - 1
max3 = nTx*nTy
act4 = ikey_dynamics - 1
ikey = (act1 + 1) + act2*max1
& + act3*max1*max2
& + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
dSigmMod(i,j) = 0. _d 0
tmpFld(i,j) = 0. _d 0
ENDDO
ENDDO
IF (GM_taper_scheme.EQ.'orig' .OR.
& GM_taper_scheme.EQ.'clipping') THEN
#ifdef GM_EXCLUDE_CLIPPING
STOP 'Need to compile without "#define GM_EXCLUDE_CLIPPING"'
#else /* GM_EXCLUDE_CLIPPING */
C- Original implementation in mitgcmuv
C (this turns out to be the same as Cox slope clipping)
C- Cox 1987 "Slope clipping"
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j)
& + dSigmaDy(i,j)*dSigmaDy(i,j)
IF ( tmpFld(i,j) .EQ. 0. ) THEN
dSigmMod(i,j) = 0. _d 0
ELSE
dSigmMod(i,j) = SQRT( tmpFld(i,j) )
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE dSigmMod(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF (dSigmMod(i,j) .NE. 0.) THEN
tmpFld(i,j) = -dSigmMod(i,j)*GM_rMaxSlope
IF ( dSigmaDr(i,j) .GE. tmpFld(i,j) )
& dSigmaDr(i,j) = tmpFld(i,j)
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF (dSigmMod(i,j) .EQ. 0.) THEN
SlopeX(i,j) = 0. _d 0
SlopeY(i,j) = 0. _d 0
ELSE
dRdSigmaLtd(i,j) = 1. _d 0/( dSigmaDr(i,j) )
SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
SlopeSqr(i,j)=SlopeX(i,j)*SlopeX(i,j)
& +SlopeY(i,j)*SlopeY(i,j)
taperFct(i,j)=1. _d 0
ENDDO
ENDDO
#endif /* GM_EXCLUDE_CLIPPING */
ELSEIF (GM_taper_scheme.EQ.'fm07' ) THEN
C-- Ferrari & Mc.Williams, 2007:
#ifdef GM_EXCLUDE_FM07_TAP
STOP 'Need to compile without "#define GM_EXCLUDE_FM07_TAP"'
#else /* GM_EXCLUDE_FM07_TAP */
C- a) Calculate separately slope magnitude (SlopeMod)
C and slope horiz. direction (Slope_X,Y : normalized)
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( k.GT.kLow(i,j) ) THEN
C- Bottom or below:
SlopeX (i,j) = 0. _d 0
SlopeY (i,j) = 0. _d 0
SlopeMod(i,j) = 0. _d 0
taperFct(i,j) = 0. _d 0
ELSE
C- Above bottom:
IF ( dSigmaDr(i,j).GE. -GM_Small_Number )
& dSigmaDr(i,j) = -GM_Small_Number
tmpFld(i,j) = dSigmaDx(i,j)*dSigmaDx(i,j)
& + dSigmaDy(i,j)*dSigmaDy(i,j)
IF ( tmpFld(i,j).GT.0. ) THEN
locVar(i,j) = SQRT( tmpFld(i,j) )
SlopeX (i,j) = dSigmaDx(i,j)/locVar(i,j)
SlopeY (i,j) = dSigmaDy(i,j)/locVar(i,j)
SlopeMod(i,j) = -locVar(i,j)/dSigmaDr(i,j)
taperFct(i,j) = 1. _d 0
ELSE
SlopeX (i,j) = 0. _d 0
SlopeY (i,j) = 0. _d 0
SlopeMod(i,j) = 0. _d 0
taperFct(i,j) = 0. _d 0
ENDIF
ENDIF
ENDDO
ENDDO
C- b) Set Transition Layer Depth:
IF ( k.EQ.1 ) THEN
minTransLay = GM_facTrL2dz*( depthZ(k) - depthZ(k+1) )
ELSE
minTransLay = GM_facTrL2dz*( depthZ(k-1) - depthZ(k) )
ENDIF
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( hTransLay(i,j).LE.0. _d 0 ) THEN
C- previously below the transition layer
tmpFld(i,j) = Lrho(i,j)*SlopeMod(i,j)
C- ensure that transition layer is at least as thick than current level and
C not thicker than the larger of the 2 : maxTransLay and facTrL2ML*hMixLay
dTransLay =
& MIN( MAX( tmpFld(i,j), minTransLay ),
& MAX( GM_facTrL2ML*hMixLay(i,j), GM_maxTransLay ) )
IF ( k.GE.kLow(i,j) ) THEN
C- bottom & below & 1rst level above the bottom :
recipLambda(i,j) = 0. _d 0
baseSlope(i,j) = SlopeMod(i,j)
C-- note: do not catch the 1rst level/interface (k=kLow) above the bottom
C since no baseSlope has yet been stored (= 0); but might fit
C well into transition layer criteria (if locally not stratified)
ELSEIF ( dTransLay+hMixLay(i,j)+depthZ(k) .GE. 0. ) THEN
C- Found the transition Layer : set depth of base of Transition layer
hTransLay(i,j) = -depthZ(k+1)
C and compute inverse length scale "1/Lambda" from slope vert. grad
IF ( baseSlope(i,j).GT.0. ) THEN
recipLambda(i,j) = recipLambda(i,j)
& / MIN( baseSlope(i,j), GM_maxSlope )
ELSE
recipLambda(i,j) = 0. _d 0
ENDIF
C slope^2 & Kwz should remain > 0 : prevent too large negative D/lambda
IF ( hMixLay(i,j)+depthZ(k+1).LT.0. ) THEN
rLambMin = 1. _d 0 /( hMixLay(i,j)+depthZ(k+1) )
recipLambda(i,j) = MAX( recipLambda(i,j), rLambMin )
ENDIF
ELSE
C- Still below Trans. layer: store slope & slope vert. grad.
recipLambda(i,j) = ( MIN( SlopeMod(i,j), GM_maxSlope )
& - MIN( baseSlope(i,j), GM_maxSlope )
& ) / ( depthZ(k) - depthZ(k+1) )
baseSlope(i,j) = SlopeMod(i,j)
ENDIF
ENDIF
ENDDO
ENDDO
C- c) Set Slope component according to vertical position
C (in Mixed-Layer / in Transition Layer / below Transition Layer)
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( hTransLay(i,j).GT.0. _d 0 ) THEN
C- already above base of transition layer:
DoverLamb = (hTransLay(i,j)-hMixLay(i,j))*recipLambda(i,j)
IF ( -depthZ(k).LE.hMixLay(i,j) ) THEN
C- compute tapering function -G(z) in the mixed Layer:
taperFctLoc =
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
& *( 2. _d 0 + DoverLamb )
& )
C- compute tapering function -G^(z) in the mixed Layer
taperFctHat =
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
& * 2. _d 0
& *( 1. _d 0 + DoverLamb )
& )
ELSE
C- compute tapering function -G(z) in the transition Layer:
taperFctLoc =
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
& *( 2. _d 0 + DoverLamb )
& )
& - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j))
& /( hTransLay(i,j)*hTransLay(i,j)
& - hMixLay(i,j)*hMixLay(i,j) )
& *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j) )
& )
C- compute tapering function -G^(z) in the transition Layer:
taperFctHat =
& ( -depthZ(k)/(hTransLay(i,j)+hMixLay(i,j))
& * 2. _d 0
& *( 1. _d 0 + DoverLamb )
& )
& - ( (depthZ(k)+hMixLay(i,j))*(depthZ(k)+hMixLay(i,j))
& /( hTransLay(i,j)*hTransLay(i,j)
& - hMixLay(i,j)*hMixLay(i,j) )
& *( 1. _d 0 + hTransLay(i,j)*recipLambda(i,j)*2. _d 0 )
& )
ENDIF
C- modify the slope (above bottom of transition layer):
c Smod = baseSlope(i,j)
C- safer to limit the slope (even if it might never exceed GM_maxSlope)
Smod = MIN( baseSlope(i,j), GM_maxSlope )
SlopeX(i,j) = SlopeX(i,j)*Smod*taperFctLoc
SlopeY(i,j) = SlopeY(i,j)*Smod*taperFctLoc
c SlopeSqr(i,j) = Smod*Smod*taperFctHat
c SlopeSqr(i,j) = baseSlope(i,j)*Smod*taperFctHat
SlopeSqr(i,j) = MIN( baseSlope(i,j), GM_bigSlope )
& *Smod*taperFctHat
ELSE
C-- Still below base of transition layer:
c Smod = SlopeMod(i,j)
C- safer to limit the slope:
Smod = MIN( SlopeMod(i,j), GM_maxSlope )
SlopeX(i,j) = SlopeX(i,j)*Smod
SlopeY(i,j) = SlopeY(i,j)*Smod
c SlopeSqr(i,j) = Smod*Smod
c SlopeSqr(i,j) = SlopeMod(i,j)*Smod
SlopeSqr(i,j) = MIN( SlopeMod(i,j), GM_bigSlope )
& *Smod
C-- end if baseSlope > 0 / else => above/below base of Trans. Layer
ENDIF
ENDDO
ENDDO
#endif /* GM_EXCLUDE_FM07_TAP */
ELSE IF (GM_taper_scheme.EQ.'ac02') THEN
#ifdef GM_EXCLUDE_AC02_TAP
STOP 'Need to compile without "#define GM_EXCLUDE_AC02_TAP"'
#else /* GM_EXCLUDE_AC02_TAP */
C- New Scheme (A. & C. 2002): relax part of the small slope approximation
C compute the true slope (no approximation)
C but still neglect Kxy & Kyx (assumed to be zero)
maxSlopeSqr = GM_maxSlope*GM_maxSlope
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
dRdSigmaLtd(i,j)=
& dSigmaDx(i,j)*dSigmaDx(i,j)
& + dSigmaDy(i,j)*dSigmaDy(i,j)
& + dSigmaDr(i,j)*dSigmaDr(i,j)
taperFct(i,j) = 1. _d 0
IF (dRdSigmaLtd(i,j).NE.0.) THEN
dRdSigmaLtd(i,j)=1. _d 0
& / ( dRdSigmaLtd(i,j) )
SlopeSqr(i,j)=(dSigmaDx(i,j)*dSigmaDx(i,j)
& +dSigmaDy(i,j)*dSigmaDy(i,j))*dRdSigmaLtd(i,j)
SlopeX(i,j)=-dSigmaDx(i,j)
& *dRdSigmaLtd(i,j)*dSigmaDr(i,j)
SlopeY(i,j)=-dSigmaDy(i,j)
& *dRdSigmaLtd(i,j)*dSigmaDr(i,j)
cph T11(i,j)=(dSigmaDr(i,j)**2)*dRdSigmaLtd(i,j)
ENDIF
#ifndef ALLOWW_AUTODIFF_TAMC
cph-- this part does not adjoint well
IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
ELSE IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN
taperFct(i,j) = 0. _d 0
ENDIF
#endif
ENDDO
ENDDO
#endif /* GM_EXCLUDE_AC02_TAP */
ELSE
#ifdef GM_EXCLUDE_TAPERING
STOP 'Need to compile without "#define GM_EXCLUDE_TAPERING"'
#else /* GM_EXCLUDE_TAPERING */
C----------------------------------------------------------------------
C- Compute the slope, no clipping, but avoid reverse slope in negatively
C stratified (Sigma_Z > 0) region :
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( dSigmaDr(i,j) .NE. 0. ) THEN
IF (dSigmaDr(i,j).GE.(-GM_Small_Number))
& dSigmaDr(i,j) = -GM_Small_Number
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE dSigmaDx(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDy(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE dSigmaDr(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( dSigmaDr(i,j) .EQ. 0. ) THEN
IF ( dSigmaDx(i,j) .NE. 0. ) THEN
SlopeX(i,j) = SIGN( GM_bigSlope, dSigmaDx(i,j) )
ELSE
SlopeX(i,j) = 0. _d 0
ENDIF
IF ( dSigmaDy(i,j) .NE. 0. ) THEN
SlopeY(i,j) = SIGN( GM_bigSlope, dSigmaDy(i,j) )
ELSE
SlopeY(i,j) = 0. _d 0
ENDIF
ELSE
dRdSigmaLtd(i,j) = 1. _d 0 / dSigmaDr(i,j)
SlopeX(i,j)=-dSigmaDx(i,j)*dRdSigmaLtd(i,j)
SlopeY(i,j)=-dSigmaDy(i,j)*dRdSigmaLtd(i,j)
c SlopeX(i,j) = -dSigmaDx(i,j)/dSigmaDr(i,j)
c SlopeY(i,j) = -dSigmaDy(i,j)/dSigmaDr(i,j)
ENDIF
ENDDO
ENDDO
#ifdef ALLOW_AUTODIFF_TAMC
cnostore CADJ STORE slopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
cnostore CADJ STORE slopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
#endif
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
& +SlopeY(i,j)*SlopeY(i,j)
taperFct(i,j) = 1. _d 0
IF ( SlopeSqr(i,j) .GT. GM_slopeSqCutoff ) THEN
slopeSqr(i,j) = GM_slopeSqCutoff
taperFct(i,j) = 0. _d 0
ENDIF
ENDDO
ENDDO
C- Compute the tapering function for the GM+Redi tensor :
IF (GM_taper_scheme.EQ.'linear') THEN
C- Simplest adiabatic tapering = Smax/Slope (linear)
maxSlopeSqr = GM_maxSlope*GM_maxSlope
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
taperFct(i,j) = 1. _d 0
ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
taperFct(i,j) = SQRT(maxSlopeSqr / SlopeSqr(i,j))
slopeSqr(i,j) = MIN( slopeSqr(i,j),GM_bigSlope*GM_bigSlope )
ENDIF
ENDDO
ENDDO
ELSEIF (GM_taper_scheme.EQ.'gkw91') THEN
C- Gerdes, Koberle and Willebrand, Clim. Dyn. 1991
maxSlopeSqr = GM_maxSlope*GM_maxSlope
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
taperFct(i,j) = 1. _d 0
ELSE IF ( SlopeSqr(i,j) .GT. maxSlopeSqr .AND.
& SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
taperFct(i,j) = maxSlopeSqr/SlopeSqr(i,j)
ENDIF
ENDDO
ENDDO
ELSEIF (GM_taper_scheme.EQ.'dm95') THEN
C- Danabasoglu and McWilliams, J. Clim. 1995
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF ( SlopeSqr(i,j) .EQ. 0. ) THEN
taperFct(i,j) = 1. _d 0
ELSE IF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
Smod=SQRT(SlopeSqr(i,j))
taperFct(i,j)=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
ENDIF
ENDDO
ENDDO
ELSEIF (GM_taper_scheme.EQ.'ldd97') THEN
C- Large, Danabasoglu and Doney, JPO 1997
DO j=1-OLy+1,sNy+OLy-1
DO i=1-OLx+1,sNx+OLx-1
IF (SlopeSqr(i,j) .EQ. 0.) THEN
taperFct(i,j) = 1. _d 0
ELSEIF ( SlopeSqr(i,j) .LT. GM_slopeSqCutoff ) THEN
Smod=SQRT(SlopeSqr(i,j))
f1=op5*( 1. _d 0 + tanh( (GM_Scrit-Smod)/GM_Sd ))
Rnondim= -depthZ(k)/(Lrho(i,j)*Smod)
IF ( Rnondim.GE.1. _d 0 ) THEN
f2 = 1. _d 0
ELSE
f2 = op5*( 1. _d 0 + SIN( fpi*(Rnondim-op5) ))
ENDIF
taperFct(i,j)=f1*f2
ENDIF
ENDDO
ENDDO
ELSEIF (GM_taper_scheme.EQ.'stableGmAdjTap') THEN
#ifndef GMREDI_WITH_STABLE_ADJOINT
STOP 'Need to compile wth "#define GMREDI_WITH_STABLE_ADJOINT"'
#else /* GMREDI_WITH_STABLE_ADJOINT */
c special choice for adjoint/optimization of parameters
c (~ strong clipping, reducing non linearity of kw=f(K))
slopeMax= 2. _d -3
CADJ STORE SlopeX(:,:) = comlev1_bibj, key=ikey, byte=isbyte
CADJ STORE SlopeY(:,:) = comlev1_bibj, key=ikey, byte=isbyte
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx
slopeSqTmp=SlopeX(i,j)*SlopeX(i,j)
& +SlopeY(i,j)*SlopeY(i,j)
IF ( slopeSqTmp .GT. slopeMax**2 ) then
slopeTmp=sqrt(slopeSqTmp)
SlopeX(i,j)=SlopeX(i,j)*slopeMax/slopeTmp
SlopeY(i,j)=SlopeY(i,j)*slopeMax/slopeTmp
ENDIF
SlopeSqr(i,j) = SlopeX(i,j)*SlopeX(i,j)
& +SlopeY(i,j)*SlopeY(i,j)
taperFct(i,j) = 1. _d 0
ENDDO
ENDDO
#endif /* GMREDI_WITH_STABLE_ADJOINT */
ELSEIF (GM_taper_scheme.NE.' ') THEN
STOP 'GMREDI_SLOPE_LIMIT: Bad GM_taper_scheme'
ENDIF
#endif /* GM_EXCLUDE_TAPERING */
ENDIF
#endif /* ALLOW_GMREDI */
RETURN
END