model/src/external_forcing_surf.F

6d54cf9ca1 Ed H*0001 #include "PACKAGES_CONFIG.h"
77af23a186 Patr*0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
1f89baba18 Patr*0006 #ifdef ALLOW_SALT_PLUME
                0007 # include "SALT_PLUME_OPTIONS.h"
                0008 #endif
cb5aad258c Jean*0009 #undef CHECK_OVERLAP_FORCING
21b8df1d84 Jean*0010
9366854e02 Chri*0011 CBOP
                0012 C     !ROUTINE: EXTERNAL_FORCING_SURF
                0013 C     !INTERFACE:
21b8df1d84 Jean*0014       SUBROUTINE EXTERNAL_FORCING_SURF(
cb5aad258c Jean*0015      I             iMin, iMax, jMin, jMax,
c61ca13fc6 Dimi*0016      I             myTime, myIter, myThid )
9366854e02 Chri*0017 C     !DESCRIPTION: \bv
                0018 C     *==========================================================*
21b8df1d84 Jean*0019 C     | SUBROUTINE EXTERNAL_FORCING_SURF
                0020 C     | o Determines forcing terms based on external fields
                0021 C     |   relaxation terms etc.
9366854e02 Chri*0022 C     *==========================================================*
                0023 C     \ev
                0024
                0025 C     !USES:
77af23a186 Patr*0026       IMPLICIT NONE
                0027 C     === Global variables ===
                0028 #include "SIZE.h"
                0029 #include "EEPARAMS.h"
                0030 #include "PARAMS.h"
                0031 #include "FFIELDS.h"
                0032 #include "DYNVARS.h"
                0033 #include "GRID.h"
fa41ecc867 Jean*0034 #include "SURFACE.h"
7c50f07931 Mart*0035 #ifdef ALLOW_AUTODIFF_TAMC
62d9359ab4 Patr*0036 # include "tamc.h"
                0037 #endif
e4775240e5 Dimi*0038
9366854e02 Chri*0039 C     !INPUT/OUTPUT PARAMETERS:
77af23a186 Patr*0040 C     === Routine arguments ===
fab93eca08 Jean*0041 C     iMin,iMax, jMin,jMax :: Range of points for calculation
                0042 C     myTime :: Current time in simulation
                0043 C     myIter :: Current iteration number in simulation
9366854e02 Chri*0044 C     myThid :: Thread no. that called this routine.
51701379d8 Ed H*0045       INTEGER iMin, iMax
                0046       INTEGER jMin, jMax
21b8df1d84 Jean*0047       _RL myTime
                0048       INTEGER myIter
                0049       INTEGER myThid
e508fdf6c2 Patr*0050
9366854e02 Chri*0051 C     !LOCAL VARIABLES:
                0052 C     === Local variables ===
cb5aad258c Jean*0053 C     bi,bj  :: tile indices
21b8df1d84 Jean*0054 C     i,j    :: loop indices
                0055 C     ks     :: index of surface interface layer
cb5aad258c Jean*0056       INTEGER bi,bj
51701379d8 Ed H*0057       INTEGER i,j
fab93eca08 Jean*0058       INTEGER ks
7c50f07931 Mart*0059 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0060       INTEGER tkey
7c50f07931 Mart*0061 #endif
faf82d94de Patr*0062       _RL recip_Cp
cb5aad258c Jean*0063 #ifdef ALLOW_BALANCE_FLUXES
                0064       _RS tmpVar(1)
                0065 #endif
                0066 #ifdef CHECK_OVERLAP_FORCING
                0067       _RS fixVal
                0068 #endif
21b8df1d84 Jean*0069 CEOP
9366854e02 Chri*0070
9669509dca Jean*0071       IF ( usingPCoords ) THEN
21b8df1d84 Jean*0072        ks        = Nr
042348d828 Jean*0073       ELSE
fab93eca08 Jean*0074        ks        = 1
042348d828 Jean*0075       ENDIF
faf82d94de Patr*0076       recip_Cp = 1. _d 0 / HeatCapacity_Cp
042348d828 Jean*0077
cb5aad258c Jean*0078 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
62d9359ab4 Patr*0079
cb5aad258c Jean*0080 C--   Apply adjustment (balancing forcing) and exchanges
                0081 C     to oceanic surface forcing
