tutorial_rotating_tank/code/apply_forcing.F

34ee56fea5 Jean*0001 #include "PACKAGES_CONFIG.h"
19780451ef Andr*0002 #include "CPP_OPTIONS.h"
e384178ec8 Jean*0003
df4a605949 Jean*0004 C--  File apply_forcing.F:
e384178ec8 Jean*0005 C--   Contents
df4a605949 Jean*0006 C--   o APPLY_FORCING_U
                0007 C--   o APPLY_FORCING_V
                0008 C--   o APPLY_FORCING_T
                0009 C--   o APPLY_FORCING_S
19780451ef Andr*0010
e384178ec8 Jean*0011 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
19780451ef Andr*0012 CBOP
df4a605949 Jean*0013 C     !ROUTINE: APPLY_FORCING_U
19780451ef Andr*0014 C     !INTERFACE:
df4a605949 Jean*0015       SUBROUTINE APPLY_FORCING_U(
                0016      U                     gU_arr,
                0017      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0018      I                     myTime, myIter, myThid )
19780451ef Andr*0019 C     !DESCRIPTION: \bv
                0020 C     *==========================================================*
df4a605949 Jean*0021 C     | S/R APPLY_FORCING_U
34ee56fea5 Jean*0022 C     | o Contains problem specific forcing for zonal velocity.
19780451ef Andr*0023 C     *==========================================================*
34ee56fea5 Jean*0024 C     | Adds terms to gU for forcing by external sources
                0025 C     | e.g. wind stress, bottom friction etc ...
19780451ef Andr*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:
df4a605949 Jean*0040 C     gU_arr    :: the tendency array
34ee56fea5 Jean*0041 C     iMin,iMax :: Working range of x-index for applying forcing.
                0042 C     jMin,jMax :: Working range of y-index for applying forcing.
df4a605949 Jean*0043 C     k         :: Current vertical level index
34ee56fea5 Jean*0044 C     bi,bj     :: Current tile indices
                0045 C     myTime    :: Current time in simulation
df4a605949 Jean*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
19780451ef Andr*0053       INTEGER myThid
                0054
                0055 C     !LOCAL VARIABLES:
34ee56fea5 Jean*0056 C     i,j       :: Loop counters
52d51b46b8 Jean*0057 C     kSurface  :: index of surface level
34ee56fea5 Jean*0058       INTEGER i, j
df4a605949 Jean*0059 #ifdef USE_OLD_EXTERNAL_FORCING
                0060       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL     tmpVar
                0062 #else
34ee56fea5 Jean*0063       INTEGER kSurface
df4a605949 Jean*0064 #endif
19780451ef Andr*0065 CEOP
                0066
df4a605949 Jean*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 )
                0077       DO j=1-OLy,sNy+OLy
                0078         DO i=1-OLx,sNx+OLx
                0079           tmpVar = gU(i,j,k,bi,bj) - locVar(i,j)
                0080           gU(i,j,k,bi,bj) = locVar(i,j)
                0081           gU_arr(i,j) = gU_arr(i,j) + tmpVar
                0082         ENDDO
                0083       ENDDO
                0084
                0085 #else  /* USE_OLD_EXTERNAL_FORCING */
                0086
34ee56fea5 Jean*0087       IF ( fluidIsAir ) THEN
                0088        kSurface = 0
                0089       ELSEIF ( usingPCoords ) THEN
                0090        kSurface = Nr
                0091       ELSE
                0092        kSurface = 1
                0093       ENDIF
                0094
19780451ef Andr*0095 C--   Forcing term
34ee56fea5 Jean*0096 #ifdef ALLOW_AIM
                0097       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_U(
df4a605949 Jean*0098      U                       gU_arr,
                0099      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0100      I                       myTime, 0, myThid )
34ee56fea5 Jean*0101 #endif /* ALLOW_AIM */
                0102
e384178ec8 Jean*0103 #ifdef ALLOW_ATM_PHYS
                0104       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_U(
df4a605949 Jean*0105      U                       gU_arr,
                0106      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0107      I                       myTime, 0, myThid )
                0108 #endif /* ALLOW_ATM_PHYS */
                0109
34ee56fea5 Jean*0110 #ifdef ALLOW_FIZHI
                0111       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_U(
df4a605949 Jean*0112      U                       gU_arr,
                0113      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0114      I                       myTime, 0, myThid )
34ee56fea5 Jean*0115 #endif /* ALLOW_FIZHI */
                0116
e384178ec8 Jean*0117 C     Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
df4a605949 Jean*0118       IF ( k .EQ. kSurface ) THEN
34ee56fea5 Jean*0119 c      DO j=1,sNy
                0120 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNy+1]
                0121        DO j=0,sNy+1
                0122         DO i=1,sNx+1
df4a605949 Jean*0123           gU_arr(i,j) = gU_arr(i,j)
e384178ec8 Jean*0124      &      +foFacMom*surfaceForcingU(i,j,bi,bj)
df4a605949 Jean*0125      &      *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
e384178ec8 Jean*0126         ENDDO
                0127        ENDDO
                0128       ELSEIF ( kSurface.EQ.-1 ) THEN
                0129        DO j=0,sNy+1
                0130         DO i=1,sNx+1
df4a605949 Jean*0131          IF ( kSurfW(i,j,bi,bj).EQ.k ) THEN
                0132           gU_arr(i,j) = gU_arr(i,j)
e384178ec8 Jean*0133      &      +foFacMom*surfaceForcingU(i,j,bi,bj)
df4a605949 Jean*0134      &      *recip_drF(k)*_recip_hFacW(i,j,k,bi,bj)
e384178ec8 Jean*0135          ENDIF
19780451ef Andr*0136         ENDDO
                0137        ENDDO
                0138       ENDIF
                0139
e384178ec8 Jean*0140 #ifdef ALLOW_EDDYPSI
                0141          CALL TAUEDDY_TENDENCY_APPLY_U(
df4a605949 Jean*0142      U                 gU_arr,
                0143      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0144      I                 myTime, 0, myThid )
34ee56fea5 Jean*0145 #endif
                0146
e384178ec8 Jean*0147 #ifdef ALLOW_RBCS
                0148       IF (useRBCS) THEN
                0149         CALL RBCS_ADD_TENDENCY(
df4a605949 Jean*0150      U                 gU_arr,
                0151      I                 k, bi, bj, -1,
e384178ec8 Jean*0152      I                 myTime, 0, myThid )
                0153
                0154       ENDIF
                0155 #endif /* ALLOW_RBCS */
                0156
