Back to home page

MITgcm
 
 

    

File indexing completed on 2025-09-19 05:09:07 UTC

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