62d9359ab4 Patr*0082
cb5aad258c Jean*0083 #ifdef ALLOW_BALANCE_FLUXES
                0084 C     balance fluxes
7e00d7e8f9 Jean*0085 # ifdef ALLOW_AUTODIFF
cb5aad258c Jean*0086       tmpVar(1) = oneRS
7448700841 Mart*0087 #  ifdef ALLOW_AUTODIFF_TAMC
                0088 CADJ INCOMPLETE tmpVar
                0089 #  endif
7e00d7e8f9 Jean*0090 # endif
                0091       IF ( selectBalanceEmPmR.GE.1 .AND.
                0092      &     (.NOT.useSeaice .OR. useThSIce) ) THEN
                0093        IF ( selectBalanceEmPmR.EQ.1 ) THEN
                0094         tmpVar(1) = oneRS
                0095         CALL REMOVE_MEAN_RS( 1, EmPmR, maskInC, maskInC, rA,
                0096      &                       tmpVar, 'EmPmR', myTime, myIter, myThid )
                0097        ELSEIF ( selectBalanceEmPmR.EQ.2 ) THEN
                0098         tmpVar(1) = -oneRS
                0099         CALL REMOVE_MEAN_RS( 1, EmPmR, weight2BalanceFlx, maskInC, rA,
                0100      &                       tmpVar, 'EmPmR', myTime, myIter, myThid )
                0101        ENDIF
cb5aad258c Jean*0102       ENDIF
                0103       IF ( balanceQnet  .AND. (.NOT.useSeaice .OR. useThSIce) ) THEN
7e00d7e8f9 Jean*0104         tmpVar(1) = oneRS
                0105         CALL REMOVE_MEAN_RS( 1, Qnet,  maskInC, maskInC, rA,
                0106      &                       tmpVar, 'Qnet ', myTime, myIter, myThid )
cb5aad258c Jean*0107       ENDIF
                0108 #endif /* ALLOW_BALANCE_FLUXES */
                0109
                0110 C-    Apply exchanges (if needed)
                0111
                0112 #ifdef CHECK_OVERLAP_FORCING
                0113 C     Put large value in overlap of forcing array to check if exch is needed
                0114 c     IF ( .NOT. useKPP ) THEN
                0115        fixVal = 1.
                0116        CALL RESET_HALO_RS ( EmPmR, fixVal, 1, myThid )
                0117        fixVal = 400.
                0118        CALL RESET_HALO_RS ( Qnet, fixVal, 1, myThid )
                0119        fixVal = -200.
                0120        CALL RESET_HALO_RS ( Qsw, fixVal, 1, myThid )
                0121        fixVal = 40.
                0122        CALL RESET_HALO_RS ( saltFlux, fixVal, 1, myThid )
                0123 c     ENDIF
                0124 #endif /* CHECK_OVERLAP_FORCING */
042348d828 Jean*0125
cb5aad258c Jean*0126 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
042348d828 Jean*0127
cb5aad258c Jean*0128 #ifdef EXACT_CONSERV
                0129 C NB: synchronous time step: PmEpR lag 1 time step behind EmPmR
                0130 C     to stay consitent with volume change (=d/dt etaH).
62d9359ab4 Patr*0131 # ifdef ALLOW_AUTODIFF_TAMC
cb5aad258c Jean*0132 CADJ STORE PmEpR = comlev1, key = ikey_dynamics,  kind = isbyte
4f48a680f8 Gael*0133 CADJ STORE EmPmR = comlev1, key = ikey_dynamics,  kind = isbyte
62d9359ab4 Patr*0134 # endif
cb5aad258c Jean*0135       DO bj=myByLo(myThid),myByHi(myThid)
                0136        DO bi=myBxLo(myThid),myBxHi(myThid)
                0137         IF ( staggerTimeStep ) THEN
                0138          DO j=1-OLy,sNy+OLy
                0139           DO i=1-OLx,sNx+OLx
                0140            PmEpR(i,j,bi,bj) = -EmPmR(i,j,bi,bj)