34ee56fea5 Jean*0157 #ifdef ALLOW_OBCS
19780451ef Andr*0158       IF (useOBCS) THEN
e384178ec8 Jean*0159         CALL OBCS_SPONGE_U(
df4a605949 Jean*0160      U                   gU_arr,
                0161      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0162      I                   myTime, 0, myThid )
19780451ef Andr*0163       ENDIF
e384178ec8 Jean*0164 #endif /* ALLOW_OBCS */
                0165
                0166 #ifdef ALLOW_MYPACKAGE
                0167       IF ( useMYPACKAGE ) THEN
                0168         CALL MYPACKAGE_TENDENCY_APPLY_U(
df4a605949 Jean*0169      U                 gU_arr,
                0170      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0171      I                 myTime, 0, myThid )
                0172       ENDIF
                0173 #endif /* ALLOW_MYPACKAGE */
19780451ef Andr*0174
df4a605949 Jean*0175 #endif /* USE_OLD_EXTERNAL_FORCING */
                0176
19780451ef Andr*0177       RETURN
                0178       END
34ee56fea5 Jean*0179
                0180 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
19780451ef Andr*0181 CBOP
df4a605949 Jean*0182 C     !ROUTINE: APPLY_FORCING_V
19780451ef Andr*0183 C     !INTERFACE:
df4a605949 Jean*0184       SUBROUTINE APPLY_FORCING_V(
                0185      U                     gV_arr,
                0186      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0187      I                     myTime, myIter, myThid )
19780451ef Andr*0188 C     !DESCRIPTION: \bv
                0189 C     *==========================================================*
df4a605949 Jean*0190 C     | S/R APPLY_FORCING_V
34ee56fea5 Jean*0191 C     | o Contains problem specific forcing for merid velocity.
19780451ef Andr*0192 C     *==========================================================*
34ee56fea5 Jean*0193 C     | Adds terms to gV for forcing by external sources
                0194 C     | e.g. wind stress, bottom friction etc ...
19780451ef Andr*0195 C     *==========================================================*
                0196 C     \ev
                0197
                0198 C     !USES:
                0199       IMPLICIT NONE
                0200 C     == Global data ==
                0201 #include "SIZE.h"
                0202 #include "EEPARAMS.h"
                0203 #include "PARAMS.h"
                0204 #include "GRID.h"
                0205 #include "DYNVARS.h"
                0206 #include "FFIELDS.h"
                0207
                0208 C     !INPUT/OUTPUT PARAMETERS:
df4a605949 Jean*0209 C     gV_arr    :: the tendency array
34ee56fea5 Jean*0210 C     iMin,iMax :: Working range of x-index for applying forcing.
                0211 C     jMin,jMax :: Working range of y-index for applying forcing.
df4a605949 Jean*0212 C     k         :: Current vertical level index
34ee56fea5 Jean*0213 C     bi,bj     :: Current tile indices
                0214 C     myTime    :: Current time in simulation
df4a605949 Jean*0215 C     myIter    :: Current iteration number
                0216 C     myThid    :: my Thread Id number
                0217       _RL     gV_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0218       INTEGER iMin, iMax, jMin, jMax
                0219       INTEGER k, bi, bj
                0220       _RL     myTime
                0221       INTEGER myIter
19780451ef Andr*0222       INTEGER myThid
                0223
                0224 C     !LOCAL VARIABLES:
34ee56fea5 Jean*0225 C     i,j       :: Loop counters
52d51b46b8 Jean*0226 C     kSurface  :: index of surface level
34ee56fea5 Jean*0227       INTEGER i, j
df4a605949 Jean*0228 #ifdef USE_OLD_EXTERNAL_FORCING
                0229       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0230       _RL     tmpVar
                0231 #else
34ee56fea5 Jean*0232       INTEGER kSurface
df4a605949 Jean*0233 #endif
19780451ef Andr*0234 CEOP
                0235
df4a605949 Jean*0236 #ifdef USE_OLD_EXTERNAL_FORCING
                0237
                0238       DO j=1-OLy,sNy+OLy
                0239         DO i=1-OLx,sNx+OLx
                0240           locVar(i,j) = gV(i,j,k,bi,bj)
                0241         ENDDO
                0242       ENDDO
                0243       CALL EXTERNAL_FORCING_V(
                0244      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0245      I              myTime, myThid )
                0246       DO j=1-OLy,sNy+OLy
                0247         DO i=1-OLx,sNx+OLx
                0248           tmpVar = gV(i,j,k,bi,bj) - locVar(i,j)
                0249           gV(i,j,k,bi,bj) = locVar(i,j)
                0250           gV_arr(i,j) = gV_arr(i,j) + tmpVar
                0251         ENDDO
                0252       ENDDO
                0253
                0254 #else  /* USE_OLD_EXTERNAL_FORCING */
                0255
34ee56fea5 Jean*0256       IF ( fluidIsAir ) THEN
                0257        kSurface = 0
                0258       ELSEIF ( usingPCoords ) THEN
                0259        kSurface = Nr
                0260       ELSE
                0261        kSurface = 1
                0262       ENDIF
                0263
19780451ef Andr*0264 C--   Forcing term
34ee56fea5 Jean*0265 #ifdef ALLOW_AIM
                0266       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_V(
df4a605949 Jean*0267      U                       gV_arr,
                0268      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0269      I                       myTime, 0, myThid )
34ee56fea5 Jean*0270 #endif /* ALLOW_AIM */
                0271
e384178ec8 Jean*0272 #ifdef ALLOW_ATM_PHYS
                0273       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_V(
df4a605949 Jean*0274      U                       gV_arr,
                0275      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0276      I                       myTime, 0, myThid )
                0277 #endif /* ALLOW_ATM_PHYS */
                0278
34ee56fea5 Jean*0279 #ifdef ALLOW_FIZHI
                0280       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_V(
df4a605949 Jean*0281      U                       gV_arr,
                0282      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0283      I                       myTime, 0, myThid )
34ee56fea5 Jean*0284 #endif /* ALLOW_FIZHI */
                0285
e384178ec8 Jean*0286 C     Ocean: Add momentum surface forcing (e.g., wind-stress) in surface level
df4a605949 Jean*0287       IF ( k .EQ. kSurface ) THEN
34ee56fea5 Jean*0288        DO j=1,sNy+1
                0289 c       DO i=1,sNx
                0290 C-jmc: Without CD-scheme, this is OK ; but with CD-scheme, needs to cover [0:sNx+1]
                0291         DO i=0,sNx+1
