pkg/ecco/ecco_phys.F

2eab8f7357 Gael*0001 #include "ECCO_OPTIONS.h"
7b8b86ab99 Timo*0002 #ifdef ALLOW_SHELFICE
                0003 # include "SHELFICE_OPTIONS.h"
                0004 #endif
35c4fdc74b Emma*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
2eab8f7357 Gael*0008
aa93ca8e85 Ciar*0009 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0010 CBOP
                0011 C     !ROUTINE: ECCO_PHYS
2eab8f7357 Gael*0012
aa93ca8e85 Ciar*0013 C     !INTERFACE:
0cd0083da8 Ou W*0014       SUBROUTINE ECCO_PHYS( myTime, myIter, myThid )
2eab8f7357 Gael*0015
aa93ca8e85 Ciar*0016 C     !DESCRIPTION:
                0017 C     Compute some derived quantities and averages
                0018 C     for GenCost and ECCO cost function.
2eab8f7357 Gael*0019
aa93ca8e85 Ciar*0020 C     !USES:
                0021       IMPLICIT NONE
2eab8f7357 Gael*0022 c     == global variables ==
                0023 #include "EEPARAMS.h"
                0024 #include "SIZE.h"
                0025 #include "PARAMS.h"
aa93ca8e85 Ciar*0026 #include "GRID.h"
2eab8f7357 Gael*0027 #include "DYNVARS.h"
                0028 #include "FFIELDS.h"
49484c0542 Gael*0029 #ifdef ALLOW_ECCO
13d362b8c1 Ou W*0030 # include "ECCO_SIZE.h"
                0031 # include "ECCO.h"
2eab8f7357 Gael*0032 #endif
81e05fa829 Gael*0033 #ifdef ALLOW_PTRACERS
                0034 # include "PTRACERS_SIZE.h"
                0035 # include "PTRACERS_FIELDS.h"
                0036 #endif
7b8b86ab99 Timo*0037 #if (defined ALLOW_GENCOST_CONTRIBUTION) && (defined ALLOW_SHELFICE)
                0038 # include "SHELFICE.h"
                0039 #endif
35c4fdc74b Emma*0040 #ifdef ALLOW_AUTODIFF_TAMC
                0041 # include "tamc.h"
                0042 #endif
2eab8f7357 Gael*0043
aa93ca8e85 Ciar*0044 C     !INPUT PARAMETERS:
0cd0083da8 Ou W*0045 C     myTime    :: Current time in simulation
                0046 C     myIter    :: Current time-step number
                0047 C     myThid    :: my Thread Id number
                0048       _RL     myTime
                0049       INTEGER myIter, myThid
2eab8f7357 Gael*0050
aa93ca8e85 Ciar*0051 C     !LOCAL VARIABLES:
0cd0083da8 Ou W*0052 C     bi, bj    :: tile indices
                0053 C     i, j, k   :: loop indices
aa93ca8e85 Ciar*0054       INTEGER bi, bj
                0055       INTEGER i, j, k
bdae1843b8 Gael*0056 #ifdef ALLOW_GENCOST_CONTRIBUTION
aa93ca8e85 Ciar*0057       INTEGER kgen, kgen3d, itr
17944dd1e8 Gael*0058       _RL areavolTile(nSx,nSy), areavolGlob
                0059       _RL tmpfld, tmpvol, tmpmsk, tmpmsk2, tmpmskW, tmpmskS