042348d828 Jean*0141           ENDDO
                0142          ENDDO
                0143         ENDIF
9e3a303fa3 Gael*0144        ENDDO
cb5aad258c Jean*0145       ENDDO
                0146 #endif /* EXACT_CONSERV */
042348d828 Jean*0147
cb5aad258c Jean*0148 C--   set surfaceForcingT,S to zero.
                0149       DO bj=myByLo(myThid),myByHi(myThid)
                0150        DO bi=myBxLo(myThid),myBxHi(myThid)
                0151         DO j=1-OLy,sNy+OLy
                0152          DO i=1-OLx,sNx+OLx
                0153            surfaceForcingT(i,j,bi,bj) = 0. _d 0
                0154            surfaceForcingS(i,j,bi,bj) = 0. _d 0
                0155          ENDDO
042348d828 Jean*0156         ENDDO
                0157        ENDDO
cb5aad258c Jean*0158       ENDDO
                0159
                0160 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
042348d828 Jean*0161
cb5aad258c Jean*0162 C--   Start with surface restoring term :
                0163       IF ( doThetaClimRelax .OR. doSaltClimRelax ) THEN
                0164        CALL FORCING_SURF_RELAX(
                0165      I              iMin, iMax, jMin, jMax,
                0166      I              myTime, myIter, myThid )
042348d828 Jean*0167       ENDIF
                0168
ea7a86e491 Jean*0169 #ifdef ALLOW_PTRACERS
cb5aad258c Jean*0170 C--   passive tracer surface forcing:
                0171 #ifdef ALLOW_AUTODIFF_TAMC
                0172 CADJ STORE surfaceForcingS = comlev1, key = ikey_dynamics,
                0173 CADJ &    kind = isbyte
abd2e981eb Matt*0174 #ifdef ALLOW_BLING
                0175 CADJ STORE EmPmR = comlev1, key = ikey_dynamics,  kind = isbyte
                0176 #endif
cb5aad258c Jean*0177 #endif
ea7a86e491 Jean*0178       IF ( usePTRACERS ) THEN
cb5aad258c Jean*0179        DO bj=myByLo(myThid),myByHi(myThid)
                0180         DO bi=myBxLo(myThid),myBxHi(myThid)
                0181          CALL PTRACERS_FORCING_SURF(
                0182      I        surfaceForcingS(1-OLx,1-OLy,bi,bj),
                0183      I        bi, bj, iMin, iMax, jMin, jMax,
e3e2f69df3 Jean*0184      I        myTime, myIter, myThid )
ea7a86e491 Jean*0185         ENDDO
                0186        ENDDO
                0187       ENDIF
                0188 #endif /* ALLOW_PTRACERS */
                0189
cb5aad258c Jean*0190 C- Notes: setting of PmEpR and pTracers surface forcing could have been
                0191 C         moved below, inside a unique bi,bj block. However this results
                0192 C         in tricky dependencies for TAF (and recomputations).
042348d828 Jean*0193 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
e305438401 Mart*0194
cb5aad258c Jean*0195       DO bj=myByLo(myThid),myByHi(myThid)
                0196        DO bi=myBxLo(myThid),myBxHi(myThid)
1d26153115 Jean*0197
4f48a680f8 Gael*0198 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0199         tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
4f48a680f8 Gael*0200 #endif /* ALLOW_AUTODIFF_TAMC */
                0201
                0202 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0203 CADJ STORE EmPmR(:,:,bi,bj) = comlev1_bibj, key=tkey,  kind = isbyte