df4a605949 Jean*0292           gV_arr(i,j) = gV_arr(i,j)
e384178ec8 Jean*0293      &      +foFacMom*surfaceForcingV(i,j,bi,bj)
df4a605949 Jean*0294      &      *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
e384178ec8 Jean*0295         ENDDO
                0296        ENDDO
                0297       ELSEIF ( kSurface.EQ.-1 ) THEN
                0298        DO j=1,sNy+1
                0299         DO i=0,sNx+1
df4a605949 Jean*0300          IF ( kSurfS(i,j,bi,bj).EQ.k ) THEN
                0301           gV_arr(i,j) = gV_arr(i,j)
e384178ec8 Jean*0302      &      +foFacMom*surfaceForcingV(i,j,bi,bj)
df4a605949 Jean*0303      &      *recip_drF(k)*_recip_hFacS(i,j,k,bi,bj)
e384178ec8 Jean*0304          ENDIF
19780451ef Andr*0305         ENDDO
                0306        ENDDO
                0307       ENDIF
                0308
e384178ec8 Jean*0309 #ifdef ALLOW_EDDYPSI
                0310          CALL TAUEDDY_TENDENCY_APPLY_V(
df4a605949 Jean*0311      U                 gV_arr,
                0312      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0313      I                 myTime, 0, myThid )
34ee56fea5 Jean*0314 #endif
                0315
e384178ec8 Jean*0316 #ifdef ALLOW_RBCS
                0317       IF (useRBCS) THEN
                0318         CALL RBCS_ADD_TENDENCY(
df4a605949 Jean*0319      U                 gV_arr,
                0320      I                 k, bi, bj, -2,
e384178ec8 Jean*0321      I                 myTime, 0, myThid )
                0322       ENDIF
                0323 #endif /* ALLOW_RBCS */
                0324
34ee56fea5 Jean*0325 #ifdef ALLOW_OBCS
19780451ef Andr*0326       IF (useOBCS) THEN
e384178ec8 Jean*0327         CALL OBCS_SPONGE_V(
df4a605949 Jean*0328      U                   gV_arr,
                0329      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0330      I                   myTime, 0, myThid )
19780451ef Andr*0331       ENDIF
e384178ec8 Jean*0332 #endif /* ALLOW_OBCS */
                0333
                0334 #ifdef ALLOW_MYPACKAGE
                0335       IF ( useMYPACKAGE ) THEN
                0336         CALL MYPACKAGE_TENDENCY_APPLY_V(
df4a605949 Jean*0337      U                 gV_arr,
                0338      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0339      I                 myTime, 0, myThid )
                0340       ENDIF
                0341 #endif /* ALLOW_MYPACKAGE */
19780451ef Andr*0342
df4a605949 Jean*0343 #endif /* USE_OLD_EXTERNAL_FORCING */
                0344
19780451ef Andr*0345       RETURN
                0346       END
34ee56fea5 Jean*0347
                0348 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
19780451ef Andr*0349 CBOP
df4a605949 Jean*0350 C     !ROUTINE: APPLY_FORCING_T
19780451ef Andr*0351 C     !INTERFACE:
df4a605949 Jean*0352       SUBROUTINE APPLY_FORCING_T(
                0353      U                     gT_arr,
                0354      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0355      I                     myTime, myIter, myThid )
19780451ef Andr*0356 C     !DESCRIPTION: \bv
                0357 C     *==========================================================*
df4a605949 Jean*0358 C     | S/R APPLY_FORCING_T
34ee56fea5 Jean*0359 C     | o Contains problem specific forcing for temperature.
19780451ef Andr*0360 C     *==========================================================*
34ee56fea5 Jean*0361 C     | Adds terms to gT for forcing by external sources
                0362 C     | e.g. heat flux, climatalogical relaxation, etc ...
19780451ef Andr*0363 C     *==========================================================*
                0364 C     \ev
                0365
                0366 C     !USES:
                0367       IMPLICIT NONE
                0368 C     == Global data ==
                0369 #include "SIZE.h"
                0370 #include "EEPARAMS.h"
                0371 #include "PARAMS.h"
                0372 #include "GRID.h"
                0373 #include "DYNVARS.h"
                0374 #include "FFIELDS.h"
e384178ec8 Jean*0375 #include "SURFACE.h"
19780451ef Andr*0376
                0377 C     !INPUT/OUTPUT PARAMETERS:
df4a605949 Jean*0378 C     gT_arr    :: the tendency array
34ee56fea5 Jean*0379 C     iMin,iMax :: Working range of x-index for applying forcing.
                0380 C     jMin,jMax :: Working range of y-index for applying forcing.
df4a605949 Jean*0381 C     k         :: Current vertical level index
34ee56fea5 Jean*0382 C     bi,bj     :: Current tile indices
                0383 C     myTime    :: Current time in simulation
df4a605949 Jean*0384 C     myIter    :: Current iteration number
                0385 C     myThid    :: my Thread Id number
                0386       _RL     gT_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0387       INTEGER iMin, iMax, jMin, jMax
                0388       INTEGER k, bi, bj
                0389       _RL     myTime
                0390       INTEGER myIter
19780451ef Andr*0391       INTEGER myThid
                0392
                0393 C     !LOCAL VARIABLES:
34ee56fea5 Jean*0394 C     i,j       :: Loop counters
52d51b46b8 Jean*0395 C     kSurface  :: index of surface level
34ee56fea5 Jean*0396       INTEGER i, j
df4a605949 Jean*0397 #ifdef USE_OLD_EXTERNAL_FORCING
                0398       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0399       _RL     tmpVar
                0400 #else
34ee56fea5 Jean*0401       INTEGER kSurface
e384178ec8 Jean*0402       INTEGER km, kc, kp
                0403       _RL tmpVar(1:sNx,1:sNy)
                0404       _RL tmpFac, delPI
52d51b46b8 Jean*0405       _RL recip_Cp
34ee56fea5 Jean*0406 CEOP
e384178ec8 Jean*0407 #ifdef SHORTWAVE_HEATING
                0408       _RL minusone
                0409       PARAMETER (minusOne=-1.)
                0410       _RL swfracb(2)
                0411       INTEGER kp1
                0412 #endif
                0413 C---- Experiment specific declaration: starts here -------------------
19780451ef Andr*0414 C     iG, jG             :: Global index temps.
                0415 C     hC, hW, hE, hN, hS :: Fractional vertical distance open to fluid temps.
                0416 C     dFlux[WENS]        :: Diffusive flux normal to each cell face.
                0417 C     faceArea           :: Temp. for holding area normal to tempurature gradient.
                0418       INTEGER iG, jG
                0419       _RL hC, hW, hE, hN, hS
                0420       _RL dFluxW, dFluxE, dFluxN, dFluxS
                0421       _RL faceArea
                0422
                0423 C--   Forcing term
