Back to home page

MITgcm

 
 

    


File indexing completed on 2022-04-25 05:09:44 UTC

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