4f48a680f8 Gael*0204 #ifdef EXACT_CONSERV
edb6656069 Mart*0205 CADJ STORE PmEpR(:,:,bi,bj) = comlev1_bibj, key=tkey,  kind = isbyte
4f48a680f8 Gael*0206 #endif
                0207 #endif /* ALLOW_AUTODIFF_TAMC */
                0208
cb5aad258c Jean*0209 C--   Surface Fluxes :
                0210         DO j = jMin, jMax
51701379d8 Ed H*0211          DO i = iMin, iMax
77af23a186 Patr*0212
042348d828 Jean*0213 C     Zonal wind stress fu:
cb5aad258c Jean*0214           surfaceForcingU(i,j,bi,bj) = fu(i,j,bi,bj)*mass2rUnit
042348d828 Jean*0215 C     Meridional wind stress fv:
cb5aad258c Jean*0216           surfaceForcingV(i,j,bi,bj) = fv(i,j,bi,bj)*mass2rUnit
042348d828 Jean*0217 C     Net heat flux Qnet:
                0218           surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
2d2cc93d4f Jean*0219      &       - ( Qnet(i,j,bi,bj)
                0220 #ifdef SHORTWAVE_HEATING
042348d828 Jean*0221      &          -Qsw(i,j,bi,bj)
2d2cc93d4f Jean*0222 #endif
faf82d94de Patr*0223      &         ) *recip_Cp*mass2rUnit
21b8df1d84 Jean*0224 C     Net Salt Flux :
042348d828 Jean*0225           surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
0b1017b546 Jean*0226      &      -saltFlux(i,j,bi,bj)*mass2rUnit
e4775240e5 Dimi*0227
                0228          ENDDO
cb5aad258c Jean*0229         ENDDO
e4775240e5 Dimi*0230
8c3259a14c Dimi*0231 #ifdef ALLOW_SALT_PLUME
                0232 C saltPlume is the amount of salt rejected by ice while freezing;
                0233 C it is here subtracted from surfaceForcingS and will be redistributed
                0234 C to multiple vertical levels later on as per Duffy et al. (GRL 1999)
1f89baba18 Patr*0235 C-- for the case of SALT_PLUME_VOLUME, need to call this S/R right
                0236 C-- before kpp in do_oceanic_phys.F due to recent moved of
                0237 C-- external_forcing_surf.F outside bi,bj loop.
                0238 #ifndef SALT_PLUME_VOLUME
cb5aad258c Jean*0239         IF ( useSALT_PLUME ) THEN
e4775240e5 Dimi*0240          CALL SALT_PLUME_FORCING_SURF(
                0241      I        bi, bj, iMin, iMax, jMin, jMax,
e3e2f69df3 Jean*0242      I        myTime, myIter, myThid )
cb5aad258c Jean*0243         ENDIF
1f89baba18 Patr*0244 #endif /* SALT_PLUME_VOLUME */
8c3259a14c Dimi*0245 #endif /* ALLOW_SALT_PLUME */
51701379d8 Ed H*0246
fab93eca08 Jean*0247 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
21b8df1d84 Jean*0248 C--   Fresh-water flux
77af23a186 Patr*0249
bf53796fe2 Jean*0250 C-    Apply mask on Fresh-Water flux (if useRealFreshWaterFlux)
                0251 C     <== removed: maskInC is applied directly in S/R SOLVE_FOR_PRESSURE
2aafdf1c7e Jean*0252
fab93eca08 Jean*0253 #ifdef EXACT_CONSERV
9669509dca Jean*0254       IF ( (nonlinFreeSurf.GT.0 .OR. usingPCoords)
51701379d8 Ed H*0255      &     .AND. useRealFreshWaterFlux ) THEN
fa41ecc867 Jean*0256
7f043fb97f Jean*0257 C--   NonLin_FrSurf and RealFreshWaterFlux : PmEpR effectively changes
                0258 C     the water column height ; temp., salt, (tracer) flux associated
                0259 C     with this input/output of water is added here to the surface tendency.
fa41ecc867 Jean*0260
51701379d8 Ed H*0261        IF (temp_EvPrRn.NE.UNSET_RL) THEN
                0262         DO j = jMin, jMax
                0263          DO i = iMin, iMax