2e53480479 Jean*0424 C     Add term which represents the diffusive flux from a circular cylinder of temperature tCylIm in the
19780451ef Andr*0425 C     interior of the simulation domain. Result is a tendency which is determined from the finite
2e53480479 Jean*0426 C     volume formulated divergence of the diffusive heat flux due to the local cylinder
19780451ef Andr*0427 C     temperature, fluid temperature difference.
                0428 C     kDiffCyl :: Diffusion coefficient
963e35159b Andr*0429 C     tCylIn     :: Temperature of the inner boundary of the cylinder
                0430 C     tCylOut     :: Temperature of the outer boundary cylinder
19780451ef Andr*0431 C     iGSource :: Index space I (global) coordinate for source center.
                0432 C     jGSource :: Index space J (global) coordinate for source center.
                0433 C     rSource  :: Extent of the source term region. Loop will skip checking points outside
                0434 C              :: this region. Within this region the source heating will be added
                0435 C              :: to any points that are at a land - fluid boundary. rSource is in grid
                0436 C              :: points, so that points checked are ophi(iGlobal,jGlobal) such that
                0437 C              :: iGlobal=iGsource +/- rSource, jGlobal = jGsource +/- rSource.
                0438 C              :: rSource, iGSource and jGSource only need to define a quadrilateral that
                0439 C              :: includes the cylinder and no other land, they do not need to be exact.
e384178ec8 Jean*0440 C afe:
                0441 C  the below appears to be an attempt to make the heat flux somewhat general in regards
                0442 C  to bathymetry, but jmc pointed out some ways it could be better.  this is not
                0443 C  an issue at this point (July 04) since all experiments are being done with straight-
                0444 C  sided tank and cyclinder walls.
                0445 C  some todo items:
                0446 C  * add terms to top and bottom -- probably critical!!!
                0447 C  * make tCyl more flexible -- maybe have as separate inner and outer variables
                0448 C       or, eventually, a forcing field
                0449 C  * think about possible problems with differing heat diffusion rates in wall materials
                0450 C    (plexiglass, air, water in the compound tank case)
19780451ef Andr*0451       _RL kDiffCyl
963e35159b Andr*0452       _RL tCyl
19780451ef Andr*0453       INTEGER rSource
                0454       INTEGER iGSource
                0455       INTEGER jGSource
e384178ec8 Jean*0456 C---- Experiment specific declaration:  ends  here -------------------
df4a605949 Jean*0457 #endif
                0458
                0459 #ifdef USE_OLD_EXTERNAL_FORCING
                0460
                0461       DO j=1-OLy,sNy+OLy
                0462         DO i=1-OLx,sNx+OLx
                0463           locVar(i,j) = gT(i,j,k,bi,bj)
                0464         ENDDO
                0465       ENDDO
                0466       CALL EXTERNAL_FORCING_T(
                0467      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0468      I              myTime, myThid )
                0469       DO j=1-OLy,sNy+OLy
                0470         DO i=1-OLx,sNx+OLx
                0471           tmpVar = gT(i,j,k,bi,bj) - locVar(i,j)
                0472           gT(i,j,k,bi,bj) = locVar(i,j)
                0473           gT_arr(i,j) = gT_arr(i,j) + tmpVar
                0474         ENDDO
                0475       ENDDO
                0476
                0477 #else  /* USE_OLD_EXTERNAL_FORCING */
34ee56fea5 Jean*0478
                0479       IF ( fluidIsAir ) THEN
                0480        kSurface = 0
e384178ec8 Jean*0481       ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
                0482        kSurface = -1
34ee56fea5 Jean*0483       ELSEIF ( usingPCoords ) THEN
                0484        kSurface = Nr
                0485       ELSE
                0486        kSurface = 1
                0487       ENDIF
52d51b46b8 Jean*0488       recip_Cp = 1. _d 0 / HeatCapacity_Cp
34ee56fea5 Jean*0489
                0490 C--   Forcing term
                0491 #ifdef ALLOW_AIM
                0492       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_T(
df4a605949 Jean*0493      U                       gT_arr,
                0494      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0495      I                       myTime, 0, myThid )
34ee56fea5 Jean*0496 #endif /* ALLOW_AIM */
                0497
e384178ec8 Jean*0498 #ifdef ALLOW_ATM_PHYS
                0499       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_T(
df4a605949 Jean*0500      U                       gT_arr,
                0501      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0502      I                       myTime, 0, myThid )
                0503 #endif /* ALLOW_ATM_PHYS */
                0504
34ee56fea5 Jean*0505 #ifdef ALLOW_FIZHI
                0506       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_T(
df4a605949 Jean*0507      U                       gT_arr,
                0508      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0509      I                       myTime, 0, myThid )
34ee56fea5 Jean*0510 #endif /* ALLOW_FIZHI */
                0511
e384178ec8 Jean*0512 #ifdef ALLOW_ADDFLUID
                0513       IF ( selectAddFluid.NE.0 .AND. temp_addMass.NE.UNSET_RL ) THEN
                0514        IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
                0515      &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
                0516          DO j=1,sNy
                0517           DO i=1,sNx
df4a605949 Jean*0518             gT_arr(i,j) = gT_arr(i,j)
                0519      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0520      &          *( temp_addMass - theta(i,j,k,bi,bj) )
e384178ec8 Jean*0521      &          *recip_rA(i,j,bi,bj)
df4a605949 Jean*0522      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0523 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
e384178ec8 Jean*0524           ENDDO
                0525          ENDDO
                0526        ELSE
                0527          DO j=1,sNy
                0528           DO i=1,sNx
df4a605949 Jean*0529             gT_arr(i,j) = gT_arr(i,j)
                0530      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0531      &          *( temp_addMass - tRef(k) )
e384178ec8 Jean*0532      &          *recip_rA(i,j,bi,bj)
df4a605949 Jean*0533      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0534 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
e384178ec8 Jean*0535           ENDDO
                0536          ENDDO
                0537        ENDIF
                0538       ENDIF
                0539 #endif /* ALLOW_ADDFLUID */
                0540
                0541 #ifdef ALLOW_FRICTION_HEATING
                0542       IF ( addFrictionHeating ) THEN
                0543         IF ( fluidIsAir ) THEN
                0544 C         conversion from in-situ Temp to Pot.Temp
df4a605949 Jean*0545           tmpFac = (atm_Po/rC(k))**atm_kappa
e384178ec8 Jean*0546 C         conversion from W/m^2/r_unit to K/s
                0547           tmpFac = (tmpFac/atm_Cp) * mass2rUnit
                0548         ELSE
                0549 C         conversion from W/m^2/r_unit to K/s
                0550           tmpFac = recip_Cp * mass2rUnit
                0551         ENDIF
                0552         DO j=1,sNy
                0553           DO i=1,sNx
