Back to home page

MITgcm

 
 

    


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

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