aa93ca8e85 Ciar*0060       _RL tmp_sigmsk, tmpsig, tmpsig_lower, tmpsig_upper
bdae1843b8 Gael*0061 #endif
2eab8f7357 Gael*0062
0cd0083da8 Ou W*0063 C- note: defined with overlap here, not needed, but more efficient
fce41e6d01 An T*0064       _RL trVolW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0065       _RL trVolS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0066       _RL trHeatW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0067       _RL trHeatS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0068       _RL trSaltW(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0069       _RL trSaltS(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0070
0cd0083da8 Ou W*0071       _RL tmpfac
                0072       _RL sIceLoadFacLoc
447bdc4b79 Gael*0073 #ifdef ATMOSPHERIC_LOADING
13d362b8c1 Ou W*0074 #ifdef ALLOW_IB_CORR
0cd0083da8 Ou W*0075       _RL ploadbar, AREAsumGlob, PLOADsumGlob
                0076       _RL AREAsumTile(nSx,nSy), PLOADsumTile(nSx,nSy)
                0077       _RL m_eta_ib(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0078       _RL sterht  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
13d362b8c1 Ou W*0079 #endif
447bdc4b79 Gael*0080 #endif
0cd0083da8 Ou W*0081       _RL rhoLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
2eab8f7357 Gael*0082 #ifdef ALLOW_PSBAR_STERIC
0cd0083da8 Ou W*0083       _RL VOLsumTile(nSx,nSy), RHOsumTile(nSx,nSy)
b0b45f2373 Ou W*0084       _RL VOLsumGlob_1, RHOsumGlob_1
0cd0083da8 Ou W*0085 c     CHARACTER*(MAX_LEN_MBUF) msgBuf
2eab8f7357 Gael*0086 #endif
0cd0083da8 Ou W*0087 C     Mload     :: total mass load (kg/m**2)
                0088 c     _RL Mload(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
35c4fdc74b Emma*0089 #if ( defined ALLOW_AUTODIFF_TAMC && defined ALLOW_GENCOST_CONTRIBUTION )
                0090 C     ikey   :: tape key (depends on kgen)
                0091       INTEGER ikey
                0092 #endif
aa93ca8e85 Ciar*0093 CEOP
2eab8f7357 Gael*0094
0cd0083da8 Ou W*0095       tmpfac = recip_rhoConst*recip_gravity
                0096       sIceLoadFacLoc = zeroRL
                0097       IF ( useRealFreshWaterFlux ) sIceLoadFacLoc = recip_rhoConst
49484c0542 Gael*0098
2eab8f7357 Gael*0099       DO bj=myByLo(myThid),myByHi(myThid)
                0100        DO bi=myBxLo(myThid),myBxHi(myThid)
0cd0083da8 Ou W*0101         IF ( myIter .EQ. -1 ) THEN
                0102           DO k = 1,Nr
                0103             CALL FIND_RHO_2D(
aa93ca8e85 Ciar*0104      I                1-OLx, sNx+OLx, 1-OLy, sNy+OLy, k,
2eab8f7357 Gael*0105      I                theta(1-OLx,1-OLy,k,bi,bj),
                0106      I                salt (1-OLx,1-OLy,k,bi,bj),
0cd0083da8 Ou W*0107      O                rhoLoc(1-OLx,1-OLy,k,bi,bj),
2eab8f7357 Gael*0108      I                k, bi, bj, myThid )
0cd0083da8 Ou W*0109           ENDDO
                0110         ELSE
                0111           DO k = 1,Nr
                0112            DO j=1-OLy,sNy+OLy
                0113             DO i=1-OLx,sNx+OLx
                0114               rhoLoc(i,j,k,bi,bj) = rhoInSitu(i,j,k,bi,bj)
                0115             ENDDO
                0116            ENDDO
                0117           ENDDO
                0118         ENDIF
aa93ca8e85 Ciar*0119        ENDDO
                0120       ENDDO
2eab8f7357 Gael*0121
0cd0083da8 Ou W*0122 #ifdef ALLOW_PSBAR_STERIC
2eab8f7357 Gael*0123       DO bj=myByLo(myThid),myByHi(myThid)
                0124        DO bi=myBxLo(myThid),myBxHi(myThid)
0cd0083da8 Ou W*0125         RHOsumTile(bi,bj) = 0. _d 0
                0126         VOLsumTile(bi,bj) = 0. _d 0
aa93ca8e85 Ciar*0127         DO k = 1,Nr
                0128          DO j = 1,sNy
                0129           DO i =  1,sNx
0cd0083da8 Ou W*0130            RHOsumTile(bi,bj) = RHOsumTile(bi,bj)
                0131      &          + ( rhoConst + rhoLoc(i,j,k,bi,bj) )
                0132      &           *hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
                0133            VOLsumTile(bi,bj) = VOLsumTile(bi,bj)
                0134      &          + hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
aa93ca8e85 Ciar*0135           ENDDO
                0136          ENDDO
                0137         ENDDO
                0138        ENDDO
                0139       ENDDO
b0b45f2373 Ou W*0140       CALL GLOBAL_SUM_TILE_RL( VOLsumTile, VOLsumGlob_1, myThid )
                0141       CALL GLOBAL_SUM_TILE_RL( RHOsumTile, RHOsumGlob_1, myThid )
                0142
3a516654c6 Jean*0143 # ifndef ALLOW_AUTODIFF
b0b45f2373 Ou W*0144       _BEGIN_MASTER(myThid)
3a516654c6 Jean*0145 # endif
b0b45f2373 Ou W*0146       VOLsumGlob = VOLsumGlob_1
                0147       RHOsumGlob = RHOsumGlob_1/VOLsumGlob
2eab8f7357 Gael*0148
aa93ca8e85 Ciar*0149       IF (RHOsumGlob_0.GT.0. _d 0) THEN
0cd0083da8 Ou W*0150         sterGloH = VOLsumGlob_0/globalArea
49484c0542 Gael*0151      &        *(1. _d 0 - RHOsumGlob/RHOsumGlob_0)
aa93ca8e85 Ciar*0152       ELSE
0cd0083da8 Ou W*0153         sterGloH = 0. _d 0
aa93ca8e85 Ciar*0154       ENDIF
49484c0542 Gael*0155
                0156 c     WRITE(msgBuf,'(A,1PE21.14)') ' sterGloH= ', sterGloH
0cd0083da8 Ou W*0157 c     CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
49484c0542 Gael*0158 c    &                       SQUEEZE_RIGHT, myThid )
3a516654c6 Jean*0159 # ifndef ALLOW_AUTODIFF
b0b45f2373 Ou W*0160       _END_MASTER(myThid)
                0161       _BARRIER
3a516654c6 Jean*0162 # endif
49484c0542 Gael*0163
0cd0083da8 Ou W*0164 #endif /* ALLOW_PSBAR_STERIC */
2eab8f7357 Gael*0165
447bdc4b79 Gael*0166 #ifdef ATMOSPHERIC_LOADING
13d362b8c1 Ou W*0167 #ifdef ALLOW_IB_CORR
                0168       DO bj=myByLo(myThid),myByHi(myThid)
                0169        DO bi=myBxLo(myThid),myBxHi(myThid)
0cd0083da8 Ou W*0170         PLOADsumTile(bi,bj) = 0. _d 0
                0171         AREAsumTile(bi,bj) = 0. _d 0
aa93ca8e85 Ciar*0172         DO j = 1,sNy
                0173          DO i =  1,sNx
0cd0083da8 Ou W*0174           PLOADsumTile(bi,bj) = PLOADsumTile(bi,bj)
                0175      &          + pLoad(i,j,bi,bj)
                0176      &           *maskC(i,j,1,bi,bj)*rA(i,j,bi,bj)
                0177           AREAsumTile(bi,bj) = AREAsumTile(bi,bj)
                0178      &          + maskC(i,j,1,bi,bj)*rA(i,j,bi,bj)
aa93ca8e85 Ciar*0179          ENDDO
                0180         ENDDO
                0181        ENDDO
                0182       ENDDO
13d362b8c1 Ou W*0183       CALL GLOBAL_SUM_TILE_RL( AREAsumTile, AREAsumGlob, myThid )
                0184       CALL GLOBAL_SUM_TILE_RL( PLOADsumTile, PLOADsumGlob, myThid )
0cd0083da8 Ou W*0185       ploadbar = PLOADsumGlob/AREAsumGlob
                0186 #endif /* ALLOW_IB_CORR */
                0187 #endif /* ATMOSPHERIC_LOADING */
2eab8f7357 Gael*0188
49484c0542 Gael*0189       DO bj=myByLo(myThid),myByHi(myThid)
                0190        DO bi=myBxLo(myThid),myBxHi(myThid)
aa93ca8e85 Ciar*0191         DO j = 1-OLy, sNy+OLy
                0192          DO i = 1-OLx, sNx+OLx
0cd0083da8 Ou W*0193
                0194 C calculate total sea level including inverse barometer (IB) effect if
13d362b8c1 Ou W*0195 C  there is air pressure forcing
0cd0083da8 Ou W*0196            m_eta(i,j,bi,bj) =
                0197      &           ( etaN(i,j,bi,bj)
                0198      &           + sIceLoad(i,j,bi,bj)*sIceLoadFacLoc
49484c0542 Gael*0199 #ifdef ALLOW_PSBAR_STERIC
0cd0083da8 Ou W*0200 C apply Greatbatch correction
                0201      &           + sterGloH
                0202 #endif /* ALLOW_PSBAR_STERIC */
                0203      &           ) * maskC(i,j,1,bi,bj)
                0204
                0205 C Model equivalent of ocean bottom pressure gauge data (in m^2/s^2)
                0206 C  that are NOT corrected for global ocean mean atmospheric pressure variations
                0207 C  = Ocean mass + sea-ice & snow load + air-pressure load + Greatbatch corr.
                0208 C  (all terms on RHS are converted to m^2/s^2). It is
                0209 C  essentially Mload (as in pkg/sbo/sbo_calc.F) plus air-pressure load
                0210            m_bp(i,j,bi,bj) =
                0211      &           ( etaN(i,j,bi,bj)
                0212 c    &           + Ro_surf(i,j,bi.bj)
                0213      &           - R_low(i,j,bi,bj)
13d362b8c1 Ou W*0214 #ifdef ALLOW_PSBAR_STERIC
                0215 C add back the correction due to the global mean steric ssh change,
0cd0083da8 Ou W*0216 C     i.e. sterGloH computed above (units converted from m to m2/s2)
                0217      &           + sterGloH
                0218 #endif /* ALLOW_PSBAR_STERIC */
                0219      &           ) * gravity
                0220 C sIceLoad in kg/m^2
                0221      &         + sIceLoad(i,j,bi,bj) * gravity * sIceLoadFacLoc
                0222 C pLoad in N/m^2
                0223      &         + pLoad(i,j,bi,bj) * recip_rhoConst
                0224          ENDDO
                0225         ENDDO
                0226 C integrate rho_anomaly through water column
                0227         DO k = 1, Nr
                0228          DO j = 1-OLy, sNy+OLy
                0229           DO i = 1-OLx, sNx+OLx
                0230            m_bp(i,j,bi,bj) = m_bp(i,j,bi,bj)
                0231      &         + rhoLoc(i,j,k,bi,bj)*drF(k)*hFacC(i,j,k,bi,bj)
                0232      &                              * gravity * recip_rhoConst
                0233           ENDDO
                0234          ENDDO
                0235         ENDDO
                0236         DO j = 1-OLy, sNy+OLy
                0237          DO i = 1-OLx, sNx+OLx
                0238            m_bp(i,j,bi,bj) = m_bp(i,j,bi,bj) * maskC(i,j,1,bi,bj)
                0239          ENDDO
                0240         ENDDO
13d362b8c1 Ou W*0241
                0242 #ifdef ATMOSPHERIC_LOADING
                0243 #ifdef ALLOW_IB_CORR
0cd0083da8 Ou W*0244         DO j = 1-OLy, sNy+OLy
                0245          DO i = 1-OLx, sNx+OLx
13d362b8c1 Ou W*0246 C calculate IB correction m_eta_ib (in m)
0cd0083da8 Ou W*0247            m_eta_ib(i,j,bi,bj) =
                0248      &           ( ploadbar - pLoad(i,j,bi,bj) )*tmpfac
13d362b8c1 Ou W*0249      &           * maskC(i,j,1,bi,bj)
                0250 C calculte dynamic sea level for comparison with altimetry data (in m)
0cd0083da8 Ou W*0251            m_eta_dyn(i,j,bi,bj) =
                0252      &           ( m_eta(i,j,bi,bj) - m_eta_ib(i,j,bi,bj) )
13d362b8c1 Ou W*0253      &           * maskC(i,j,1,bi,bj)
                0254
                0255 C calculate GRACE-equvivalent ocean bottom pressure (in m2/s2)
0cd0083da8 Ou W*0256 C  by removing global ocean mean atmospheric pressure variations
                0257            m_bp_nopabar(i,j,bi,bj) =
                0258      &           ( m_bp(i,j,bi,bj)
                0259      &           - ploadbar * recip_rhoConst
                0260      &           ) * maskC(i,j,1,bi,bj)
                0261 C calculate steric height
                0262 C (in m; = m_eta_dyn - (m_bp_nopabr * recip_gravity + R_low))
                0263 C R_low (<0) is Depth in m.
                0264            sterht(i,j,bi,bj) = m_eta_dyn(i,j,bi,bj)
                0265      &          - ( m_bp_nopabar(i,j,bi,bj) * recip_gravity
                0266 c    &            - Ro_surf(i,j,bi,bj)
                0267      &            + R_low(i,j,bi,bj)  )
aa93ca8e85 Ciar*0268          ENDDO
                0269         ENDDO
0cd0083da8 Ou W*0270 #endif /* ALLOW_IB_CORR */
                0271 #endif /* ATMOSPHERIC_LOADING */
aa93ca8e85 Ciar*0272        ENDDO
                0273       ENDDO
49484c0542 Gael*0274
0cd0083da8 Ou W*0275 #ifdef ALLOW_DIAGNOSTICS
                0276       IF ( useDiagnostics .AND. myIter.GE.0 ) THEN
                0277         CALL DIAGNOSTICS_FILL( m_eta, 'SSHNOIBC', 0,1, 0,1,1, myThid )
                0278         CALL DIAGNOSTICS_SCALE_FILL( m_bp, recip_gravity, 1,
                0279      &                         'OBPGMAP ', 0,1, 0,1,1, myThid )
                0280 #ifdef ATMOSPHERIC_LOADING
                0281 #ifdef ALLOW_IB_CORR
                0282         CALL DIAGNOSTICS_FILL( m_eta_ib,
                0283      &                         'SSHIBC  ', 0,1, 0,1,1, myThid )
                0284         CALL DIAGNOSTICS_FILL( m_eta_dyn,
                0285      &                         'SSH     ', 0,1, 0,1,1, myThid )
                0286         CALL DIAGNOSTICS_FILL( sterht,
                0287      &                         'STERICHT', 0,1, 0,1,1, myThid )
                0288         CALL DIAGNOSTICS_SCALE_FILL( m_bp_nopabar, recip_gravity, 1,
                0289      &                         'OBP     ', 0,1, 0,1,1, myThid )
                0290 #endif /* ALLOW_IB_CORR */
                0291 #endif /* ATMOSPHERIC_LOADING */
                0292       ENDIF
                0293 #endif /* ALLOW_DIAGNOSTICS */
                0294
556841ad42 Gael*0295       DO bj=myByLo(myThid),myByHi(myThid)
                0296        DO bi=myBxLo(myThid),myBxHi(myThid)
aa93ca8e85 Ciar*0297         DO k = 1,Nr
0cd0083da8 Ou W*0298          DO j = 1-OLy,sNy+OLy
                0299           DO i = 1-OLx,sNx+OLx
                0300            m_UE(i,j,k,bi,bj) = 0. _d 0
                0301            m_VN(i,j,k,bi,bj) = 0. _d 0
aa93ca8e85 Ciar*0302           ENDDO
                0303          ENDDO
                0304         ENDDO
                0305        ENDDO
                0306       ENDDO
556841ad42 Gael*0307
                0308       CALL ROTATE_UV2EN_RL(
                0309      U          uVel, vVel, m_UE, m_VN,
13d362b8c1 Ou W*0310      I          .TRUE., .TRUE., .FALSE., Nr, myThid )
0761692d75 An T*0311
                0312 c--   trVol : volume flux    --- [m^3/sec] (order of 10^6 = 1 Sv)
                0313 c--   trHeat: heat transport --- [Watt] (order of 1.E15 = PW)
                0314 c--   trSalt: salt transport --- [kg/sec] (order 1.E9 equiv. 1 Sv in vol.)
                0315 c--       convert from [ppt*m^3/sec] via rhoConst/1000.
                0316 c--       ( 1ppt = 1000*[mass(salt)]/[mass(seawater)] )
                0317
                0318 c-- init
0cd0083da8 Ou W*0319       CALL ECCO_ZERO( trVol,  Nr, zeroRL, myThid )
                0320       CALL ECCO_ZERO( trHeat, Nr, zeroRL, myThid )
                0321       CALL ECCO_ZERO( trSalt, Nr, zeroRL, myThid )
0761692d75 An T*0322
bdae1843b8 Gael*0323 #ifdef ALLOW_GENCOST_CONTRIBUTION
bcdcbe969d Gael*0324
7b8b86ab99 Timo*0325 cts ---
                0326 c First: Fill the following SCALAR masks & weights for each (i,j,k,bi,bj) grid cell
                0327 c   tmpvol - 3D cell volume
                0328 c   tmpmsk - mask for the gencost_barfile field (e.g. theta)
                0329 c            Either: expand from 2D mask gencost_mskCsurf across nonzero
                0330 c            entries of gencost_mskVertical (Nr x NGENCOST array)
                0331 c            or
                0332 c            copy from 3D mask gencost_mskC
                0333 cts ---
aa93ca8e85 Ciar*0334       DO kgen=1,NGENCOST
bdae1843b8 Gael*0335
aa93ca8e85 Ciar*0336        itr = gencost_itracer(kgen)
81e05fa829 Gael*0337
0cd0083da8 Ou W*0338        CALL ECCO_ZERO( gencost_storefld(1-OLx,1-OLy,1,1,kgen),
                0339      &                 1, zeroRL, myThid )
bdae1843b8 Gael*0340
aa93ca8e85 Ciar*0341        DO bj=myByLo(myThid),myByHi(myThid)
                0342         DO bi=myBxLo(myThid),myBxHi(myThid)
bcdcbe969d Gael*0343          areavolTile(bi,bj)=0. _d 0
aa93ca8e85 Ciar*0344         ENDDO
                0345        ENDDO
                0346        areavolGlob=0. _d 0
                0347
                0348        DO bj=myByLo(myThid),myByHi(myThid)
                0349         DO bi=myBxLo(myThid),myBxHi(myThid)
                0350          DO j = 1,sNy
                0351           DO i =  1,sNx
17944dd1e8 Gael*0352 c---------
aa93ca8e85 Ciar*0353            DO k = 1,Nr
17944dd1e8 Gael*0354             tmpvol=hFacC(i,j,k,bi,bj)*drF(k)*rA(i,j,bi,bj)
aa93ca8e85 Ciar*0355
17944dd1e8 Gael*0356             tmpmsk=0. _d 0
aa93ca8e85 Ciar*0357             IF (.NOT.gencost_msk_is3d(kgen)) THEN
                0358              tmpmsk=gencost_mskCsurf(i,j,bi,bj,kgen)*
                0359      &              gencost_mskVertical(k,kgen)
17944dd1e8 Gael*0360 #ifdef ALLOW_GENCOST3D
aa93ca8e85 Ciar*0361             ELSE
                0362              kgen3d=gencost_msk_pointer3d(kgen)
                0363              tmpmsk=gencost_mskC(i,j,k,bi,bj,kgen3d)
0cd0083da8 Ou W*0364 #endif /* ALLOW_GENCOST3D */
aa93ca8e85 Ciar*0365             ENDIF
                0366
                0367 C ---- If density mask is enabled, use it here ----
                0368             IF ( maskC(i,j,k,bi,bj).EQ.oneRS .AND.
                0369      &           gencost_useDensityMask(kgen) ) THEN
                0370 C            - first, calculate the scalar density
                0371              CALL FIND_RHO_SCALAR(
                0372      I              theta(i,j,k,bi,bj),
                0373      I              salt(i,j,k,bi,bj),
                0374      I              gencost_refPressure(kgen),
                0375      O              tmpsig,
                0376      I              myThid )
                0377 C            - subtract 1000 to get sigma
                0378              tmpsig = tmpsig - 1000. _d 0
                0379 C            - now, tmpmsk is sigmoid times this value
                0380              tmpsig_lower = 0.5 + 0.5*tanh(gencost_tanhScale(kgen)
                0381      &           *(tmpsig-gencost_sigmaLow(kgen)))
                0382              tmpsig_upper = 0.5 - 0.5*tanh(gencost_tanhScale(kgen)
                0383      &           *(tmpsig-gencost_sigmaHigh(kgen)))
                0384 C             - update mask value based on the sigmoid function
                0385              tmp_sigmsk = tmpsig_lower*tmpsig_upper
                0386              tmpmsk = tmpmsk*tmp_sigmsk
                0387             ENDIF
                0388 C ---- end of density mask (but tmpmsk is used below)
                0389
7b8b86ab99 Timo*0390 cts ---
                0391 c Now: at each (i,j,k,bi,bj) fill the SCALAR variables
                0392 c   tmpfld - from 3D field theta, salt, ptracer
                0393 c            or
                0394 c            from 2D field with eta, shelfice
                0395 c
                0396 c   tmpmsk2 - 1 or 0 weighting for areavolTile
                0397 cts ---
17944dd1e8 Gael*0398             tmpfld=0. _d 0
                0399             tmpmsk2=0. _d 0
aa93ca8e85 Ciar*0400             IF (gencost_barfile(kgen)(1:15).EQ.'m_boxmean_theta') THEN
                0401              tmpfld=theta(i,j,k,bi,bj)
                0402              IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
                0403             ELSEIF (gencost_barfile(kgen)(1:14).EQ.'m_boxmean_salt')
                0404      &        THEN
                0405              tmpfld=salt(i,j,k,bi,bj)
                0406              IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
38d78826c7 Ciar*0407             ELSEIF (gencost_barfile(kgen)(1:13).EQ.'m_boxmean_vol') THEN
                0408              tmpfld=1. _d 0
                0409              IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
81e05fa829 Gael*0410 #ifdef ALLOW_PTRACERS
aa93ca8e85 Ciar*0411             ELSEIF (gencost_barfile(kgen)(1:17).EQ.'m_boxmean_ptracer')
                0412      &        THEN
                0413              tmpfld=pTracer(i,j,k,bi,bj,itr)
                0414              IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
0cd0083da8 Ou W*0415 #endif /* ALLOW_PTRACERS */
aa93ca8e85 Ciar*0416             ENDIF
                0417
7b8b86ab99 Timo*0418 cts ---
                0419 c Fill 3D field
                0420 c       gencost_store - masked field of interest * grid cell volume
                0421 c                       note: this accumulates along z dim
                0422 c
                0423 c Fill tile field (1 val per tile)
                0424 c       areavolTile - volume of each tile, this gets summed to a global
                0425 c                     value
                0426 cts ---
17944dd1e8 Gael*0427             gencost_storefld(i,j,bi,bj,kgen) =
                0428      &          gencost_storefld(i,j,bi,bj,kgen)
                0429      &          +tmpmsk*tmpfld*tmpvol
                0430             areavolTile(bi,bj)=areavolTile(bi,bj)
35c4fdc74b Emma*0431 #ifdef ECCO_VARIABLE_AREAVOLGLOB
                0432      &          +tmpmsk2*tmpvol
                0433 #else
17944dd1e8 Gael*0434      &          +tmpmsk2*eccoVol_0(i,j,k,bi,bj)
35c4fdc74b Emma*0435 #endif
7b8b86ab99 Timo*0436
aa93ca8e85 Ciar*0437            ENDDO ! Ends do k=1,Nr
                0438
                0439            tmpmsk  = 0. _d 0
                0440            tmpfld  = 0. _d 0
                0441            tmpmsk2 = 0. _d 0
                0442            IF (gencost_barfile(kgen)(1:13).EQ.'m_boxmean_eta') THEN
7b8b86ab99 Timo*0443             tmpmsk=maskC(i,j,1,bi,bj)*gencost_mskCsurf(i,j,bi,bj,kgen)
13d362b8c1 Ou W*0444             tmpfld = m_eta(i,j,bi,bj)
                0445 #if (defined ATMOSPHERIC_LOADING) && (defined ALLOW_IB_CORR)
aa93ca8e85 Ciar*0446             IF (gencost_barfile(kgen)(1:17).EQ.'m_boxmean_eta_dyn') THEN
                0447              tmpfld = m_eta_dyn(i,j,bi,bj)
                0448             ENDIF
13d362b8c1 Ou W*0449 #endif
aa93ca8e85 Ciar*0450             IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
                0451            ENDIF
7b8b86ab99 Timo*0452 #ifdef ALLOW_SHELFICE
                0453 cts ---
                0454 c Shelfice:
                0455 c   Simply accumulate shelfice FWF or HF into tmpfld here
                0456 c   This will fill gencost_storefld with this value *rA
                0457 c   For FreshWaterFlux
                0458 c           gencost_storefld = shelficefreshwaterflux / rho * rA
                0459 c                            = [kg/m^2/s] / [kg/m^3] * [m^2]
                0460 c                            = [m^3/s]
17944dd1e8 Gael*0461 c
7b8b86ab99 Timo*0462 c   For heatflux
                0463 c           gencost_storefld = shelficeheatflux * rA
                0464 c                            = [W/m^2] *[m^2]
                0465 c                            = [W]
                0466 cts ---
aa93ca8e85 Ciar*0467            IF((gencost_barfile(kgen)(1:16).EQ.'m_boxmean_shifwf').OR.
                0468      &        (gencost_barfile(kgen)(1:16).EQ.'m_boxmean_shihtf')) THEN
7b8b86ab99 Timo*0469
                0470             tmpmsk=maskSHI(i,j,1,bi,bj)*
                0471      &             gencost_mskCsurf(i,j,bi,bj,kgen)
                0472
aa93ca8e85 Ciar*0473             IF (gencost_barfile(kgen)(11:16).EQ.'shifwf') THEN
                0474              tmpfld=shelficeFreshWaterFlux(i,j,bi,bj) / rhoConstFresh
                0475             ELSEIF (gencost_barfile(kgen)(11:16).EQ.'shihtf') THEN
                0476              tmpfld=shelficeHeatFlux(i,j,bi,bj)
                0477             ENDIF
                0478             IF (tmpmsk.NE.0. _d 0) tmpmsk2=1. _d 0
                0479            ENDIF
7b8b86ab99 Timo*0480 #endif /* ALLOW_SHELFICE */
                0481
                0482 cts ---
                0483 c Fill 2D field
                0484 c   gencost_store - masked field of interest * rA
                0485 c
                0486 c Fill tile field (1 val per tile)
                0487 c       areavolTile - total rA on each tile for mskC != 0
                0488 cts ---
aa93ca8e85 Ciar*0489            gencost_storefld(i,j,bi,bj,kgen) =
c17b89ca05 Gael*0490      &        gencost_storefld(i,j,bi,bj,kgen)
17944dd1e8 Gael*0491      &        +tmpmsk*tmpfld*rA(i,j,bi,bj)
aa93ca8e85 Ciar*0492            areavolTile(bi,bj)=areavolTile(bi,bj)
17944dd1e8 Gael*0493      &        +tmpmsk2*rA(i,j,bi,bj)
                0494 c---------
aa93ca8e85 Ciar*0495            DO k = 1,Nr
                0496
17944dd1e8 Gael*0497             tmpmskW=0. _d 0
                0498             tmpmskS=0. _d 0
aa93ca8e85 Ciar*0499             IF (.NOT.gencost_msk_is3d(kgen)) THEN
17944dd1e8 Gael*0500               tmpmskW=gencost_mskWsurf(i,j,bi,bj,kgen)
                0501      &          *gencost_mskVertical(k,kgen)
                0502               tmpmskS=gencost_mskSsurf(i,j,bi,bj,kgen)
                0503      &          *gencost_mskVertical(k,kgen)
                0504 #ifdef ALLOW_GENCOST3D
aa93ca8e85 Ciar*0505             ELSE
17944dd1e8 Gael*0506               kgen3d=gencost_msk_pointer3d(kgen)
                0507               tmpmskW=gencost_mskW(i,j,k,bi,bj,kgen3d)
                0508               tmpmskS=gencost_mskS(i,j,k,bi,bj,kgen3d)
0cd0083da8 Ou W*0509 #endif /* ALLOW_GENCOST3D */
aa93ca8e85 Ciar*0510             ENDIF
17944dd1e8 Gael*0511             tmpmskW=tmpmskW*hFacW(i,j,k,bi,bj)*dyG(i,j,bi,bj)*drF(k)
                0512             tmpmskS=tmpmskS*hFacS(i,j,k,bi,bj)*dxG(i,j,bi,bj)*drF(k)
aa93ca8e85 Ciar*0513
                0514             IF (gencost_barfile(kgen)(1:13).EQ.'m_horflux_vol') THEN
17944dd1e8 Gael*0515               gencost_storefld(i,j,bi,bj,kgen) =
                0516      &          gencost_storefld(i,j,bi,bj,kgen)
                0517      &          +uVel(i,j,k,bi,bj)*tmpmskW
                0518      &          +vVel(i,j,k,bi,bj)*tmpmskS
736e27304c Timo*0519
                0520             ! Only compute tr[Vol,Heat,Salt] if necessary, use
7b8b86ab99 Timo*0521             ! gencost_mask[W/S] rather than old msktrVol
aa93ca8e85 Ciar*0522             ELSEIF ( gencost_barfile(kgen)(1:7).EQ.'m_trVol' .OR.
                0523      &               gencost_barfile(kgen)(1:8).EQ.'m_trHeat'.OR.
                0524      &               gencost_barfile(kgen)(1:8).EQ.'m_trSalt'    ) THEN
736e27304c Timo*0525
aa93ca8e85 Ciar*0526              trVolW(i,j,k) =
736e27304c Timo*0527      &                 uVel(i,j,k,bi,bj)*tmpmskW
                0528      &                *maskInW(i,j,bi,bj)
aa93ca8e85 Ciar*0529              trVolS(i,j,k) =
736e27304c Timo*0530      &                 vVel(i,j,k,bi,bj)*tmpmskS
                0531      &                *maskInS(i,j,bi,bj)
                0532
aa93ca8e85 Ciar*0533              trHeatW(i,j,k) = trVolW(i,j,k)
736e27304c Timo*0534      &                *(theta(i,j,k,bi,bj)+theta(i-1,j,k,bi,bj))*halfRL
                0535      &                *HeatCapacity_Cp*rhoConst
aa93ca8e85 Ciar*0536              trHeatS(i,j,k) = trVolS(i,j,k)
736e27304c Timo*0537      &                *(theta(i,j,k,bi,bj)+theta(i,j-1,k,bi,bj))*halfRL
                0538      &                *HeatCapacity_Cp*rhoConst
                0539
aa93ca8e85 Ciar*0540              trSaltW(i,j,k) = trVolW(i,j,k)
736e27304c Timo*0541      &                *(salt(i,j,k,bi,bj)+salt(i-1,j,k,bi,bj))*halfRL
                0542      &                *rhoConst/1000.
aa93ca8e85 Ciar*0543              trSaltS(i,j,k) = trVolS(i,j,k)
736e27304c Timo*0544      &                *(salt(i,j,k,bi,bj)+salt(i,j-1,k,bi,bj))*halfRL
                0545      &                *rhoConst/1000.
                0546 c now summing
aa93ca8e85 Ciar*0547              trVol(i,j,k,bi,bj)=trVolW(i,j,k)+trVolS(i,j,k)
                0548              trHeat(i,j,k,bi,bj)=trHeatW(i,j,k)+trHeatS(i,j,k)
                0549              trSalt(i,j,k,bi,bj)=trSaltW(i,j,k)+trSaltS(i,j,k)
736e27304c Timo*0550
aa93ca8e85 Ciar*0551             ENDIF
                0552 C     end k-loop
                0553            ENDDO
17944dd1e8 Gael*0554 c---------
aa93ca8e85 Ciar*0555           ENDDO
                0556          ENDDO
                0557         ENDDO
                0558        ENDDO
bdae1843b8 Gael*0559
7b8b86ab99 Timo*0560 cts ---
                0561 c Divide all values in gencost_storefld by
                0562 c   areavolGlob: scalar representing global volume of
                0563 c                quantity of interest.
                0564 c
11c3150c71 Mart*0565 c Note: for shelfice, do not take this final average to make
7b8b86ab99 Timo*0566 c       comparable to shelfice_cost_final.
                0567 cts ---
38d78826c7 Ciar*0568        IF ( gencost_barfile(kgen)(1:9).EQ.'m_boxmean' .AND.
                0569      &      gencost_barfile(kgen)(11:13).NE.'shi' .AND.
                0570      &      gencost_barfile(kgen)(11:13).NE.'vol' ) THEN
17944dd1e8 Gael*0571         CALL GLOBAL_SUM_TILE_RL( areavolTile, areavolGlob, myThid )
35c4fdc74b Emma*0572 #ifdef ALLOW_AUTODIFF_TAMC
                0573 C     avoid recomputing areavolGlob multiple times each time this
                0574 C     routine is called
                0575         ikey = kgen + (ikey_dynamics-1)*NGENCOST
                0576 CADJ STORE areavolGlob = comlev1_ngencost, key = ikey
                0577 #endif
0cd0083da8 Ou W*0578         CALL ECCO_DIV( gencost_storefld(1-OLx,1-OLy,1,1,kgen),
11c3150c71 Mart*0579      &                 areavolGlob, 1, 1, myThid )
aa93ca8e85 Ciar*0580        ENDIF
bcdcbe969d Gael*0581
aa93ca8e85 Ciar*0582 C     end kgen-loop
                0583       ENDDO
bdae1843b8 Gael*0584
                0585 #endif /* ALLOW_GENCOST_CONTRIBUTION */
                0586
aa93ca8e85 Ciar*0587       RETURN
                0588       END