df4a605949 Jean*0554             gT_arr(i,j) = gT_arr(i,j)
                0555      &         + frictionHeating(i,j,k,bi,bj)
e384178ec8 Jean*0556      &          *tmpFac*recip_rA(i,j,bi,bj)
df4a605949 Jean*0557      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0558           ENDDO
                0559         ENDDO
                0560       ENDIF
                0561 #endif /* ALLOW_FRICTION_HEATING */
                0562
                0563       IF ( fluidIsAir .AND. atm_Rq.NE.zeroRL .AND. Nr.NE.1 ) THEN
                0564 C--   Compressible fluid: account for difference between moist and dry air
                0565 C     specific volume in Enthalpy equation (+ V.dP term), since only the
                0566 C     dry air part is accounted for in the (dry) Pot.Temp formulation.
                0567 C     Used centered averaging from interface to center (consistent with
                0568 C     conversion term in KE eq) and same discretisation ( [T*Q]_bar_k )
                0569 C     as for Theta_v in CALC_PHI_HYD
                0570
                0571 C     conversion from in-situ Temp to Pot.Temp
df4a605949 Jean*0572         tmpFac = (atm_Po/rC(k))**atm_kappa
e384178ec8 Jean*0573 C     conversion from W/kg to K/s
                0574         tmpFac = tmpFac/atm_Cp
df4a605949 Jean*0575         km = k-1
                0576         kc = k
                0577         kp = k+1
                0578         IF ( k.EQ.1 ) THEN
e384178ec8 Jean*0579           DO j=1,sNy
                0580            DO i=1,sNx
                0581             tmpVar(i,j) = 0.
                0582            ENDDO
                0583           ENDDO
                0584         ELSE
                0585           delPI = atm_Cp*( (rC(km)/atm_Po)**atm_kappa
                0586      &                   - (rC(kc)/atm_Po)**atm_kappa )
                0587           DO j=1,sNy
                0588            DO i=1,sNx
                0589             tmpVar(i,j) = wVel(i,j,kc,bi,bj)*delPI*atm_Rq
                0590      &                  *( theta(i,j,km,bi,bj)*salt(i,j,km,bi,bj)
                0591      &                   + theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
                0592      &                   )*maskC(i,j,km,bi,bj)*0.25 _d 0
                0593            ENDDO
                0594           ENDDO
                0595         ENDIF
df4a605949 Jean*0596         IF ( k.LT.Nr ) THEN
e384178ec8 Jean*0597           delPI = atm_Cp*( (rC(kc)/atm_Po)**atm_kappa
                0598      &                   - (rC(kp)/atm_Po)**atm_kappa )
                0599           DO j=1,sNy
                0600            DO i=1,sNx
                0601             tmpVar(i,j) = tmpVar(i,j)
                0602      &                  + wVel(i,j,kp,bi,bj)*delPI*atm_Rq
                0603      &                  *( theta(i,j,kc,bi,bj)*salt(i,j,kc,bi,bj)
                0604      &                   + theta(i,j,kp,bi,bj)*salt(i,j,kp,bi,bj)
                0605      &                   )*maskC(i,j,kp,bi,bj)*0.25 _d 0
                0606            ENDDO
                0607           ENDDO
                0608         ENDIF
                0609         DO j=1,sNy
                0610           DO i=1,sNx
df4a605949 Jean*0611             gT_arr(i,j) = gT_arr(i,j)
e384178ec8 Jean*0612      &         + tmpVar(i,j)*tmpFac
df4a605949 Jean*0613      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0614           ENDDO
                0615         ENDDO
                0616 #ifdef ALLOW_DIAGNOSTICS
                0617         IF ( useDiagnostics ) THEN
                0618 C     conversion to W/m^2
                0619           tmpFac = rUnit2mass
                0620           CALL DIAGNOSTICS_SCALE_FILL( tmpVar, tmpFac, 1,
                0621      &                     'MoistCor', kc, 1, 3, bi,bj,myThid )
                0622         ENDIF
                0623 #endif /* ALLOW_DIAGNOSTICS */
                0624       ENDIF
                0625
                0626 C     Ocean: Add temperature surface forcing (e.g., heat-flux) in surface level
df4a605949 Jean*0627       IF ( k .EQ. kSurface ) THEN
34ee56fea5 Jean*0628        DO j=1,sNy
                0629         DO i=1,sNx
df4a605949 Jean*0630           gT_arr(i,j) = gT_arr(i,j)
e384178ec8 Jean*0631      &      +surfaceForcingT(i,j,bi,bj)
df4a605949 Jean*0632      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0633         ENDDO
                0634        ENDDO
                0635       ELSEIF ( kSurface.EQ.-1 ) THEN
                0636        DO j=1,sNy
                0637         DO i=1,sNx
df4a605949 Jean*0638          IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
                0639           gT_arr(i,j) = gT_arr(i,j)
e384178ec8 Jean*0640      &      +surfaceForcingT(i,j,bi,bj)
df4a605949 Jean*0641      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0642          ENDIF
                0643         ENDDO
                0644        ENDDO
                0645       ENDIF
                0646
                0647       IF (linFSConserveTr) THEN
                0648        DO j=1,sNy
                0649         DO i=1,sNx
df4a605949 Jean*0650           IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
                0651             gT_arr(i,j) = gT_arr(i,j)
                0652      &        +TsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0653           ENDIF
34ee56fea5 Jean*0654         ENDDO
                0655        ENDDO
                0656       ENDIF
                0657
                0658 #ifdef SHORTWAVE_HEATING
                0659 C Penetrating SW radiation
                0660 c     IF ( usePenetratingSW ) THEN
df4a605949 Jean*0661        swfracb(1)=abs(rF(k))
                0662        swfracb(2)=abs(rF(k+1))
34ee56fea5 Jean*0663        CALL SWFRAC(
e384178ec8 Jean*0664      I             2, minusOne,
e872d3901a Jean*0665      U             swfracb,
                0666      I             myTime, 1, myThid )
df4a605949 Jean*0667        kp1 = k+1
                0668        IF (k.EQ.Nr) THEN
                0669         kp1 = k
34ee56fea5 Jean*0670         swfracb(2)=0. _d 0
                0671        ENDIF
                0672        DO j=1,sNy
                0673         DO i=1,sNx
