** Warning **

Issuing rollback() due to DESTROY without explicit disconnect() of DBD::mysql::db handle dbname=MITgcm at /usr/local/share/lxr/lib/LXR/Common.pm line 1224.

Last-Modified: Thu, 15 May 2024 05:11:27 GMT Content-Type: text/html; charset=utf-8 MITgcm/MITgcm/model/src/apply_forcing.F
Back to home page

MITgcm

 
 

    


File indexing completed on 2021-09-28 05:14:31 UTC

view on githubraw file Latest commit 44d39862 on 2021-09-27 19:44:02 UTC
e874fa47e5 Jean*0001 #include "PACKAGES_CONFIG.h"
                0002 #include "CPP_OPTIONS.h"
                0003 
                0004 C--  File apply_forcing.F:
                0005 C--   Contents
                0006 C--   o APPLY_FORCING_U
                0007 C--   o APPLY_FORCING_V
                0008 C--   o APPLY_FORCING_T
                0009 C--   o APPLY_FORCING_S
                0010 
                0011 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0012 CBOP
                0013 C     !ROUTINE: APPLY_FORCING_U
                0014 C     !INTERFACE:
                0015       SUBROUTINE APPLY_FORCING_U(
                0016      U                     gU_arr,
                0017      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0018      I                     myTime, myIter, myThid )
                0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
                0021 C     | S/R APPLY_FORCING_U
                0022 C     | o Contains problem specific forcing for zonal velocity.
                0023 C     *==========================================================*
                0024 C     | Adds terms to gU for forcing by external sources
                0025 C     | e.g. wind stress, bottom friction etc ...
                0026 C     *==========================================================*
                0027 C     \ev
                0028 
                0029 C     !USES:
                0030       IMPLICIT NONE
                0031 C     == Global data ==
                0032 #include "SIZE.h"
                0033 #include "EEPARAMS.h"
                0034 #include "PARAMS.h"
                0035 #include "GRID.h"
                0036 #include "DYNVARS.h"
                0037 #include "FFIELDS.h"
                0038 
                0039 C     !INPUT/OUTPUT PARAMETERS:
                0040 C     gU_arr    :: the tendency array
                0041 C     iMin,iMax :: Working range of x-index for applying forcing.
                0042 C     jMin,jMax :: Working range of y-index for applying forcing.
                0043 C     k         :: Current vertical level index
                0044 C     bi,bj     :: Current tile indices
                0045 C     myTime    :: Current time in simulation
                0046 C     myIter    :: Current iteration number
                0047 C     myThid    :: my Thread Id number
                0048       _RL     gU_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049       INTEGER iMin, iMax, jMin, jMax
                0050       INTEGER k, bi, bj
                0051       _RL     myTime
                0052       INTEGER myIter
                0053       INTEGER myThid
                0054 
                0055 C     !LOCAL VARIABLES:
                0056 C     i,j       :: Loop counters
                0057 C     kSurface  :: index of surface level
                0058       INTEGER i, j
                0059 #ifdef USE_OLD_EXTERNAL_FORCING
                0060       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
202e12438b Jean*0061       _RL     tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e874fa47e5 Jean*0062 #else
                0063       INTEGER kSurface
5a705ed756 Jean*0064 #endif /* USE_OLD_EXTERNAL_FORCING */
e874fa47e5 Jean*0065 CEOP
                0066 
                0067 #ifdef USE_OLD_EXTERNAL_FORCING
                0068 
                0069       DO j=1-OLy,sNy+OLy
                0070         DO i=1-OLx,sNx+OLx
                0071           locVar(i,j) = gU(i,j,k,bi,bj)
                0072         ENDDO
                0073       ENDDO
                0074       CALL EXTERNAL_FORCING_U(
                0075      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0076      I              myTime, myThid )
202e12438b Jean*0077 C-    Use 2-d local array tmpVar and split loop in 2 parts to avoid compiler
                0078 C     to mess-up this part by re-arranging the order of instructions (wrong
                0079 C     when gU and gU_arr are the same array, i.e., called with argument gU).
e874fa47e5 Jean*0080       DO j=1-OLy,sNy+OLy
                0081         DO i=1-OLx,sNx+OLx
202e12438b Jean*0082           tmpVar(i,j) = gU(i,j,k,bi,bj) - locVar(i,j)
e874fa47e5 Jean*0083           gU(i,j,k,bi,bj) = locVar(i,j)
202e12438b Jean*0084         ENDDO
                0085       ENDDO
aa25968b23 Jean*0086 C-    not needed since APPLY_FORCING_U is no longer called with argument gU
                0087 c     CALL FOOL_THE_COMPILER_RL( tmpVar(1,1) )
202e12438b Jean*0088       DO j=1-OLy,sNy+OLy
                0089         DO i=1-OLx,sNx+OLx
                0090           gU_arr(i,j) = gU_arr(i,j) + tmpVar(i,j)
e874fa47e5 Jean*0091         ENDDO
                0092       ENDDO
                0093 
                0094 #else  /* USE_OLD_EXTERNAL_FORCING */
                0095 
                0096       IF ( fluidIsAir ) THEN
                0097        kSurface = 0
                0098       ELSEIF ( usingPCoords ) THEN
                0099        kSurface = Nr
                0100       ELSE
                0101        kSurface = 1
                0102       ENDIF
                0103 
                0104 C--   Forcing term
                0105 #ifdef ALLOW_AIM
                0106       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
                0107      U                       gU_arr,
                0108      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0109      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0110 #endif /* ALLOW_AIM */
                0111 
                0112 #ifdef ALLOW_ATM_PHYS
                0113       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
                0114      U                       gU_arr,
                0115      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0116      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0117 #endif /* ALLOW_ATM_PHYS */
                0118 
                0119 #ifdef ALLOW_FIZHI
                0120       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
                0121      U                       gU_arr,
                0122      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0123      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0124 #endif /* ALLOW_FIZHI */
                0125 
