C $Header: /u/gcmpack/MITgcm/model/src/forcing_surf_relax.F,v 1.1 2013/07/04 23:02:48 jmc Exp $
C $Name: $
#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: FORCING_SURF_RELAX
C !INTERFACE:
SUBROUTINE FORCING_SURF_RELAX(
I iMin, iMax, jMin, jMax,
I myTime, myIter, myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE FORCING_SURF_RELAX
C | o Calculate relaxation surface forcing terms
C | for temperature and salinity
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "SURFACE.h"
#ifdef ALLOW_SEAICE
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# include "SEAICE.h"
#endif /* ALLOW_SEAICE */
C !INPUT/OUTPUT PARAMETERS:
C === Routine arguments ===
C iMin,iMax, jMin,jMax :: Range of points for calculation
C myTime :: Current time in simulation
C myIter :: Current iteration number in simulation
C myThid :: Thread no. that called this routine.
INTEGER iMin, iMax
INTEGER jMin, jMax
_RL myTime
INTEGER myIter
INTEGER myThid
C !LOCAL VARIABLES:
C === Local variables ===
C bi,bj :: tile indices
C i,j :: loop indices
C ks :: index of surface interface layer
INTEGER bi,bj
INTEGER i,j
INTEGER ks
CEOP
#ifdef ALLOW_DIAGNOSTICS
_RL tmpFac
#endif /* ALLOW_DIAGNOSTICS */
#ifdef ALLOW_BALANCE_RELAX
CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL sumTile(nSx,nSy), sumGlob, globAver
#endif /* ALLOW_BALANCE_RELAX */
IF ( usingPCoords ) THEN
ks = Nr
ELSE
ks = 1
ENDIF
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_SEAICE
IF ( useSEAICE .AND. (.NOT. SEAICErestoreUnderIce) ) THEN
C Do not restore under sea-ice
DO j = jMin, jMax
DO i = iMin, iMax
C Heat Flux (restoring term) :
surfaceForcingT(i,j,bi,bj) =
& -lambdaThetaClimRelax(i,j,bi,bj)*(1.-AREA(i,j,bi,bj))
& *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
& *drF(ks)*_hFacC(i,j,ks,bi,bj)
C Salt Flux (restoring term) :
surfaceForcingS(i,j,bi,bj) =
& -lambdaSaltClimRelax(i,j,bi,bj) *(1.-AREA(i,j,bi,bj))
& *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
& *drF(ks)*_hFacC(i,j,ks,bi,bj)
ENDDO
ENDDO
ELSE
#endif /* ALLOW_SEAICE */
DO j = jMin, jMax
DO i = iMin, iMax
C Heat Flux (restoring term) :
surfaceForcingT(i,j,bi,bj) =
& -lambdaThetaClimRelax(i,j,bi,bj)
& *(theta(i,j,ks,bi,bj)-SST(i,j,bi,bj))
& *drF(ks)*_hFacC(i,j,ks,bi,bj)
C Salt Flux (restoring term) :
surfaceForcingS(i,j,bi,bj) =
& -lambdaSaltClimRelax(i,j,bi,bj)
& *(salt(i,j,ks,bi,bj)-SSS(i,j,bi,bj))
& *drF(ks)*_hFacC(i,j,ks,bi,bj)
ENDDO
ENDDO
#ifdef ALLOW_SEAICE
ENDIF
#endif /* ALLOW_SEAICE */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef NONLIN_FRSURF
C- T,S surface forcing will be applied (thermodynamics) after the update
C of surf.thickness (hFac): account for change in surf.thickness
IF (staggerTimeStep.AND.nonlinFreeSurf.GT.0) THEN
IF ( select_rStar.GT.0 ) THEN
# ifndef DISABLE_RSTAR_CODE
DO j=jMin,jMax
DO i=iMin,iMax
surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
& * rStarExpC(i,j,bi,bj)
surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
& * rStarExpC(i,j,bi,bj)
ENDDO
ENDDO
# endif /* DISABLE_RSTAR_CODE */
ELSEIF ( selectSigmaCoord.NE.0 ) THEN
# ifndef DISABLE_SIGMA_CODE
DO j=jMin,jMax
DO i=iMin,iMax
surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
& *( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
& *dBHybSigF(ks)*recip_drF(ks)
& *recip_hFacC(i,j,ks,bi,bj)
& )
surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
& *( 1. _d 0 + dEtaHdt(i,j,bi,bj)*deltaTFreeSurf
& *dBHybSigF(ks)*recip_drF(ks)
& *recip_hFacC(i,j,ks,bi,bj)
& )
ENDDO
ENDDO
# endif /* DISABLE_SIGMA_CODE */
ELSE
DO j=jMin,jMax
DO i=iMin,iMax
IF (ks.EQ.kSurfC(i,j,bi,bj)) THEN
surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
& *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
& *_recip_hFacC(i,j,ks,bi,bj)*hFac_surfC(i,j,bi,bj)
ENDIF
ENDDO
ENDDO
ENDIF
ENDIF
#endif /* NONLIN_FRSURF */
C-- end bi,bj loops.
ENDDO
ENDDO
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_BALANCE_RELAX
IF ( balanceThetaClimRelax ) THEN
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
sumTile(bi,bj) = 0. _d 0
DO j=1,sNy
DO i=1,sNx
sumTile(bi,bj) = sumTile(bi,bj)
& + surfaceForcingT(i,j,bi,bj)
& *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
globAver = sumGlob
IF ( globalArea.GT.zeroRL ) globAver = globAver / globalArea
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
& - globAver
ENDDO
ENDDO
ENDDO
ENDDO
IF ( balancePrintMean ) THEN
_BEGIN_MASTER( myThid )
WRITE(msgBuf,'(A,E24.17)')
& 'rm Global mean of SSTrelax =', globAver
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
_END_MASTER( myThid )
ENDIF
ENDIF
IF ( balanceSaltClimRelax ) THEN
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
sumTile(bi,bj) = 0. _d 0
DO j=1,sNy
DO i=1,sNx
sumTile(bi,bj) = sumTile(bi,bj)
& + surfaceForcingS(i,j,bi,bj)
& *rA(i,j,bi,bj)*maskInC(i,j,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
CALL GLOBAL_SUM_TILE_RL( sumTile, sumGlob, myThid )
globAver = sumGlob
IF ( globalArea.GT.zeroRL ) globAver = globAver / globalArea
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
& - globAver
ENDDO
ENDDO
ENDDO
ENDDO
IF ( balancePrintMean ) THEN
_BEGIN_MASTER( myThid )
WRITE(msgBuf,'(A,E24.17)')
& 'rm Global mean of SSSrelax =', globAver
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
& SQUEEZE_RIGHT, myThid )
_END_MASTER( myThid )
ENDIF
ENDIF
#endif /* ALLOW_BALANCE_RELAX */
#ifdef ALLOW_DIAGNOSTICS
IF ( useDiagnostics ) THEN
C tRelax (temperature relaxation [W/m2], positive <-> increasing Theta)
tmpFac = HeatCapacity_Cp*rUnit2mass
CALL DIAGNOSTICS_SCALE_FILL( surfaceForcingT, tmpFac, 1,
& 'TRELAX ', 0, 1, 0,1,1, myThid )
C sRelax (salt relaxation [g/m2/s], positive <-> increasing Salt)
tmpFac = rUnit2mass
CALL DIAGNOSTICS_SCALE_FILL( surfaceForcingS, tmpFac, 1,
& 'SRELAX ', 0, 1, 0,1,1, myThid )
ENDIF
#endif /* ALLOW_DIAGNOSTICS */
RETURN
END