df4a605949 Jean*0674          gT_arr(i,j) = gT_arr(i,j)
                0675      &   -Qsw(i,j,bi,bj)*(swfracb(1)*maskC(i,j,k,bi,bj)
34ee56fea5 Jean*0676      &                   -swfracb(2)*maskC(i,j,kp1, bi,bj))
e384178ec8 Jean*0677      &    *recip_Cp*mass2rUnit
df4a605949 Jean*0678      &    *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
34ee56fea5 Jean*0679         ENDDO
                0680        ENDDO
                0681 c     ENDIF
                0682 #endif
                0683
e384178ec8 Jean*0684 C---- Experiment specific code: starts here --------------------------
34ee56fea5 Jean*0685       kDiffCyl = 3. _d -7
19780451ef Andr*0686       rSource  = 3
                0687       iGSource = 30
                0688       jGSource = 8
34ee56fea5 Jean*0689       DO j=1,sNy
                0690        DO i=1,sNx
19780451ef Andr*0691         dFluxW = 0.
                0692         dFluxE = 0.
                0693         dFluxN = 0.
                0694         dFluxS = 0.
                0695         jG = myYGlobalLo-1+(bj-1)*sNy+J
                0696         iG = myXGlobalLo-1+(bi-1)*sNx+I
                0697 c      IF(jG .GE. jGSource-rSource .AND. jG .LE. jGSource+rSource) THEN
963e35159b Andr*0698 c     the following bites the big one
19780451ef Andr*0699       IF(jG .LE. 10) THEN
963e35159b Andr*0700          tCyl = tCylIn
2e53480479 Jean*0701       ELSE
963e35159b Andr*0702          tCyl = tCylOut
19780451ef Andr*0703       ENDIF
                0704 c      IF(iG .GE. iGSource-rSource .AND. iG .LE. iGSource+rSource) THEN
df4a605949 Jean*0705           hC = hFacC(i  ,j  ,k,bi,bj)
                0706           hW = hFacW(i  ,j  ,k,bi,bj)
                0707           hE = hFacW(i+1,j  ,k,bi,bj)
                0708           hN = hFacS(i  ,j+1,k,bi,bj)
                0709           hS = hFacS(i  ,j  ,k,bi,bj)
19780451ef Andr*0710           IF ( hC .NE. 0. .AND .hW .EQ. 0. ) THEN
                0711 C          Cylinder to west
df4a605949 Jean*0712            faceArea = drF(k)*dyG(i,j,bi,bj)
2e53480479 Jean*0713            dFluxW =
df4a605949 Jean*0714      &      -faceArea*kDiffCyl*(theta(i,j,k,bi,bj) - tCyl)
19780451ef Andr*0715      &      *recip_dxC(i,j,bi,bj)
                0716           ENDIF
                0717           IF ( hC .NE. 0. .AND .hE .EQ. 0. ) THEN
                0718 C          Cylinder to east
df4a605949 Jean*0719            faceArea = drF(k)*dyG(i+1,j,bi,bj)
19780451ef Andr*0720            dFluxE =
df4a605949 Jean*0721      &      -faceArea*kDiffCyl*(tCyl - theta(i,j,k,bi,bj))
19780451ef Andr*0722      &      *recip_dxC(i,j,bi,bj)
                0723           ENDIF
                0724           IF ( hC .NE. 0. .AND .hN .EQ. 0. ) THEN
                0725 C          Cylinder to north
df4a605949 Jean*0726            faceArea = drF(k)*dxG(i,j+1,bi,bj)
19780451ef Andr*0727            dFluxN =
df4a605949 Jean*0728      &      -faceArea*kDiffCyl*(tCyl-theta(i,j,k,bi,bj))
19780451ef Andr*0729      &      *recip_dyC(i,j,bi,bj)
                0730           ENDIF
                0731           IF ( hC .NE. 0. .AND .hS .EQ. 0. ) THEN
                0732 C          Cylinder to south
df4a605949 Jean*0733            faceArea = drF(k)*dxG(i,j,bi,bj)
19780451ef Andr*0734            dFluxS =
df4a605949 Jean*0735      &      -faceArea*kDiffCyl*(theta(i,j,k,bi,bj) - tCyl)
19780451ef Andr*0736      &      *recip_dyC(i,j,bi,bj)
                0737           ENDIF
                0738 c      ENDIF
                0739 c      ENDIF
                0740 C       Net tendency term is minus flux divergence divided by the volume.
df4a605949 Jean*0741         gT_arr(i,j) = gT_arr(i,j)
                0742      &  -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)
19780451ef Andr*0743      &  *recip_rA(i,j,bi,bj)
                0744      &  *(
                0745      &    dFluxE-dFluxW
                0746      &   +dFluxN-dFluxS
                0747      &   )
                0748
                0749        ENDDO
                0750       ENDDO
e384178ec8 Jean*0751 C---- Experiment specific code:  ends  here --------------------------
                0752
                0753 #ifdef ALLOW_FRAZIL
                0754       IF ( useFRAZIL )
                0755      &     CALL FRAZIL_TENDENCY_APPLY_T(
df4a605949 Jean*0756      U                 gT_arr,
                0757      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0758      I                 myTime, 0, myThid )
                0759 #endif /* ALLOW_FRAZIL */
                0760
                0761 #ifdef ALLOW_SHELFICE
                0762       IF ( useShelfIce )
                0763      &     CALL SHELFICE_FORCING_T(
df4a605949 Jean*0764      U                   gT_arr,
                0765      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0766      I                   myTime, 0, myThid )
                0767 #endif /* ALLOW_SHELFICE */
                0768
                0769 #ifdef ALLOW_ICEFRONT
                0770       IF ( useICEFRONT )
                0771      &     CALL ICEFRONT_TENDENCY_APPLY_T(
df4a605949 Jean*0772      U                   gT_arr,
                0773      I                   k, bi, bj, myTime, 0, myThid )
e384178ec8 Jean*0774 #endif /* ALLOW_ICEFRONT */
                0775
                0776 #ifdef ALLOW_SALT_PLUME
                0777       IF ( useSALT_PLUME )
                0778      &     CALL SALT_PLUME_TENDENCY_APPLY_T(
df4a605949 Jean*0779      U                     gT_arr,
                0780      I                     iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0781      I                     myTime, 0, myThid )
                0782 #endif /* ALLOW_SALT_PLUME */
                0783
                0784 #ifdef ALLOW_RBCS
                0785       IF (useRBCS) THEN
                0786         CALL RBCS_ADD_TENDENCY(
df4a605949 Jean*0787      U                 gT_arr,
                0788      I                 k, bi, bj, 1,
e384178ec8 Jean*0789      I                 myTime, 0, myThid )
                0790       ENDIF
                0791 #endif /* ALLOW_RBCS */
19780451ef Andr*0792
34ee56fea5 Jean*0793 #ifdef ALLOW_OBCS
                0794       IF (useOBCS) THEN
