C $Header: /u/gcmpack/MITgcm/model/src/update_sigma.F,v 1.2 2014/04/04 20:56:32 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
c#ifdef ALLOW_AUTODIFF
c# include "AUTODIFF_OPTIONS.h"
c#endif
CBOP
C !ROUTINE: UPDATE_SIGMA
C !INTERFACE:
SUBROUTINE UPDATE_SIGMA( etaHc, myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE UPDATE_SIGMA
C | o Update the thickness fractions (hFacC,W,S)
C | according to the surface r-position = Non-Linear FrSurf
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
c #include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
c#ifdef ALLOW_AUTODIFF_TAMC
c# include "tamc.h"
c# include "tamc_keys.h"
c#endif
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C etaHc :: surface r-anomaly at grid cell center
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: my Thread Id. number
_RL etaHc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL myTime
INTEGER myIter
INTEGER myThid
#ifdef NONLIN_FRSURF
c#ifndef DISABLE_SIGMA_CODE
C !LOCAL VARIABLES:
C Local variables
C bi, bj :: tile indices
C i, j, k :: Loop counters
C rEmpty :: empty column r-position
C rFullDepth :: maximum depth of a full column
C tmpFld :: Temporary array used to compute & write Total Depth
C msgBuf :: Informational/error message buffer
INTEGER bi, bj
INTEGER i, j, k
_RL rFullDepth
_RL tmpFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c _RL hFactmp
c CHARACTER*(MAX_LEN_MBUF) msgBuf
CEOP
rFullDepth = rF(1)-rF(Nr+1)
DO bj=myByLo(myThid), myByHi(myThid)
DO bi=myBxLo(myThid), myBxHi(myThid)
c#ifdef ALLOW_AUTODIFF_TAMC
c act1 = bi - myBxLo(myThid)
c max1 = myBxHi(myThid) - myBxLo(myThid) + 1
c act2 = bj - myByLo(myThid)
c max2 = myByHi(myThid) - myByLo(myThid) + 1
c act3 = myThid - 1
c max3 = nTx*nTy
c act4 = ikey_dynamics - 1
c idynkey = (act1 + 1) + act2*max1
c & + act3*max1*max2
c & + act4*max1*max2*max3
c#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
c#ifdef ALLOW_OBCS
c# ifdef ALLOW_AUTODIFF_TAMC
cCADJ STORE rStarFacC(:,:,bi,bj) =
cCADJ & comlev1_bibj, key = idynkey, byte = isbyte
cCADJ STORE rStarFacS(:,:,bi,bj) =
cCADJ & comlev1_bibj, key = idynkey, byte = isbyte
cCADJ STORE rStarFacW(:,:,bi,bj) =
cCADJ & comlev1_bibj, key = idynkey, byte = isbyte
c# endif /* ALLOW_AUTODIFF_TAMC */
cC-- Apply OBC to rStar_Factor_W,S before updating hFacW,S
c IF (useOBCS) THEN
c CALL OBCS_APPLY_R_STAR(
c I bi, bj,
c U rStarFacC, rStarFacW, rStarFacS,
c I myTime, myIter, myThid )
c ENDIF
c#endif /* ALLOW_OBCS */
C-- Update the fractional thickness hFacC (& "recip_hFac") :
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF ( kSurfC(i,j,bi,bj).LE.Nr ) THEN
tmpFld(i,j) = etaHc(i,j,bi,bj)
& + ( Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj) )
ELSE
tmpFld(i,j) = rFullDepth
ENDIF
ENDDO
ENDDO
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hFacC(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)
& *( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )*recip_drF(k)
recip_hFacC(i,j,k,bi,bj) = maskC(i,j,k,bi,bj)*drF(k)
& /( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )
ENDDO
ENDDO
ENDDO
C-- Update the fractional thickness hFacW (& "recip_hFac") :
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF ( kSurfW(i,j,bi,bj).LE.Nr ) THEN
tmpFld(i,j) = etaHw(i,j,bi,bj)
& + ( rSurfW(i,j,bi,bj)-rLowW(i,j,bi,bj) )
ELSE
tmpFld(i,j) = rFullDepth
ENDIF
ENDDO
ENDDO
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hFacW(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)
& *( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )*recip_drF(k)
recip_hFacW(i,j,k,bi,bj) = maskW(i,j,k,bi,bj)*drF(k)
& /( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )
ENDDO
ENDDO
ENDDO
C-- Update the fractional thickness hFacS (& "recip_hFac") :
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
IF ( kSurfS(i,j,bi,bj).LE.Nr ) THEN
tmpFld(i,j) = etaHs(i,j,bi,bj)
& + ( rSurfS(i,j,bi,bj)-rLowS(i,j,bi,bj) )
ELSE
tmpFld(i,j) = rFullDepth
ENDIF
ENDDO
ENDDO
DO k=1,Nr
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
hFacS(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)
& *( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )*recip_drF(k)
recip_hFacS(i,j,k,bi,bj) = maskS(i,j,k,bi,bj)*drF(k)
& /( dAHybSigF(k)*rFullDepth
& +dBHybSigF(k)*tmpFld(i,j)
& )
ENDDO
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C- end bi,bj loop
ENDDO
ENDDO
c _EXCH_XYZ_RS( hFacC, myThid )
c _EXCH_XYZ_RS( recip_hFacC, myThid )
c CALL EXCH_UV_XYZ_RS(hFacW,hFacS,.FALSE.,myThid)
c CALL EXCH_UV_XYZ_RS(recip_hFacW,recip_hFacS,.FALSE.,myThid)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
c#endif /* ndef DISABLE_SIGMA_CODE */
#endif /* NONLIN_FRSURF */
RETURN
END