C $Header: /u/gcmpack/MITgcm/pkg/ecco/ecco_phys.F,v 1.18 2017/04/03 23:16:38 ou.wang Exp $
C $Name: $
#include "ECCO_OPTIONS.h"
subroutine ECCO_PHYS( mythid )
c ==================================================================
c SUBROUTINE ecco_phys
c ==================================================================
c
c ==================================================================
c SUBROUTINE ecco_phys
c ==================================================================
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "GRID.h"
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
#ifdef ALLOW_PTRACERS
# include "PTRACERS_SIZE.h"
# include "PTRACERS_FIELDS.h"
#endif
c == routine arguments ==
integer mythid
c == local variables ==
integer bi,bj
integer i,j,k
integer jmin,jmax
integer imin,imax
#ifdef ALLOW_GENCOST_CONTRIBUTION
integer kgen, kgen3d, itr
_RL areavolTile(nSx,nSy), areavolGlob
_RL tmpfld, tmpvol, tmpmsk, tmpmsk2, tmpmskW, tmpmskS
#endif
c- note defined with overlap here, not needed, but more efficient
_RL trVolW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL trVolS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL trHeatW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL trHeatS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL trSaltW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL trSaltS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
#ifdef ATMOSPHERIC_LOADING
_RL sIceLoadFac
#endif
#ifdef ALLOW_PSBAR_STERIC
_RL RHOInSituLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
_RL VOLsumTile(nSx,nSy),RHOsumTile(nSx,nSy)
#endif
c need to include halos for find_rho_2d
iMin = 1-OLx
iMax = sNx+OLx
jMin = 1-OLy
jMax = sNy+OLy
#ifdef ALLOW_PSBAR_STERIC
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
do k = 1,nr
CALL FIND_RHO_2D(
I iMin, iMax, jMin, jMax, k,
I theta(1-OLx,1-OLy,k,bi,bj),
I salt (1-OLx,1-OLy,k,bi,bj),
O RHOInSituLoc(1-OLx,1-OLy,k,bi,bj),
I k, bi, bj, myThid )
enddo
enddo
enddo
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
RHOsumTile(bi,bj)=0. _d 0
VOLsumTile(bi,bj)=0. _d 0
VOLsumGlob=0. _d 0
RHOsumGlob=0. _d 0
do k = 1,nr
do j = 1,sNy
do i = 1,sNx
RHOsumTile(bi,bj)=RHOsumTile(bi,bj)+
& (rhoConst+RHOInSituLoc(i,j,k,bi,bj))*
& hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
VOLsumTile(bi,bj)=VOLsumTile(bi,bj)+
& hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
enddo
enddo
enddo
enddo
enddo
CALL GLOBAL_SUM_TILE_RL( VOLsumTile, VOLsumGlob, myThid )
CALL GLOBAL_SUM_TILE_RL( RHOsumTile, RHOsumGlob, myThid )
RHOsumGlob=RHOsumGlob/VOLsumGlob
if (RHOsumGlob_0.GT.0. _d 0) then
sterGloH=VOLsumGlob_0/globalArea
& *(1. _d 0 - RHOsumGlob/RHOsumGlob_0)
else
sterGloH=0. _d 0
endif
c WRITE(msgBuf,'(A,1PE21.14)') ' sterGloH= ', sterGloH
c CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
c & SQUEEZE_RIGHT, myThid )
#endif
#ifdef ATMOSPHERIC_LOADING
sIceLoadFac=zeroRL
IF ( useRealFreshWaterFlux ) sIceLoadFac=recip_rhoConst
#endif
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
do j = jmin,jmax
do i = imin,imax
m_eta(i,j,bi,bj)=
& etan(i,j,bi,bj)
#ifdef ATMOSPHERIC_LOADING
& +sIceLoad(i,j,bi,bj)*sIceLoadFac
#endif
#ifdef ALLOW_PSBAR_STERIC
& +sterGloH
#endif
enddo
enddo
enddo
enddo
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
do k = 1,nr
do j = 1,sNy
do i = 1,sNx
m_UE(i,j,k,bi,bj)=0. _d 0
m_VN(i,j,k,bi,bj)=0. _d 0
enddo
enddo
enddo
enddo
enddo
CALL ROTATE_UV2EN_RL(
U uVel, vVel, m_UE, m_VN,
I .TRUE., .TRUE., .FALSE., Nr, mythid )
c-- trVol : volume flux --- [m^3/sec] (order of 10^6 = 1 Sv)
c-- trHeat: heat transport --- [Watt] (order of 1.E15 = PW)
c-- trSalt: salt transport --- [kg/sec] (order 1.E9 equiv. 1 Sv in vol.)
c-- convert from [ppt*m^3/sec] via rhoConst/1000.
c-- ( 1ppt = 1000*[mass(salt)]/[mass(seawater)] )
c-- init
call ECCO_ZERO(trVol,Nr,zeroRL,myThid)
call ECCO_ZERO(trHeat,Nr,zeroRL,myThid)
call ECCO_ZERO(trSalt,Nr,zeroRL,myThid)
DO bj=myByLo(myThid),myByHi(myThid)
DO bi=myBxLo(myThid),myBxHi(myThid)
c-- init: if done with overlap, more efficient. But since overwritten
c immediately afterward, init is probably not needed.
do k = 1,nr
do j = 1-OLy,sNy+Oly
do i = 1-OLx,sNx+OLx
trVolW(i,j,k)=0. _d 0
trVolS(i,j,k)=0. _d 0
trHeatW(i,j,k)=0. _d 0
trHeatS(i,j,k)=0. _d 0
trSaltW(i,j,k)=0. _d 0
trSaltS(i,j,k)=0. _d 0
enddo
enddo
enddo
do k = 1,nr
do j = 1,sNy
do i = 1,sNx
trVolW(i,j,k) =
& uVel(i,j,k,bi,bj)*hFacW(i,j,k,bi,bj)
& *dyG(i,j,bi,bj)*drF(k)*msktrVolW(i,j,bi,bj)
& *maskInW(i,j,bi,bj)
trVolS(i,j,k) =
& vVel(i,j,k,bi,bj)*hFacS(i,j,k,bi,bj)
& *dxG(i,j,bi,bj)*drF(k)*msktrVolS(i,j,bi,bj)
& *maskInS(i,j,bi,bj)
trHeatW(i,j,k) = trVolW(i,j,k)
& *(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*halfRL
& *HeatCapacity_Cp*rhoConst
trHeatS(i,j,k) = trVolS(i,j,k)
& *(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*halfRL
& *HeatCapacity_Cp*rhoConst
trSaltW(i,j,k) = trVolW(i,j,k)
& *(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*halfRL
& *rhoConst/1000.
trSaltS(i,j,k) = trVolS(i,j,k)
& *(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*halfRL
& *rhoConst/1000.
c now summing
trVol(i,j,k,bi,bj)=trVolW(i,j,k)+trVolS(i,j,k)
trHeat(i,j,k,bi,bj)=trHeatW(i,j,k)+trHeatS(i,j,k)
trSalt(i,j,k,bi,bj)=trSaltW(i,j,k)+trSaltS(i,j,k)
enddo
enddo
enddo
enddo
enddo
#ifdef ALLOW_GENCOST_CONTRIBUTION
do kgen=1,NGENCOST
itr = gencost_itracer(kgen)
call ECCO_ZERO(gencost_storefld(1-OLx,1-OLy,1,1,kgen),
& 1,zeroRL,myThid)
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
areavolTile(bi,bj)=0. _d 0
enddo
enddo
areavolGlob=0. _d 0
do bj=myByLo(myThid),myByHi(myThid)
do bi=myBxLo(myThid),myBxHi(myThid)
do j = 1,sNy
do i = 1,sNx
c---------
do k = 1,nr
tmpvol=hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
c
tmpmsk=0. _d 0
if (.NOT.gencost_msk_is3d(kgen)) then
tmpmsk=gencost_mskCsurf(i,j,bi,bj,kgen)*
& gencost_mskVertical(k,kgen)
#ifdef ALLOW_GENCOST3D
else
kgen3d=gencost_msk_pointer3d(kgen)
tmpmsk=gencost_mskC(i,j,k,bi,bj,kgen3d)
#endif
endif
c
tmpfld=0. _d 0
tmpmsk2=0. _d 0
if (gencost_barfile(kgen)(1:15).EQ.'m_boxmean_theta') then
tmpfld=theta(i,j,k,bi,bj)
if (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
elseif (gencost_barfile(kgen)(1:14).EQ.'m_boxmean_salt')
& then
tmpfld=salt(i,j,k,bi,bj)
if (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
#ifdef ALLOW_PTRACERS
elseif (gencost_barfile(kgen)(1:17).EQ.'m_boxmean_ptracer')
& then
tmpfld=pTracer(i,j,k,bi,bj,itr)
if (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
#endif
endif
c
gencost_storefld(i,j,bi,bj,kgen) =
& gencost_storefld(i,j,bi,bj,kgen)
& +tmpmsk*tmpfld*tmpvol
areavolTile(bi,bj)=areavolTile(bi,bj)
& +tmpmsk2*eccoVol_0(i,j,k,bi,bj)
c
enddo
c---------
tmpmsk=maskC(i,j,1,bi,bj)*gencost_mskCsurf(i,j,bi,bj,kgen)
tmpfld=0. _d 0
tmpmsk2=0. _d 0
if (gencost_barfile(kgen)(1:13).EQ.'m_boxmean_eta') then
tmpfld=m_eta(i,j,bi,bj)
if (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
endif
c
gencost_storefld(i,j,bi,bj,kgen) =
& gencost_storefld(i,j,bi,bj,kgen)
& +tmpmsk*tmpfld*rA(i,j,bi,bj)
areavolTile(bi,bj)=areavolTile(bi,bj)
& +tmpmsk2*rA(i,j,bi,bj)
c---------
do k = 1,nr
c
tmpmskW=0. _d 0
tmpmskS=0. _d 0
if (.NOT.gencost_msk_is3d(kgen)) then
tmpmskW=gencost_mskWsurf(i,j,bi,bj,kgen)
& *gencost_mskVertical(k,kgen)
tmpmskS=gencost_mskSsurf(i,j,bi,bj,kgen)
& *gencost_mskVertical(k,kgen)
#ifdef ALLOW_GENCOST3D
else
kgen3d=gencost_msk_pointer3d(kgen)
tmpmskW=gencost_mskW(i,j,k,bi,bj,kgen3d)
tmpmskS=gencost_mskS(i,j,k,bi,bj,kgen3d)
#endif
endif
tmpmskW=tmpmskW*hFacW(i,j,k,bi,bj)*dyG(i,j,bi,bj)*drF(k)
tmpmskS=tmpmskS*hFacS(i,j,k,bi,bj)*dxG(i,j,bi,bj)*drF(k)
c
if (gencost_barfile(kgen)(1:13).EQ.'m_horflux_vol') then
gencost_storefld(i,j,bi,bj,kgen) =
& gencost_storefld(i,j,bi,bj,kgen)
& +uVel(i,j,k,bi,bj)*tmpmskW
& +vVel(i,j,k,bi,bj)*tmpmskS
endif
enddo
c---------
enddo
enddo
enddo
enddo
if (gencost_barfile(kgen)(1:9).EQ.'m_boxmean') then
CALL GLOBAL_SUM_TILE_RL( areavolTile, areavolGlob, myThid )
CALL ECCO_DIV( gencost_storefld(1-OLx,1-OLy,1,1,kgen),
& 1, areavolGlob, myThid )
endif
enddo
#endif /* ALLOW_GENCOST_CONTRIBUTION */
return
end