e384178ec8 Jean*0795         CALL OBCS_SPONGE_T(
df4a605949 Jean*0796      U                   gT_arr,
                0797      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0798      I                   myTime, 0, myThid )
34ee56fea5 Jean*0799       ENDIF
e384178ec8 Jean*0800 #endif /* ALLOW_OBCS */
                0801
                0802 #ifdef ALLOW_BBL
                0803       IF ( useBBL ) CALL BBL_TENDENCY_APPLY_T(
df4a605949 Jean*0804      U                       gT_arr,
                0805      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0806      I                       myTime, 0, myThid )
                0807 #endif /* ALLOW_BBL */
                0808
                0809 #ifdef ALLOW_MYPACKAGE
                0810       IF ( useMYPACKAGE ) THEN
                0811         CALL MYPACKAGE_TENDENCY_APPLY_T(
df4a605949 Jean*0812      U                 gT_arr,
                0813      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0814      I                 myTime, 0, myThid )
                0815       ENDIF
                0816 #endif /* ALLOW_MYPACKAGE */
34ee56fea5 Jean*0817
df4a605949 Jean*0818 #endif /* USE_OLD_EXTERNAL_FORCING */
                0819
19780451ef Andr*0820       RETURN
                0821       END
34ee56fea5 Jean*0822
                0823 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
19780451ef Andr*0824 CBOP
df4a605949 Jean*0825 C     !ROUTINE: APPLY_FORCING_S
19780451ef Andr*0826 C     !INTERFACE:
df4a605949 Jean*0827       SUBROUTINE APPLY_FORCING_S(
                0828      U                     gS_arr,
                0829      I                     iMin,iMax,jMin,jMax, k, bi, bj,
                0830      I                     myTime, myIter, myThid )
19780451ef Andr*0831 C     !DESCRIPTION: \bv
                0832 C     *==========================================================*
df4a605949 Jean*0833 C     | S/R APPLY_FORCING_S
34ee56fea5 Jean*0834 C     | o Contains problem specific forcing for merid velocity.
19780451ef Andr*0835 C     *==========================================================*
34ee56fea5 Jean*0836 C     | Adds terms to gS for forcing by external sources
                0837 C     | e.g. fresh-water flux, climatalogical relaxation, etc ...
19780451ef Andr*0838 C     *==========================================================*
                0839 C     \ev
                0840
                0841 C     !USES:
                0842       IMPLICIT NONE
                0843 C     == Global data ==
                0844 #include "SIZE.h"
                0845 #include "EEPARAMS.h"
                0846 #include "PARAMS.h"
                0847 #include "GRID.h"
                0848 #include "DYNVARS.h"
                0849 #include "FFIELDS.h"
e384178ec8 Jean*0850 #include "SURFACE.h"
19780451ef Andr*0851
                0852 C     !INPUT/OUTPUT PARAMETERS:
df4a605949 Jean*0853 C     gS_arr    :: the tendency array
34ee56fea5 Jean*0854 C     iMin,iMax :: Working range of x-index for applying forcing.
                0855 C     jMin,jMax :: Working range of y-index for applying forcing.
df4a605949 Jean*0856 C     k         :: Current vertical level index
34ee56fea5 Jean*0857 C     bi,bj     :: Current tile indices
                0858 C     myTime    :: Current time in simulation
df4a605949 Jean*0859 C     myIter    :: Current iteration number
                0860 C     myThid    :: my Thread Id number
                0861       _RL     gS_arr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0862       INTEGER iMin, iMax, jMin, jMax
                0863       INTEGER k, bi, bj
                0864       _RL     myTime
                0865       INTEGER myIter
19780451ef Andr*0866       INTEGER myThid
                0867
                0868 C     !LOCAL VARIABLES:
34ee56fea5 Jean*0869 C     i,j       :: Loop counters
52d51b46b8 Jean*0870 C     kSurface  :: index of surface level
34ee56fea5 Jean*0871       INTEGER i, j
df4a605949 Jean*0872 #ifdef USE_OLD_EXTERNAL_FORCING
                0873       _RL     locVar(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0874       _RL     tmpVar
                0875 #else
34ee56fea5 Jean*0876       INTEGER kSurface
df4a605949 Jean*0877 #endif
19780451ef Andr*0878 CEOP
                0879
df4a605949 Jean*0880 #ifdef USE_OLD_EXTERNAL_FORCING
                0881
                0882       DO j=1-OLy,sNy+OLy
                0883         DO i=1-OLx,sNx+OLx
                0884           locVar(i,j) = gS(i,j,k,bi,bj)
                0885         ENDDO
                0886       ENDDO
                0887       CALL EXTERNAL_FORCING_S(
                0888      I              iMin, iMax, jMin, jMax, bi, bj, k,
                0889      I              myTime, myThid )
                0890       DO j=1-OLy,sNy+OLy
                0891         DO i=1-OLx,sNx+OLx
                0892           tmpVar = gS(i,j,k,bi,bj) - locVar(i,j)
                0893           gS(i,j,k,bi,bj) = locVar(i,j)
                0894           gS_arr(i,j) = gS_arr(i,j) + tmpVar
                0895         ENDDO
                0896       ENDDO
                0897
                0898 #else  /* USE_OLD_EXTERNAL_FORCING */
                0899
34ee56fea5 Jean*0900       IF ( fluidIsAir ) THEN
                0901        kSurface = 0
e384178ec8 Jean*0902       ELSEIF ( usingZCoords .AND. useShelfIce ) THEN
                0903        kSurface = -1
34ee56fea5 Jean*0904       ELSEIF ( usingPCoords ) THEN
                0905        kSurface = Nr
                0906       ELSE
                0907        kSurface = 1
                0908       ENDIF
                0909
19780451ef Andr*0910 C--   Forcing term
34ee56fea5 Jean*0911 #ifdef ALLOW_AIM
                0912       IF ( useAIM ) CALL AIM_TENDENCY_APPLY_S(
df4a605949 Jean*0913      U                       gS_arr,
                0914      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0915      I                       myTime, 0, myThid )
34ee56fea5 Jean*0916 #endif /* ALLOW_AIM */
                0917
e384178ec8 Jean*0918 #ifdef ALLOW_ATM_PHYS
                0919       IF ( useAtm_Phys ) CALL ATM_PHYS_TENDENCY_APPLY_S(
df4a605949 Jean*0920      U                       gS_arr,
                0921      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0922      I                       myTime, 0, myThid )
                0923 #endif /* ALLOW_ATM_PHYS */
                0924
34ee56fea5 Jean*0925 #ifdef ALLOW_FIZHI
                0926       IF ( useFIZHI ) CALL FIZHI_TENDENCY_APPLY_S(
df4a605949 Jean*0927      U                       gS_arr,
                0928      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0929      I                       myTime, 0, myThid )
34ee56fea5 Jean*0930 #endif /* ALLOW_FIZHI */
                0931
e384178ec8 Jean*0932 #ifdef ALLOW_ADDFLUID
                0933       IF ( selectAddFluid.NE.0 .AND. salt_addMass.NE.UNSET_RL ) THEN
                0934        IF ( ( selectAddFluid.GE.1 .AND. nonlinFreeSurf.GT.0 )
                0935      &      .OR. convertFW2Salt.EQ.-1. _d 0 ) THEN
                0936          DO j=1,sNy
                0937           DO i=1,sNx
df4a605949 Jean*0938             gS_arr(i,j) = gS_arr(i,j)
                0939      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0940      &          *( salt_addMass - salt(i,j,k,bi,bj) )
e384178ec8 Jean*0941      &          *recip_rA(i,j,bi,bj)
df4a605949 Jean*0942      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0943 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
e384178ec8 Jean*0944           ENDDO
                0945          ENDDO
                0946        ELSE
                0947          DO j=1,sNy
                0948           DO i=1,sNx
df4a605949 Jean*0949             gS_arr(i,j) = gS_arr(i,j)
                0950      &        + addMass(i,j,k,bi,bj)*mass2rUnit
                0951      &          *( salt_addMass - sRef(k) )
e384178ec8 Jean*0952      &          *recip_rA(i,j,bi,bj)
df4a605949 Jean*0953      &          *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
                0954 C    &          *recip_deepFac2C(k)*recip_rhoFacC(k)
e384178ec8 Jean*0955           ENDDO
                0956          ENDDO
                0957        ENDIF
                0958       ENDIF
                0959 #endif /* ALLOW_ADDFLUID */
                0960
                0961 C     Ocean: Add salinity surface forcing (e.g., fresh-water) in surface level
df4a605949 Jean*0962       IF ( k .EQ. kSurface ) THEN
34ee56fea5 Jean*0963        DO j=1,sNy
                0964         DO i=1,sNx
df4a605949 Jean*0965           gS_arr(i,j) = gS_arr(i,j)
e384178ec8 Jean*0966      &      +surfaceForcingS(i,j,bi,bj)
df4a605949 Jean*0967      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0968         ENDDO
                0969        ENDDO
                0970       ELSEIF ( kSurface.EQ.-1 ) THEN
                0971        DO j=1,sNy
                0972         DO i=1,sNx
df4a605949 Jean*0973          IF ( kSurfC(i,j,bi,bj).EQ.k ) THEN
                0974           gS_arr(i,j) = gS_arr(i,j)
e384178ec8 Jean*0975      &      +surfaceForcingS(i,j,bi,bj)
df4a605949 Jean*0976      &      *recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0977          ENDIF
19780451ef Andr*0978         ENDDO
                0979        ENDDO
                0980       ENDIF
                0981
e384178ec8 Jean*0982       IF (linFSConserveTr) THEN
                0983        DO j=1,sNy
                0984         DO i=1,sNx
df4a605949 Jean*0985           IF (k .EQ. kSurfC(i,j,bi,bj)) THEN
                0986             gS_arr(i,j) = gS_arr(i,j)
                0987      &        +SsurfCor*recip_drF(k)*_recip_hFacC(i,j,k,bi,bj)
e384178ec8 Jean*0988           ENDIF
                0989         ENDDO
                0990        ENDDO
                0991       ENDIF
                0992
                0993 #ifdef ALLOW_SHELFICE
                0994       IF ( useShelfIce )
                0995      &     CALL SHELFICE_FORCING_S(
df4a605949 Jean*0996      U                   gS_arr,
                0997      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*0998      I                   myTime, 0, myThid )
                0999 #endif /* ALLOW_SHELFICE */
                1000
                1001 #ifdef ALLOW_ICEFRONT
                1002       IF ( useICEFRONT )
                1003      &     CALL ICEFRONT_TENDENCY_APPLY_S(
df4a605949 Jean*1004      U                   gS_arr,
                1005      I                   k, bi, bj, myTime, 0, myThid )
e384178ec8 Jean*1006 #endif /* ALLOW_ICEFRONT */
                1007
                1008 #ifdef ALLOW_SALT_PLUME
                1009       IF ( useSALT_PLUME )
                1010      &     CALL SALT_PLUME_TENDENCY_APPLY_S(
df4a605949 Jean*1011      U                     gS_arr,
                1012      I                     iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*1013      I                     myTime, 0, myThid )
                1014 #endif /* ALLOW_SALT_PLUME */
                1015
                1016 #ifdef ALLOW_RBCS
                1017       IF (useRBCS) THEN
                1018         CALL RBCS_ADD_TENDENCY(
df4a605949 Jean*1019      U                 gS_arr,
                1020      I                 k, bi, bj, 2,
e384178ec8 Jean*1021      I                 myTime, 0, myThid )
                1022       ENDIF
                1023 #endif /* ALLOW_RBCS */
                1024
34ee56fea5 Jean*1025 #ifdef ALLOW_OBCS
19780451ef Andr*1026       IF (useOBCS) THEN
e384178ec8 Jean*1027         CALL OBCS_SPONGE_S(
df4a605949 Jean*1028      U                   gS_arr,
                1029      I                   iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*1030      I                   myTime, 0, myThid )
19780451ef Andr*1031       ENDIF
e384178ec8 Jean*1032 #endif /* ALLOW_OBCS */
                1033
                1034 #ifdef ALLOW_BBL
                1035       IF ( useBBL ) CALL BBL_TENDENCY_APPLY_S(
df4a605949 Jean*1036      U                       gS_arr,
                1037      I                       iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*1038      I                       myTime, 0, myThid )
                1039 #endif /* ALLOW_BBL */
                1040
                1041 #ifdef ALLOW_MYPACKAGE
                1042       IF ( useMYPACKAGE ) THEN
                1043         CALL MYPACKAGE_TENDENCY_APPLY_S(
df4a605949 Jean*1044      U                 gS_arr,
                1045      I                 iMin,iMax,jMin,jMax, k, bi,bj,
e384178ec8 Jean*1046      I                 myTime, 0, myThid )
                1047       ENDIF
                1048 #endif /* ALLOW_MYPACKAGE */
19780451ef Andr*1049
df4a605949 Jean*1050 #endif /* USE_OLD_EXTERNAL_FORCING */
                1051
19780451ef Andr*1052       RETURN
                1053       END