614bfb3d2a Jean*0126 C     Add Tidal momentum forcing from 2-d geopotential anomaly
                0127       IF ( momTidalForcing ) THEN
                0128        DO j=0,sNy+1
                0129         DO i=1,sNx+1
                0130           gU_arr(i,j) = gU_arr(i,j)
                0131      &      - recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0132      &      * ( phiTide2d(i,j,bi,bj) - phiTide2d(i-1,j,bi,bj) )
                0133      &      *_maskW(i,j,k,bi,bj)
                0134         ENDDO
                0135        ENDDO
                0136       ENDIF
                0137 
e874fa47e5 Jean*0138 C     Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
                0139       IF ( k .EQ. kSurface ) THEN
                0140 c      DO j=1,sNy
                0141 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
                0142        DO j=0,sNy+1
                0143         DO i=1,sNx+1
                0144           gU_arr(i,j) = gU_arr(i,j)
                0145      &      +foFacMom*surfaceForcingU(i,j,bi,bj)
                0146      &      *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
                0147         ENDDO
                0148        ENDDO
                0149       ELSEIF ( kSurface.EQ.-1 ) THEN
                0150        DO j=0,sNy+1
                0151         DO i=1,sNx+1
                0152          IF ( kSurfW(i,j,bi,bj).EQ.k ) THEN
                0153           gU_arr(i,j) = gU_arr(i,j)
                0154      &      +foFacMom*surfaceForcingU(i,j,bi,bj)
                0155      &      *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
                0156          ENDIF
                0157         ENDDO
                0158        ENDDO
                0159       ENDIF
                0160 
                0161 #ifdef ALLOW_EDDYPSI
                0162          CALL TAUEDDY_TENDENCY_APPLY_U(
                0163      U                 gU_arr,
                0164      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0165      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0166 #endif
                0167 
                0168 #ifdef ALLOW_RBCS
                0169       IF (useRBCS) THEN
                0170         CALL RBCS_ADD_TENDENCY(
                0171      U                 gU_arr,
                0172      I                 k, bi, bj, -1,
2c160c3ab4 Jean*0173      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0174 
                0175       ENDIF
                0176 #endif /* ALLOW_RBCS */
                0177 
                0178 #ifdef ALLOW_OBCS
                0179       IF (useOBCS) THEN
                0180         CALL OBCS_SPONGE_U(
                0181      U                   gU_arr,
                0182      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0183      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0184       ENDIF
                0185 #endif /* ALLOW_OBCS */
                0186 
                0187 #ifdef ALLOW_MYPACKAGE
                0188       IF ( useMYPACKAGE ) THEN
                0189         CALL MYPACKAGE_TENDENCY_APPLY_U(
                0190      U                 gU_arr,
                0191      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0192      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0193       ENDIF
                0194 #endif /* ALLOW_MYPACKAGE */
                0195 
                0196 #endif /* USE_OLD_EXTERNAL_FORCING */
                0197 
                0198       RETURN
                0199       END
                0200 
                0201 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0202 CBOP
                0203 C     !ROUTINE: APPLY_FORCING_V
                0204 C     !INTERFACE:
                0205       SUBROUTINE APPLY_FORCING_V(
                0206      U                     gV_arr,
                0207      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0208      I                     myTime, myIter, myThid )
                0209 C     !DESCRIPTION: \bv
                0210 C     *==========================================================*
                0211 C     | S/R APPLY_FORCING_V
                0212 C     | o Contains problem specific forcing for merid velocity.
                0213 C     *==========================================================*
                0214 C     | Adds terms to gV for forcing by external sources
                0215 C     | e.g. wind stress, bottom friction etc ...
                0216 C     *==========================================================*
                0217 C     \ev
                0218 
                0219 C     !USES:
                0220       IMPLICIT NONE
                0221 C     == Global data ==
                0222 #include "SIZE.h"
                0223 #include "EEPARAMS.h"
                0224 #include "PARAMS.h"
                0225 #include "GRID.h"
                0226 #include "DYNVARS.h"
                0227 #include "FFIELDS.h"
                0228 
                0229 C     !INPUT/OUTPUT PARAMETERS:
                0230 C     gV_arr    :: the tendency array
                0231 C     iMin,iMax :: Working range of x-index for applying forcing.
                0232 C     jMin,jMax :: Working range of y-index for applying forcing.
                0233 C     k         :: Current vertical level index
                0234 C     bi,bj     :: Current tile indices
                0235 C     myTime    :: Current time in simulation
                0236 C     myIter    :: Current iteration number
                0237 C     myThid    :: my Thread Id number
                0238       _RL     gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0239       INTEGER iMin, iMax, jMin, jMax
                0240       INTEGER k, bi, bj
                0241       _RL     myTime
                0242       INTEGER myIter
                0243       INTEGER myThid
                0244 
                0245 C     !LOCAL VARIABLES:
                0246 C     i,j       :: Loop counters
                0247 C     kSurface  :: index of surface level
                0248       INTEGER i, j
                0249 #ifdef USE_OLD_EXTERNAL_FORCING
                0250       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
202e12438b Jean*0251       _RL     tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e874fa47e5 Jean*0252 #else
                0253       INTEGER kSurface
5a705ed756 Jean*0254 #endif /* USE_OLD_EXTERNAL_FORCING */
e874fa47e5 Jean*0255 CEOP
                0256 
                0257 #ifdef USE_OLD_EXTERNAL_FORCING
                0258 
                0259       DO j=1-OLy,sNy+OLy
                0260         DO i=1-OLx,sNx+OLx
                0261           locVar(i,j) = gV(i,j,k,bi,bj)
                0262         ENDDO
                0263       ENDDO
                0264       CALL EXTERNAL_FORCING_V(
                0265      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0266      I              myTime, myThid )
202e12438b Jean*0267 C-    Use 2-d local array tmpVar and split loop in 2 parts to avoid compiler
                0268 C     to mess-up this part by re-arranging the order of instructions (wrong
                0269 C     when gV and gV_arr are the same array, i.e., called with argument gV).
e874fa47e5 Jean*0270       DO j=1-OLy,sNy+OLy
                0271         DO i=1-OLx,sNx+OLx
202e12438b Jean*0272           tmpVar(i,j) = gV(i,j,k,bi,bj) - locVar(i,j)
e874fa47e5 Jean*0273           gV(i,j,k,bi,bj) = locVar(i,j)
202e12438b Jean*0274         ENDDO
                0275       ENDDO
aa25968b23 Jean*0276 C-    not needed since APPLY_FORCING_V is no longer called with argument gV
                0277 c     CALL FOOL_THE_COMPILER_RL( tmpVar(1,1) )
202e12438b Jean*0278       DO j=1-OLy,sNy+OLy
                0279         DO i=1-OLx,sNx+OLx
                0280           gV_arr(i,j) = gV_arr(i,j) + tmpVar(i,j)
e874fa47e5 Jean*0281         ENDDO
                0282       ENDDO
                0283 
                0284 #else  /* USE_OLD_EXTERNAL_FORCING */
                0285 
                0286       IF ( fluidIsAir ) THEN
                0287        kSurface = 0
                0288       ELSEIF ( usingPCoords ) THEN
                0289        kSurface = Nr
                0290       ELSE
                0291        kSurface = 1
                0292       ENDIF
                0293 
                0294 C--   Forcing term
                0295 #ifdef ALLOW_AIM
                0296       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
                0297      U                       gV_arr,
                0298      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0299      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0300 #endif /* ALLOW_AIM */
                0301 
                0302 #ifdef ALLOW_ATM_PHYS
                0303       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
                0304      U                       gV_arr,
                0305      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0306      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0307 #endif /* ALLOW_ATM_PHYS */
                0308 
                0309 #ifdef ALLOW_FIZHI
                0310       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
                0311      U                       gV_arr,
                0312      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0313      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0314 #endif /* ALLOW_FIZHI */
                0315 
614bfb3d2a Jean*0316 C     Add Tidal momentum forcing from 2-d geopotential anomaly
                0317       IF ( momTidalForcing ) THEN
                0318        DO j=1,sNy+1
                0319         DO i=0,sNx+1
                0320           gV_arr(i,j) = gV_arr(i,j)
                0321      &      - recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0322      &      *( phiTide2d(i,j,bi,bj) - phiTide2d(i,j-1,bi,bj) )
                0323      &      *_maskS(i,j,k,bi,bj)
                0324         ENDDO
                0325        ENDDO
                0326       ENDIF
                0327 
e874fa47e5 Jean*0328 C     Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
                0329       IF ( k .EQ. kSurface ) THEN
                0330        DO j=1,sNy+1
                0331 c       DO i=1,sNx
                0332 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
                0333         DO i=0,sNx+1
                0334           gV_arr(i,j) = gV_arr(i,j)
                0335      &      +foFacMom*surfaceForcingV(i,j,bi,bj)
                0336      &      *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
                0337         ENDDO
                0338        ENDDO
                0339       ELSEIF ( kSurface.EQ.-1 ) THEN
                0340        DO j=1,sNy+1
                0341         DO i=0,sNx+1
                0342          IF ( kSurfS(i,j,bi,bj).EQ.k ) THEN
                0343           gV_arr(i,j) = gV_arr(i,j)
                0344      &      +foFacMom*surfaceForcingV(i,j,bi,bj)
                0345      &      *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
                0346          ENDIF
                0347         ENDDO
                0348        ENDDO
                0349       ENDIF
                0350 
                0351 #ifdef ALLOW_EDDYPSI
                0352          CALL TAUEDDY_TENDENCY_APPLY_V(
                0353      U                 gV_arr,
                0354      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0355      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0356 #endif
                0357 
                0358 #ifdef ALLOW_RBCS
                0359       IF (useRBCS) THEN
                0360         CALL RBCS_ADD_TENDENCY(
                0361      U                 gV_arr,
                0362      I                 k, bi, bj, -2,
2c160c3ab4 Jean*0363      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0364       ENDIF
                0365 #endif /* ALLOW_RBCS */
                0366 
                0367 #ifdef ALLOW_OBCS
                0368       IF (useOBCS) THEN
                0369         CALL OBCS_SPONGE_V(
                0370      U                   gV_arr,
                0371      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0372      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0373       ENDIF
                0374 #endif /* ALLOW_OBCS */
                0375 
                0376 #ifdef ALLOW_MYPACKAGE
                0377       IF ( useMYPACKAGE ) THEN
                0378         CALL MYPACKAGE_TENDENCY_APPLY_V(
                0379      U                 gV_arr,
                0380      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0381      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0382       ENDIF
                0383 #endif /* ALLOW_MYPACKAGE */
                0384 
                0385 #endif /* USE_OLD_EXTERNAL_FORCING */
                0386 
                0387       RETURN
                0388       END
                0389 
                0390 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0391 CBOP
                0392 C     !ROUTINE: APPLY_FORCING_T
                0393 C     !INTERFACE:
                0394       SUBROUTINE APPLY_FORCING_T(
                0395      U                     gT_arr,
                0396      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0397      I                     myTime, myIter, myThid )
                0398 C     !DESCRIPTION: \bv
                0399 C     *==========================================================*
                0400 C     | S/R APPLY_FORCING_T
                0401 C     | o Contains problem specific forcing for temperature.
                0402 C     *==========================================================*
                0403 C     | Adds terms to gT for forcing by external sources
                0404 C     | e.g. heat flux, climatalogical relaxation, etc ...
                0405 C     *==========================================================*
                0406 C     \ev
                0407 
                0408 C     !USES:
                0409       IMPLICIT NONE
                0410 C     == Global data ==
                0411 #include "SIZE.h"
                0412 #include "EEPARAMS.h"
                0413 #include "PARAMS.h"
                0414 #include "GRID.h"
                0415 #include "DYNVARS.h"
                0416 #include "FFIELDS.h"
                0417 #include "SURFACE.h"
                0418 
                0419 C     !INPUT/OUTPUT PARAMETERS:
                0420 C     gT_arr    :: the tendency array
                0421 C     iMin,iMax :: Working range of x-index for applying forcing.
                0422 C     jMin,jMax :: Working range of y-index for applying forcing.
                0423 C     k         :: Current vertical level index
                0424 C     bi,bj     :: Current tile indices
                0425 C     myTime    :: Current time in simulation
                0426 C     myIter    :: Current iteration number
                0427 C     myThid    :: my Thread Id number
                0428       _RL     gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0429       INTEGER iMin, iMax, jMin, jMax
                0430       INTEGER k, bi, bj
                0431       _RL     myTime
                0432       INTEGER myIter
                0433       INTEGER myThid
                0434 
                0435 C     !LOCAL VARIABLES:
                0436 C     i,j       :: Loop counters
                0437 C     kSurface  :: index of surface level
                0438       INTEGER i, j
5a705ed756 Jean*0439 #ifndef USE_OLD_EXTERNAL_FORCING
e874fa47e5 Jean*0440       INTEGER kSurface
                0441       INTEGER km, kc, kp
5a705ed756 Jean*0442       _RL     tmpVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
e874fa47e5 Jean*0443       _RL tmpFac, delPI
                0444       _RL recip_Cp
                0445 #ifdef SHORTWAVE_HEATING
                0446       _RL swfracb(2)
                0447 #endif
5a705ed756 Jean*0448 #endif /* USE_OLD_EXTERNAL_FORCING */
                0449 CEOP
e874fa47e5 Jean*0450 
                0451 #ifdef USE_OLD_EXTERNAL_FORCING
                0452 
                0453       DO j=1-OLy,sNy+OLy
                0454         DO i=1-OLx,sNx+OLx
5a705ed756 Jean*0455           gT(i,j,k,bi,bj) = 0. _d 0
e874fa47e5 Jean*0456         ENDDO
                0457       ENDDO
                0458       CALL EXTERNAL_FORCING_T(
                0459      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0460      I              myTime, myThid )
202e12438b Jean*0461       DO j=1-OLy,sNy+OLy
                0462         DO i=1-OLx,sNx+OLx
5a705ed756 Jean*0463           gT_arr(i,j) = gT_arr(i,j) + gT(i,j,k,bi,bj)
e874fa47e5 Jean*0464         ENDDO
                0465       ENDDO
                0466 
                0467 #else  /* USE_OLD_EXTERNAL_FORCING */
                0468 
                0469       IF ( fluidIsAir ) THEN
                0470        kSurface = 0
                0471       ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
                0472        kSurface = -1
                0473       ELSEIF ( usingPCoords ) THEN
                0474        kSurface = Nr
                0475       ELSE
                0476        kSurface = 1
                0477       ENDIF
                0478       recip_Cp = 1. _d 0 / HeatCapacity_Cp
                0479 
5a705ed756 Jean*0480 C--   Note on loop range: For model dynamics, only needs to get correct
                0481 C     forcing (update gT_arr) in tile interior (1:sNx,1:sNy);
                0482 C     However, for some diagnostics, we may want to get valid forcing
                0483 C     extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
                0484 
e874fa47e5 Jean*0485 C--   Forcing term
                0486 #ifdef ALLOW_AIM
                0487       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
                0488      U                       gT_arr,
                0489      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0490      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0491 #endif /* ALLOW_AIM */
                0492 
                0493 #ifdef ALLOW_ATM_PHYS
                0494       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
                0495      U                       gT_arr,
                0496      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0497      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0498 #endif /* ALLOW_ATM_PHYS */
                0499 
                0500 #ifdef ALLOW_FIZHI
                0501       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
                0502      U                       gT_arr,
                0503      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0504      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0505 #endif /* ALLOW_FIZHI */
                0506 
                0507 #ifdef ALLOW_ADDFLUID
                0508       IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
                0509        IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
                0510      &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
5a705ed756 Jean*0511          DO j=0,sNy+1
                0512           DO i=0,sNx+1
e874fa47e5 Jean*0513             gT_arr(i,j) = gT_arr(i,j)
                0514      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0515      &          *( temp_addMass - theta(i,j,k,bi,bj) )
                0516      &          *recip_rA(i,j,bi,bj)
                0517      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0518 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
                0519           ENDDO
                0520          ENDDO
                0521        ELSE
5a705ed756 Jean*0522          DO j=0,sNy+1
                0523           DO i=0,sNx+1
e874fa47e5 Jean*0524             gT_arr(i,j) = gT_arr(i,j)
                0525      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0526      &          *( temp_addMass - tRef(k) )
                0527      &          *recip_rA(i,j,bi,bj)
                0528      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0529 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
                0530           ENDDO
                0531          ENDDO
                0532        ENDIF
                0533       ENDIF
                0534 #endif /* ALLOW_ADDFLUID */
                0535 
                0536 #ifdef ALLOW_FRICTION_HEATING
                0537       IF ( addFrictionHeating ) THEN
                0538         IF ( fluidIsAir ) THEN
                0539 C         conversion from in-situ Temp to Pot.Temp
                0540           tmpFac = (atm_Po/rC(k))**atm_kappa
                0541 C         conversion from W/m^2/r_unit to K/s
                0542           tmpFac = (tmpFac/atm_Cp) * mass2rUnit
                0543         ELSE
                0544 C         conversion from W/m^2/r_unit to K/s
                0545           tmpFac = recip_Cp * mass2rUnit
                0546         ENDIF
5a705ed756 Jean*0547         DO j=0,sNy+1
                0548           DO i=0,sNx+1
e874fa47e5 Jean*0549             gT_arr(i,j) = gT_arr(i,j)
e24c9bfc82 Jean*0550      &         + frictionHeating(i,j,k,bi,bj)*tmpFac
e874fa47e5 Jean*0551      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0552           ENDDO
                0553         ENDDO
                0554       ENDIF
                0555 #endif /* ALLOW_FRICTION_HEATING */
                0556 
                0557       IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
                0558 C--   Compressible fluid: account for difference between moist and dry air
                0559 C     specific volume in Enthalpy equation (+ V.dP term), since only the
                0560 C     dry air part is accounted for in the (dry) Pot.Temp formulation.
                0561 C     Used centered averaging from interface to center (consistent with
                0562 C     conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
                0563 C     as for Theta_v in CALC_PHI_HYD
                0564 
                0565 C     conversion from in-situ Temp to Pot.Temp
                0566         tmpFac = (atm_Po/rC(k))**atm_kappa
                0567 C     conversion from W/kg to K/s
                0568         tmpFac = tmpFac/atm_Cp
                0569         km = k-1
                0570         kc = k
                0571         kp = k+1
                0572         IF ( k.EQ.1 ) THEN
5a705ed756 Jean*0573           DO j=0,sNy+1
                0574            DO i=0,sNx+1
e874fa47e5 Jean*0575             tmpVar(i,j) = 0.
                0576            ENDDO
                0577           ENDDO
                0578         ELSE
                0579           delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
                0580      &                   - (rC(kc)/atm_Po)**atm_kappa )
5a705ed756 Jean*0581           DO j=0,sNy+1
                0582            DO i=0,sNx+1
e874fa47e5 Jean*0583             tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
                0584      &                  *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
                0585      &                   + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
                0586      &                   )*maskC(i,j,km,bi,bj)*0.25 _d 0
                0587            ENDDO
                0588           ENDDO
                0589         ENDIF
                0590         IF ( k.LT.Nr ) THEN
                0591           delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
                0592      &                   - (rC(kp)/atm_Po)**atm_kappa )
5a705ed756 Jean*0593           DO j=0,sNy+1
                0594            DO i=0,sNx+1
e874fa47e5 Jean*0595             tmpVar(i,j) = tmpVar(i,j)
                0596      &                  + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
                0597      &                  *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
                0598      &                   + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
                0599      &                   )*maskC(i,j,kp,bi,bj)*0.25 _d 0
                0600            ENDDO
                0601           ENDDO
                0602         ENDIF
5a705ed756 Jean*0603         DO j=0,sNy+1
                0604           DO i=0,sNx+1
e874fa47e5 Jean*0605             gT_arr(i,j) = gT_arr(i,j)
                0606      &         + tmpVar(i,j)*tmpFac
                0607      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0608           ENDDO
                0609         ENDDO
                0610 #ifdef ALLOW_DIAGNOSTICS
                0611         IF ( useDiagnostics ) THEN
                0612 C     conversion to W/m^2
                0613           tmpFac = rUnit2mass
                0614           CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
5a705ed756 Jean*0615      &                     'MoistCor', kc, 1, 2, bi,bj,myThid )
e874fa47e5 Jean*0616         ENDIF
                0617 #endif /* ALLOW_DIAGNOSTICS */
                0618       ENDIF
                0619 
                0620 C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
                0621       IF ( k .EQ. kSurface ) THEN
5a705ed756 Jean*0622        DO j=0,sNy+1
                0623         DO i=0,sNx+1
e874fa47e5 Jean*0624           gT_arr(i,j) = gT_arr(i,j)
                0625      &      +surfaceForcingT(i,j,bi,bj)
                0626      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0627         ENDDO
                0628        ENDDO
                0629       ELSEIF ( kSurface.EQ.-1 ) THEN
5a705ed756 Jean*0630        DO j=0,sNy+1
                0631         DO i=0,sNx+1
e874fa47e5 Jean*0632          IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
                0633           gT_arr(i,j) = gT_arr(i,j)
                0634      &      +surfaceForcingT(i,j,bi,bj)
                0635      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0636          ENDIF
                0637         ENDDO
                0638        ENDDO
                0639       ENDIF
                0640 
                0641       IF (linFSConserveTr) THEN
5a705ed756 Jean*0642        DO j=0,sNy+1
                0643         DO i=0,sNx+1
0320e25227 Mart*0644          IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
                0645           gT_arr(i,j) = gT_arr(i,j)
e874fa47e5 Jean*0646      &        +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0320e25227 Mart*0647          ENDIF
e874fa47e5 Jean*0648         ENDDO
                0649        ENDDO
                0650       ENDIF
                0651 
90929f8806 Patr*0652 #ifdef ALLOW_GEOTHERMAL_FLUX
                0653       IF ( usingZCoords ) THEN
5a705ed756 Jean*0654        DO j=0,sNy+1
                0655         DO i=0,sNx+1
90929f8806 Patr*0656          IF ( k.EQ.kLowC(i,j,bi,bj) ) THEN
                0657           gT_arr(i,j)=gT_arr(i,j)
                0658      &      + geothermalFlux(i,j,bi,bj)
                0659      &        *recip_Cp*mass2rUnit
                0660      &        *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0661          ENDIF
                0662         ENDDO
                0663        ENDDO
0320e25227 Mart*0664       ELSEIF ( kSurface .EQ. Nr ) THEN
                0665 C     this is oceanic pressure coordinate case
                0666 C     where the flux at the bottom is applied as kSurfC
                0667        DO j=0,sNy+1
                0668         DO i=0,sNx+1
                0669          IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
                0670           gT_arr(i,j)=gT_arr(i,j)
                0671      &      + geothermalFlux(i,j,bi,bj)
                0672      &        *recip_Cp*mass2rUnit
                0673      &        *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0674          ENDIF
                0675         ENDDO
                0676        ENDDO
                0677       ELSE
44d3986245 Jean*0678 C-    Neither Z-Coords nor kSurface=Nr : not implemented
                0679        STOP 'ABNORMAL END: S/R APPLY_FORCING_T (geothermal-flux)'
90929f8806 Patr*0680       ENDIF
                0681 #endif /* ALLOW_GEOTHERMAL_FLUX */
                0682 
e874fa47e5 Jean*0683 #ifdef SHORTWAVE_HEATING
                0684 C Penetrating SW radiation
                0685 c     IF ( usePenetratingSW ) THEN
0320e25227 Mart*0686 C     Put depth of interface above & below current level in swfracb (1 & 2)
                0687 C     and SWFRAC returns fraction of light that crosses these interfaces ;
                0688 C     note: here km is for the mask below current level k
                0689        IF ( usingZCoords ) THEN
                0690         swfracb(1) = rF(k)
                0691         swfracb(2) = rF(k+1)
                0692         km = MIN(k+1,Nr)
                0693        ELSE
                0694 C     this is the oceanic pressure coordinate case
                0695         swfracb(1) = -rF(k+1)*recip_rhoConst*recip_gravity
                0696         swfracb(2) = -rF(k)*recip_rhoConst*recip_gravity
                0697         km = MAX(k-1,1)
                0698        ENDIF
e874fa47e5 Jean*0699        CALL SWFRAC(
0320e25227 Mart*0700      I             2, oneRL,
e874fa47e5 Jean*0701      U             swfracb,
0320e25227 Mart*0702      I             myTime, myIter, myThid )
                0703        IF ( k.EQ.km ) swfracb(2) = 0. _d 0
5a705ed756 Jean*0704        DO j=0,sNy+1
                0705         DO i=0,sNx+1
e874fa47e5 Jean*0706          gT_arr(i,j) = gT_arr(i,j)
0320e25227 Mart*0707      &        - Qsw(i,j,bi,bj)*( swfracb(1)*maskC(i,j,k,bi,bj)
                0708      &                         - swfracb(2)*maskC(i,j,km,bi,bj) )
                0709      &        *recip_Cp*mass2rUnit
                0710      &        *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e874fa47e5 Jean*0711         ENDDO
                0712        ENDDO
                0713 c     ENDIF
0320e25227 Mart*0714 #endif /* SHORTWAVE_HEATING */
e874fa47e5 Jean*0715 
                0716 #ifdef ALLOW_FRAZIL
                0717       IF ( useFRAZIL )
                0718      &     CALL FRAZIL_TENDENCY_APPLY_T(
                0719      U                 gT_arr,
                0720      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0721      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0722 #endif /* ALLOW_FRAZIL */
                0723 
                0724 #ifdef ALLOW_SHELFICE
                0725       IF ( useShelfIce )
                0726      &     CALL SHELFICE_FORCING_T(
                0727      U                   gT_arr,
                0728      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0729      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0730 #endif /* ALLOW_SHELFICE */
                0731 
                0732 #ifdef ALLOW_ICEFRONT
                0733       IF ( useICEFRONT )
                0734      &     CALL ICEFRONT_TENDENCY_APPLY_T(
                0735      U                   gT_arr,
2c160c3ab4 Jean*0736      I                   k, bi, bj, myTime, myIter, myThid )
e874fa47e5 Jean*0737 #endif /* ALLOW_ICEFRONT */
                0738 
                0739 #ifdef ALLOW_SALT_PLUME
                0740       IF ( useSALT_PLUME )
                0741      &     CALL SALT_PLUME_TENDENCY_APPLY_T(
                0742      U                     gT_arr,
                0743      I                     iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0744      I                     myTime, myIter, myThid )
e874fa47e5 Jean*0745 #endif /* ALLOW_SALT_PLUME */
                0746 
                0747 #ifdef ALLOW_RBCS
                0748       IF (useRBCS) THEN
                0749         CALL RBCS_ADD_TENDENCY(
                0750      U                 gT_arr,
                0751      I                 k, bi, bj, 1,
2c160c3ab4 Jean*0752      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0753       ENDIF
                0754 #endif /* ALLOW_RBCS */
                0755 
                0756 #ifdef ALLOW_OBCS
                0757       IF (useOBCS) THEN
                0758         CALL OBCS_SPONGE_T(
                0759      U                   gT_arr,
                0760      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0761      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0762       ENDIF
                0763 #endif /* ALLOW_OBCS */
                0764 
                0765 #ifdef ALLOW_BBL
                0766       IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
                0767      U                       gT_arr,
                0768      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0769      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0770 #endif /* ALLOW_BBL */
                0771 
                0772 #ifdef ALLOW_MYPACKAGE
                0773       IF ( useMYPACKAGE ) THEN
                0774         CALL MYPACKAGE_TENDENCY_APPLY_T(
                0775      U                 gT_arr,
                0776      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0777      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0778       ENDIF
                0779 #endif /* ALLOW_MYPACKAGE */
                0780 
                0781 #endif /* USE_OLD_EXTERNAL_FORCING */
                0782 
                0783       RETURN
                0784       END
                0785 
                0786 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0787 CBOP
                0788 C     !ROUTINE: APPLY_FORCING_S
                0789 C     !INTERFACE:
                0790       SUBROUTINE APPLY_FORCING_S(
                0791      U                     gS_arr,
                0792      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0793      I                     myTime, myIter, myThid )
                0794 C     !DESCRIPTION: \bv
                0795 C     *==========================================================*
                0796 C     | S/R APPLY_FORCING_S
                0797 C     | o Contains problem specific forcing for merid velocity.
                0798 C     *==========================================================*
                0799 C     | Adds terms to gS for forcing by external sources
                0800 C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
                0801 C     *==========================================================*
                0802 C     \ev
                0803 
                0804 C     !USES:
                0805       IMPLICIT NONE
                0806 C     == Global data ==
                0807 #include "SIZE.h"
                0808 #include "EEPARAMS.h"
                0809 #include "PARAMS.h"
                0810 #include "GRID.h"
                0811 #include "DYNVARS.h"
                0812 #include "FFIELDS.h"
                0813 #include "SURFACE.h"
                0814 
                0815 C     !INPUT/OUTPUT PARAMETERS:
                0816 C     gS_arr    :: the tendency array
                0817 C     iMin,iMax :: Working range of x-index for applying forcing.
                0818 C     jMin,jMax :: Working range of y-index for applying forcing.
                0819 C     k         :: Current vertical level index
                0820 C     bi,bj     :: Current tile indices
                0821 C     myTime    :: Current time in simulation
                0822 C     myIter    :: Current iteration number
                0823 C     myThid    :: my Thread Id number
                0824       _RL     gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0825       INTEGER iMin, iMax, jMin, jMax
                0826       INTEGER k, bi, bj
                0827       _RL     myTime
                0828       INTEGER myIter
                0829       INTEGER myThid
                0830 
                0831 C     !LOCAL VARIABLES:
                0832 C     i,j       :: Loop counters
                0833 C     kSurface  :: index of surface level
                0834       INTEGER i, j
5a705ed756 Jean*0835 #ifndef USE_OLD_EXTERNAL_FORCING
e874fa47e5 Jean*0836       INTEGER kSurface
5a705ed756 Jean*0837 #endif /* USE_OLD_EXTERNAL_FORCING */
e874fa47e5 Jean*0838 CEOP
                0839 
                0840 #ifdef USE_OLD_EXTERNAL_FORCING
                0841 
                0842       DO j=1-OLy,sNy+OLy
                0843         DO i=1-OLx,sNx+OLx
5a705ed756 Jean*0844           gS(i,j,k,bi,bj) = 0. _d 0
e874fa47e5 Jean*0845         ENDDO
                0846       ENDDO
                0847       CALL EXTERNAL_FORCING_S(
                0848      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0849      I              myTime, myThid )
                0850       DO j=1-OLy,sNy+OLy
                0851         DO i=1-OLx,sNx+OLx
5a705ed756 Jean*0852           gS_arr(i,j) = gS_arr(i,j) + gS(i,j,k,bi,bj)
e874fa47e5 Jean*0853         ENDDO
                0854       ENDDO
                0855 
                0856 #else  /* USE_OLD_EXTERNAL_FORCING */
                0857 
                0858       IF ( fluidIsAir ) THEN
                0859        kSurface = 0
                0860       ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
                0861        kSurface = -1
                0862       ELSEIF ( usingPCoords ) THEN
                0863        kSurface = Nr
                0864       ELSE
                0865        kSurface = 1
                0866       ENDIF
                0867 
5a705ed756 Jean*0868 C--   Note on loop range: For model dynamics, only needs to get correct
                0869 C     forcing (update gS_arr) in tile interior (1:sNx,1:sNy);
                0870 C     However, for some diagnostics, we may want to get valid forcing
                0871 C     extended over 1 point in tile halo region (0:sNx+1,0:sNy=1).
                0872 
e874fa47e5 Jean*0873 C--   Forcing term
                0874 #ifdef ALLOW_AIM
                0875       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
                0876      U                       gS_arr,
                0877      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0878      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0879 #endif /* ALLOW_AIM */
                0880 
                0881 #ifdef ALLOW_ATM_PHYS
                0882       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
                0883      U                       gS_arr,
                0884      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0885      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0886 #endif /* ALLOW_ATM_PHYS */
                0887 
                0888 #ifdef ALLOW_FIZHI
                0889       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
                0890      U                       gS_arr,
                0891      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0892      I                       myTime, myIter, myThid )
e874fa47e5 Jean*0893 #endif /* ALLOW_FIZHI */
                0894 
                0895 #ifdef ALLOW_ADDFLUID
                0896       IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
                0897        IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
                0898      &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
5a705ed756 Jean*0899          DO j=0,sNy+1
                0900           DO i=0,sNx+1
e874fa47e5 Jean*0901             gS_arr(i,j) = gS_arr(i,j)
                0902      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0903      &          *( salt_addMass - salt(i,j,k,bi,bj) )
                0904      &          *recip_rA(i,j,bi,bj)
                0905      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0906 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
                0907           ENDDO
                0908          ENDDO
                0909        ELSE
5a705ed756 Jean*0910          DO j=0,sNy+1
                0911           DO i=0,sNx+1
e874fa47e5 Jean*0912             gS_arr(i,j) = gS_arr(i,j)
                0913      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0914      &          *( salt_addMass - sRef(k) )
                0915      &          *recip_rA(i,j,bi,bj)
                0916      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0917 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
                0918           ENDDO
                0919          ENDDO
                0920        ENDIF
                0921       ENDIF
                0922 #endif /* ALLOW_ADDFLUID */
                0923 
                0924 C     Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
                0925       IF ( k .EQ. kSurface ) THEN
5a705ed756 Jean*0926        DO j=0,sNy+1
                0927         DO i=0,sNx+1
e874fa47e5 Jean*0928           gS_arr(i,j) = gS_arr(i,j)
                0929      &      +surfaceForcingS(i,j,bi,bj)
                0930      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0931         ENDDO
                0932        ENDDO
                0933       ELSEIF ( kSurface.EQ.-1 ) THEN
5a705ed756 Jean*0934        DO j=0,sNy+1
                0935         DO i=0,sNx+1
e874fa47e5 Jean*0936          IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
                0937           gS_arr(i,j) = gS_arr(i,j)
                0938      &      +surfaceForcingS(i,j,bi,bj)
                0939      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0940          ENDIF
                0941         ENDDO
                0942        ENDDO
                0943       ENDIF
                0944 
                0945       IF (linFSConserveTr) THEN
5a705ed756 Jean*0946        DO j=0,sNy+1
                0947         DO i=0,sNx+1
0320e25227 Mart*0948          IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
                0949           gS_arr(i,j) = gS_arr(i,j)
e874fa47e5 Jean*0950      &        +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
0320e25227 Mart*0951          ENDIF
e874fa47e5 Jean*0952         ENDDO
                0953        ENDDO
                0954       ENDIF
                0955 
                0956 #ifdef ALLOW_SHELFICE
                0957       IF ( useShelfIce )
                0958      &     CALL SHELFICE_FORCING_S(
                0959      U                   gS_arr,
                0960      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0961      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0962 #endif /* ALLOW_SHELFICE */
                0963 
                0964 #ifdef ALLOW_ICEFRONT
                0965       IF ( useICEFRONT )
                0966      &     CALL ICEFRONT_TENDENCY_APPLY_S(
                0967      U                   gS_arr,
2c160c3ab4 Jean*0968      I                   k, bi, bj, myTime, myIter, myThid )
e874fa47e5 Jean*0969 #endif /* ALLOW_ICEFRONT */
                0970 
                0971 #ifdef ALLOW_SALT_PLUME
                0972       IF ( useSALT_PLUME )
                0973      &     CALL SALT_PLUME_TENDENCY_APPLY_S(
                0974      U                     gS_arr,
                0975      I                     iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0976      I                     myTime, myIter, myThid )
e874fa47e5 Jean*0977 #endif /* ALLOW_SALT_PLUME */
                0978 
                0979 #ifdef ALLOW_RBCS
                0980       IF (useRBCS) THEN
                0981         CALL RBCS_ADD_TENDENCY(
                0982      U                 gS_arr,
                0983      I                 k, bi, bj, 2,
2c160c3ab4 Jean*0984      I                 myTime, myIter, myThid )
e874fa47e5 Jean*0985       ENDIF
                0986 #endif /* ALLOW_RBCS */
                0987 
                0988 #ifdef ALLOW_OBCS
                0989       IF (useOBCS) THEN
                0990         CALL OBCS_SPONGE_S(
                0991      U                   gS_arr,
                0992      I                   iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*0993      I                   myTime, myIter, myThid )
e874fa47e5 Jean*0994       ENDIF
                0995 #endif /* ALLOW_OBCS */
                0996 
                0997 #ifdef ALLOW_BBL
                0998       IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
                0999      U                       gS_arr,
                1000      I                       iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*1001      I                       myTime, myIter, myThid )
e874fa47e5 Jean*1002 #endif /* ALLOW_BBL */
                1003 
                1004 #ifdef ALLOW_MYPACKAGE
                1005       IF ( useMYPACKAGE ) THEN
                1006         CALL MYPACKAGE_TENDENCY_APPLY_S(
                1007      U                 gS_arr,
                1008      I                 iMin,iMax,jMin,jMax, k, bi,bj,
2c160c3ab4 Jean*1009      I                 myTime, myIter, myThid )
e874fa47e5 Jean*1010       ENDIF
                1011 #endif /* ALLOW_MYPACKAGE */
                1012 
                1013 #endif /* USE_OLD_EXTERNAL_FORCING */
                1014 
                1015       RETURN
                1016       END