042348d828 Jean*0264           surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
51701379d8 Ed H*0265      &      + PmEpR(i,j,bi,bj)
7f043fb97f Jean*0266      &          *( temp_EvPrRn - theta(i,j,ks,bi,bj) )
62fd6ae4e5 Jean*0267      &          *mass2rUnit
51701379d8 Ed H*0268          ENDDO
                0269         ENDDO
                0270        ENDIF
                0271
                0272        IF (salt_EvPrRn.NE.UNSET_RL) THEN
                0273         DO j = jMin, jMax
                0274          DO i = iMin, iMax
042348d828 Jean*0275           surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
7f043fb97f Jean*0276      &      + PmEpR(i,j,bi,bj)
                0277      &          *( salt_EvPrRn - salt(i,j,ks,bi,bj) )
62fd6ae4e5 Jean*0278      &          *mass2rUnit
51701379d8 Ed H*0279          ENDDO
                0280         ENDDO
                0281        ENDIF
fa41ecc867 Jean*0282
463053c692 Jean*0283 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
51701379d8 Ed H*0284       ELSE
fab93eca08 Jean*0285 #else /* EXACT_CONSERV */
51701379d8 Ed H*0286       IF (.TRUE.) THEN
fab93eca08 Jean*0287 #endif /* EXACT_CONSERV */
463053c692 Jean*0288
7f043fb97f Jean*0289 C--   EmPmR does not really affect the water column height (for tracer budget)
                0290 C     and is converted to a salt tendency.
463053c692 Jean*0291
51701379d8 Ed H*0292        IF (convertFW2Salt .EQ. -1.) THEN
7f043fb97f Jean*0293 C-    use local surface tracer field to calculate forcing term:
                0294
                0295         IF (temp_EvPrRn.NE.UNSET_RL) THEN
                0296 C     account for Rain/Evap heat content (temp_EvPrRn) using local SST
                0297          DO j = jMin, jMax
                0298           DO i = iMin, iMax
                0299            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0300      &       + EmPmR(i,j,bi,bj)
                0301      &           *( theta(i,j,ks,bi,bj) - temp_EvPrRn )
62fd6ae4e5 Jean*0302      &           *mass2rUnit
7f043fb97f Jean*0303           ENDDO
51701379d8 Ed H*0304          ENDDO
7f043fb97f Jean*0305         ENDIF
                0306         IF (salt_EvPrRn.NE.UNSET_RL) THEN
                0307 C     converts EmPmR to salinity tendency using surface local salinity
                0308          DO j = jMin, jMax
                0309           DO i = iMin, iMax
                0310            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0311      &       + EmPmR(i,j,bi,bj)
                0312      &           *( salt(i,j,ks,bi,bj) - salt_EvPrRn )
62fd6ae4e5 Jean*0313      &           *mass2rUnit
7f043fb97f Jean*0314           ENDDO
                0315          ENDDO
                0316         ENDIF
                0317
21b8df1d84 Jean*0318        ELSE
7f043fb97f Jean*0319 C-    use uniform tracer value to calculate forcing term:
                0320
                0321         IF (temp_EvPrRn.NE.UNSET_RL) THEN
                0322 C     account for Rain/Evap heat content (temp_EvPrRn) assuming uniform SST (=tRef)
                0323          DO j = jMin, jMax
                0324           DO i = iMin, iMax
                0325            surfaceForcingT(i,j,bi,bj) = surfaceForcingT(i,j,bi,bj)
                0326      &       + EmPmR(i,j,bi,bj)
                0327      &           *( tRef(ks) - temp_EvPrRn )
