Back to home page

MITgcm
 
 

    

File indexing completed on 2024-10-09 05:10:45 UTC

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

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!