C $Header: /u/gcmpack/MITgcm/pkg/layers/layers_thermodynamics.F,v 1.5 2015/06/15 21:34:07 jmc Exp $
C $Name: $
#include "LAYERS_OPTIONS.h"
C-- File layers_thermodynamics.F:
C-- Contents
C-- o LAYERS_CALC_RHS
CBOP 0
C !ROUTINE: LAYERS_CALC_RHS
C !INTERFACE:
SUBROUTINE LAYERS_CALC_RHS(
I myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE LAYERS_CALC_RHS
C | Recalculate the divergence of the RHS terms in T and S eqns.
C | Replaces the values of layers_surfflux, layers_df? IN PLACE
C | with the corresponding tendencies (same units as GT and GS)
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "FFIELDS.h"
#include "LAYERS_SIZE.h"
#include "LAYERS.h"
C !INPUT PARAMETERS:
C myThid :: my Thread Id number
INTEGER myThid
CEOP
#ifdef ALLOW_LAYERS
#ifdef LAYERS_THERMODYNAMICS
C !LOCAL VARIABLES:
C bi, bj :: tile indices
C i,j :: horizontal indices
C k :: vertical index for model grid
C kdown :: temporary placeholder
C fluxfac :: scaling factor for converting surface flux to tendency
C fluxfac :: scaling factor for converting diffusive flux to tendency
C downfac :: mask for lower point
INTEGER bi, bj
INTEGER i,j,k,kdown,iTracer
_RL fluxfac(2), downfac, tmpfac
c CHARACTER*(MAX_LEN_MBUF) msgBuf
_RL minusone
PARAMETER (minusOne=-1.)
#ifdef SHORTWAVE_HEATING
_RL swfracb(2)
#endif
C -- These factors convert the units of TFLUX and SFLUX diagnostics
C -- back to surfaceForcingT and surfaceForcingS units
fluxfac(1) = 1.0/(HeatCapacity_Cp*rUnit2mass)
fluxfac(2) = 1.0/rUnit2mass
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
DO iTracer = 1,2
k = 1
C -- Loop for surface fluxes
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
#ifdef SHORTWAVE_HEATING
C -- Have to remove the shortwave from the surface flux because it is added later
IF (iTracer.EQ.1) THEN
layers_surfflux(i,j,k,iTracer,bi,bj) =
& layers_surfflux(i,j,k,iTracer,bi,bj)
C -- Sign convention for Qsw means we have to add it to subtract it
& +Qsw(i,j,bi,bj)
ENDIF
#endif /* SHORTWAVE_HEATING */
layers_surfflux(i,j,k,iTracer,bi,bj) =
& layers_surfflux(i,j,k,iTracer,bi,bj)
& *recip_drF(1)*_recip_hFacC(i,j,1,bi,bj)
& *fluxfac(iTracer)
ENDDO
ENDDO
C -- Loop for diffusive fluxes
C -- If done correctly, we can overwrite the flux array in place
C -- with its own divergence
DO k=1,Nr
kdown= MIN(k+1,Nr)
IF (k.EQ.Nr) THEN
downfac = 0. _d 0
ELSE
downfac = 1. _d 0
ENDIF
DO j=1-OLy,sNy+OLy-1
DO i=1-OLx,sNx+OLx-1
C -- Diffusion
tmpfac = -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
& *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
layers_dfx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
& tmpfac * ( layers_dfx(i+1,j,k,iTracer,bi,bj) -
& layers_dfx(i,j,k,iTracer,bi,bj) )
layers_dfy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
& tmpfac * ( layers_dfy(i,j+1,k,iTracer,bi,bj) -
& layers_dfy(i,j,k,iTracer,bi,bj) )
layers_dfr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
& ( layers_dfr(i,j,kdown,iTracer,bi,bj)*downfac -
& layers_dfr(i,j,k,iTracer,bi,bj) )
C -- Advection
layers_afx(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
& tmpfac * ( layers_afx(i+1,j,k,iTracer,bi,bj) -
& layers_afx(i,j,k,iTracer,bi,bj) )
layers_afy(i,j,k,iTracer,bi,bj) = maskInC(i,j,bi,bj) *
& tmpfac * ( layers_afy(i,j+1,k,iTracer,bi,bj) -
& layers_afy(i,j,k,iTracer,bi,bj) )
layers_afr(i,j,k,iTracer,bi,bj) = tmpfac * rkSign *
& ( layers_afr(i,j,kdown,iTracer,bi,bj)*downfac -
& layers_afr(i,j,k,iTracer,bi,bj) )
#ifdef SHORTWAVE_HEATING
IF (iTracer.EQ.1) THEN
swfracb(1)=abs(rF(k))
swfracb(2)=abs(rF(k+1))
CALL SWFRAC(
I 2, minusOne,
U swfracb,
I 1.0, 1, myThid )
C ----- debuggin
C IF ((i.EQ.0).AND.(j.EQ.0)) THEN
C WRITE(msgBuf,'(2A,I3,A,F6.2,A,F6.2)')
C & 'S/R LAYERS_THERMODYNAMICS:',
C & ' k=', k,
C & ' swfracb(1)=', swfracb(1),
C & ' swfracb(2)=', swfracb(2)
C CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
C & SQUEEZE_RIGHT, myThid )
C ENDIF
C --- kdown == kp1
C kp1 = k+1
IF (k.EQ.Nr) THEN
C kp1 = k
swfracb(2)=0. _d 0
ENDIF
layers_sw(i,j,k,iTracer,bi,bj) =
& layers_sw(i,j,k,iTracer,bi,bj)
& -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
& -swfracb(2)*maskC(i,j,kdown,bi,bj))
& *fluxfac(1)
& *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
ENDIF
#endif /* SHORTWAVE_HEATING */
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
C- TFLUX (=total heat flux, match heat-content variations, [W/m2])
C IF ( fluidIsWater .AND.
C & DIAGNOSTICS_IS_ON('TFLUX ',myThid) ) THEN
C DO bj = myByLo(myThid), myByHi(myThid)
C DO bi = myBxLo(myThid), myBxHi(myThid)
C DO j = 1,sNy
C DO i = 1,sNx
C tmp1k(i,j,bi,bj) =
C #ifdef SHORTWAVE_HEATING
C & -Qsw(i,j,bi,bj)+
C #endif
C & (surfaceForcingT(i,j,bi,bj)+surfaceForcingTice(i,j,bi,bj))
C & *HeatCapacity_Cp*rUnit2mass
C ENDDO
C ENDDO
C #ifdef NONLIN_FRSURF
C IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
C & .AND. useRealFreshWaterFlux ) THEN
C DO j=1,sNy
C DO i=1,sNx
C tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
C & + PmEpR(i,j,bi,bj)*theta(i,j,ks,bi,bj)*HeatCapacity_Cp
C ENDDO
C ENDDO
C ENDIF
C #endif /* NONLIN_FRSURF */
C ENDDO
C ENDDO
C CALL DIAGNOSTICS_FILL( tmp1k,'TFLUX ',0,1,0,1,1,myThid )
C ENDIF
C
C C- SFLUX (=total salt flux, match salt-content variations [g/m2/s])
C IF ( fluidIsWater .AND.
C & DIAGNOSTICS_IS_ON('SFLUX ',myThid) ) THEN
C DO bj = myByLo(myThid), myByHi(myThid)
C DO bi = myBxLo(myThid), myBxHi(myThid)
C DO j = 1,sNy
C DO i = 1,sNx
C tmp1k(i,j,bi,bj) =
C & surfaceForcingS(i,j,bi,bj)*rUnit2mass
C ENDDO
C ENDDO
C
C #ifdef NONLIN_FRSURF
C IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
C & .AND. useRealFreshWaterFlux ) THEN
C DO j=1,sNy
C DO i=1,sNx
C tmp1k(i,j,bi,bj) = tmp1k(i,j,bi,bj)
C & + PmEpR(i,j,bi,bj)*salt(i,j,ks,bi,bj)
C ENDDO
C ENDDO
C ENDIF
C #endif /* NONLIN_FRSURF */
C
C ENDDO
C ENDDO
C CALL DIAGNOSTICS_FILL( tmp1k,'SFLUX ',0,1,0,1,1,myThid )
C ENDIF
C Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
C IF ( kLev .EQ. kSurface ) THEN
C DO j=1,sNy
C DO i=1,sNx
C gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
C & +surfaceForcingT(i,j,bi,bj)
C & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
C ENDDO
C ENDDO
C ELSEIF ( kSurface.EQ.-1 ) THEN
C DO j=1,sNy
C DO i=1,sNx
C IF ( kSurfC(i,j,bi,bj).EQ.kLev ) THEN
C gT(i,j,kLev,bi,bj)=gT(i,j,kLev,bi,bj)
C & +surfaceForcingT(i,j,bi,bj)
C & *recip_drF(kLev)*_recip_hFacC(i,j,kLev,bi,bj)
C ENDIF
C ENDDO
C ENDDO
C ENDIF
C-- Divergence of fluxes
C Anelastic: scale vertical fluxes by rhoFac and leave Horizontal fluxes unchanged
C for Stevens OBC: keep only vertical diffusive contribution on boundaries
C DO j=1-OLy,sNy+OLy-1
C DO i=1-OLx,sNx+OLx-1
C gTracer(i,j,k,bi,bj)=gTracer(i,j,k,bi,bj)
C & -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
C & *recip_rA(i,j,bi,bj)*recip_deepFac2C(k)*recip_rhoFacC(k)
C & *( (fZon(i+1,j)-fZon(i,j))*maskInC(i,j,bi,bj)
C & +(fMer(i,j+1)-fMer(i,j))*maskInC(i,j,bi,bj)
C & +(fVerT(i,j,kDown)-fVerT(i,j,kUp))*rkSign
C & -localT(i,j)*( (uTrans(i+1,j)-uTrans(i,j))*advFac
C & +(vTrans(i,j+1)-vTrans(i,j))*advFac
C & +(rTransKp1(i,j)-rTrans(i,j))*rAdvFac
C & )*maskInC(i,j,bi,bj)
C & )
C ENDDO
C ENDDO
#endif /* LAYERS_THERMODYNAMICS */
#endif /* USE_LAYERS */
RETURN
END
C -- end of S/R LAYERS_CALC_RHS