62fd6ae4e5 Jean*0328      &           *mass2rUnit
7f043fb97f Jean*0329           ENDDO
51701379d8 Ed H*0330          ENDDO
7f043fb97f Jean*0331         ENDIF
                0332         IF (salt_EvPrRn.NE.UNSET_RL) THEN
                0333 C     converts EmPmR to virtual salt flux using uniform salinity (default=35)
                0334          DO j = jMin, jMax
                0335           DO i = iMin, iMax
                0336            surfaceForcingS(i,j,bi,bj) = surfaceForcingS(i,j,bi,bj)
                0337      &       + EmPmR(i,j,bi,bj)
                0338      &           *( convertFW2Salt - salt_EvPrRn )
62fd6ae4e5 Jean*0339      &           *mass2rUnit
7f043fb97f Jean*0340           ENDDO
                0341          ENDDO
                0342         ENDIF
                0343
                0344 C-    end local-surface-tracer / uniform-value distinction
51701379d8 Ed H*0345        ENDIF
463053c692 Jean*0346
51701379d8 Ed H*0347       ENDIF
d8206d87ee Patr*0348
463053c692 Jean*0349 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
25eac65bee Jean*0350
463053c692 Jean*0351 #ifdef ATMOSPHERIC_LOADING
a5e98ae15f Jean*0352 C-- Atmospheric surface Pressure loading : added to phi0surf when using Z-coord;
                0353 C   Not yet implemented for Ocean in P: would need to be applied to the other end
                0354 C   of the column, as a vertical velocity (omega); (meaningless for Atmos in P).
                0355 C- Note:
                0356 C   Using P-coord., a hack (now directly applied from S/R INI_FORCING)
                0357 C   is sometime used to read phi0surf from a file (pLoadFile) instead
                0358 C   of computing it from bathymetry & density ref. profile.
463053c692 Jean*0359
cb5aad258c Jean*0360         IF ( usingZCoords ) THEN
a5e98ae15f Jean*0361 C   The true atmospheric P-loading is not yet implemented for P-coord
                0362 C   (requires time varying dP(Nr) like dP(k-bottom) with NonLin FS).
cb5aad258c Jean*0363          IF ( useRealFreshWaterFlux ) THEN
                0364           DO j = jMin, jMax
                0365            DO i = iMin, iMax
                0366             phi0surf(i,j,bi,bj) = ( pLoad(i,j,bi,bj)
0320e25227 Mart*0367      &                          +sIceLoad(i,j,bi,bj)*gravity*sIceLoadFac
cb5aad258c Jean*0368      &                            )*recip_rhoConst
                0369            ENDDO
                0370           ENDDO
                0371          ELSE
                0372           DO j = jMin, jMax
                0373            DO i = iMin, iMax
                0374             phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)*recip_rhoConst
                0375            ENDDO
                0376           ENDDO
                0377          ENDIF
                0378 c       ELSEIF ( usingPCoords ) THEN
a5e98ae15f Jean*0379 C-- This is a hack used to read phi0surf from a file (pLoadFile)
463053c692 Jean*0380 C   instead of computing it from bathymetry & density ref. profile.
a5e98ae15f Jean*0381 C   ==> now done only once, in S/R INI_FORCING
cb5aad258c Jean*0382 c         DO j = jMin, jMax
                0383 c          DO i = iMin, iMax
                0384 c           phi0surf(i,j,bi,bj) = pLoad(i,j,bi,bj)
                0385 c          ENDDO
                0386 c         ENDDO
                0387         ENDIF
463053c692 Jean*0388 #endif /* ATMOSPHERIC_LOADING */
fa41ecc867 Jean*0389
03759a535c Mart*0390 #ifdef ALLOW_SHELFICE
e3e2f69df3 Jean*0391         IF ( useSHELFICE) THEN
                0392           CALL SHELFICE_FORCING_SURF(
                0393      I                  bi, bj, iMin, iMax, jMin, jMax,
                0394      I                  myTime, myIter, myThid )
cb5aad258c Jean*0395         ENDIF
03759a535c Mart*0396 #endif /* ALLOW_SHELFICE */
                0397
cb5aad258c Jean*0398 C--   end bi,bj loops.
                0399        ENDDO
                0400       ENDDO
                0401
77af23a186 Patr*0402       RETURN
                0403       END