pkg/thsice/thsice_step_fwd.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.33 2012/02/04 17:21:23 jmc Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"
#ifdef ALLOW_ATM2D
# include "ctrparam.h"
#endif
#define THSICE_OLD_STEP_FWD

CBOP
C     !ROUTINE: THSICE_STEP_FWD
C     !INTERFACE:
      SUBROUTINE THSICE_STEP_FWD(
     I             bi, bj, iMin, iMax, jMin, jMax,
     I             prcAtm,
     I             myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R  THSICE_STEP_FWD
C     | o Step Forward Therm-SeaIce model.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#ifdef  ALLOW_ATM2D
# include "ATMSIZE.h"
# include "ATM2D_VARS.h"
#endif
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#include "THSICE_TAVE.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

      INTEGER siLo, siHi, sjLo, sjHi
      PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
      PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C- input:
C     bi,bj   :: tile indices
C   iMin,iMax :: computation domain: 1rst index range
C   jMin,jMax :: computation domain: 2nd  index range
C     prcAtm  :: total precip from the atmosphere [kg/m2/s]
C     myTime  :: current Time of simulation [s]
C     myIter  :: current Iteration number in simulation
C     myThid  :: my Thread Id number
C-- Use fluxes hold in commom blocks
C- input:
C     icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
C     icFlxAtm  :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
C     icFrwAtm  :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
C- output
C     icFlxAtm  :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
C     icFrwAtm  :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL  myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_THSICE
C     !LOCAL VARIABLES:
C     === Local variables ===
C     iceFrac   :: fraction of grid area covered in ice
C     flx2oc    :: net heat flux from the ice to the ocean (+=down) [W/m2]
C     frw2oc    :: fresh-water flux from the ice to the ocean
C     fsalt     :: mass salt flux to the ocean
C     frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
C     tFrzOce   :: sea-water freezing temperature [oC] (function of S)
C     isIceFree :: true for ice-free grid-cell that remains ice-free
C     ageFac    :: snow aging factor [1]
C     snowFac   :: snowing refreshing-age factor [units of 1/snowPr]
      LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     iceFrac  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     flx2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     frw2oc   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     fsalt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     tFrzOce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL     frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL ageFac
      _RL snowFac
      _RL cphm
      _RL opFrac, icFrac
#ifdef ALLOW_DIAGNOSTICS
      _RL tmpFac
#endif
      INTEGER i,j
      LOGICAL dBugFlag

C-    define grid-point location where to print debugging values
#include "THSICE_DEBUG.h"

 1010 FORMAT(A,1P4E14.6)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_AUTODIFF_TAMC
      act1 = bi - myBxLo(myThid)
      max1 = myBxHi(myThid) - myBxLo(myThid) + 1
      act2 = bj - myByLo(myThid)
      max2 = myByHi(myThid) - myByLo(myThid) + 1
      act3 = myThid - 1
      max3 = nTx*nTy
      act4 = ikey_dynamics - 1
      ticekey = (act1 + 1) + act2*max1
     &                     + act3*max1*max2
     &                     + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

C-    Initialise
      dBugFlag = debugLevel.GE.debLevC
      DO j = 1-OLy, sNy+OLy
       DO i = 1-OLx, sNx+OLx
          isIceFree(i,j) = .FALSE.
#ifdef ALLOW_ATM2D
          sFluxFromIce(i,j) = 0. _d 0
#else
          saltFlux(i,j,bi,bj) = 0. _d 0
#endif
#ifdef ALLOW_AUTODIFF_TAMC
          iceFrac(i,j) = 0.
          flx2oc(i,j) = 0. _d 0
          frw2oc(i,j) = 0. _d 0
          fsalt (i,j) = 0. _d 0
#endif
       ENDDO
      ENDDO

      ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
      snowFac = thSIce_deltaT/(rhos*hNewSnowAge)

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
      DO j = jMin, jMax
       DO i = iMin, iMax
        IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
C--  Snow aging :
          snowAge(i,j,bi,bj) = thSIce_deltaT
     &                       + snowAge(i,j,bi,bj)*ageFac
          IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
     &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
     &          * EXP( - snowFac*snowPrc(i,j,bi,bj) )
C-------
C note: Any flux of mass (here fresh water) that enter or leave the system
C       with a non zero energy HAS TO be counted: add snow precip.
          icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
     &                        - Lfresh*snowPrc(i,j,bi,bj)
C--
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
        tmpFac = 1. _d 0
        CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
        CALL DIAGNOSTICS_FRACT_FILL(
     I                   snowPrc,   iceMask,tmpFac,1,'SIsnwPrc',
     I                   0,1,1,bi,bj,myThid)
        CALL DIAGNOSTICS_FRACT_FILL(
     I                   siceAlb,   iceMask,tmpFac,1,'SIalbedo',
     I                   0,1,1,bi,bj,myThid)
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */
      DO j = jMin, jMax
       DO i = iMin, iMax
          siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
       ENDDO
      ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C    part.2 : ice-covered fraction ;
C     change in ice/snow thickness and ice-fraction
C     note: can only reduce the ice-fraction but not increase it.
C-------
      DO j = jMin, jMax
       DO i = iMin, iMax

        tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
        cphm    = cpwater*rhosw*hOceMxL(i,j,bi,bj)
        frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
     &                 * cphm/ocean_deltaT
        iceFrac(i,j) = iceMask(i,j,bi,bj)
        flx2oc(i,j)  = icFlxSW(i,j,bi,bj)
C-------
#ifdef ALLOW_DBUG_THSICE
        IF ( dBug(i,j,bi,bj) ) THEN
         IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
          WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
          WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf  =',
     &                  iceFrac(i,j), iceHeight(i,j,bi,bj),
     &                  snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
          WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
     &                     tOceMxL(i,j,bi,bj), tFrzOce(i,j),
     &                     frzmltMxL(i,j), Qnet(i,j,bi,bj)
         ENDIF
         IF (iceFrac(i,j).GT.0.)
     &    WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
     &      iceFrac(i,j), icFlxAtm(i,j,bi,bj),
     &      icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
        ENDIF
#endif
       ENDDO
      ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
#endif

      CALL THSICE_CALC_THICKN(
     I          bi, bj,
     I          iMin,iMax, jMin,jMax, dBugFlag,
     I          iceMask(siLo,sjLo,bi,bj), tFrzOce,
     I          tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
     I          snowPrc(siLo,sjLo,bi,bj), prcAtm,
     I          sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
     U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
     U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
     U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
     U          icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
     O          frw2oc, fsalt,
     I          myTime, myIter, myThid )

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE fsalt(:,:)  = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE flx2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
CADJ STORE frw2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
#endif
C--    Net fluxes :
      DO j = jMin, jMax
       DO i = iMin, iMax
c#ifdef ALLOW_AUTODIFF_TAMC
c         ikey_1 = i
c    &         + sNx*(j-1)
c    &         + sNx*sNy*act1
c    &         + sNx*sNy*max1*act2
c    &         + sNx*sNy*max1*max2*act3
c    &         + sNx*sNy*max1*max2*max3*act4
c#endif /* ALLOW_AUTODIFF_TAMC */
c#ifdef ALLOW_AUTODIFF_TAMC
cCADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
c#endif
        IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
C-     weighted average net fluxes:
c#ifdef ALLOW_AUTODIFF_TAMC
cCADJ STORE  fsalt(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE  flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE  frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
cCADJ STORE  icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
c#endif
          icFrac = iceMask(i,j,bi,bj)
          opFrac= 1. _d 0-icFrac
#ifdef ALLOW_ATM2D
          pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
          pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
          sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
#else
          icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
     &                        - opFrac*Qnet(i,j,bi,bj)
          icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
     &                        + opFrac*EmPmR(i,j,bi,bj)
          Qnet(i,j,bi,bj) = -icFrac*flx2oc(i,j) + opFrac*Qnet(i,j,bi,bj)
          EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
     &                    +  opFrac*EmPmR(i,j,bi,bj)
          saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
#endif

#ifdef ALLOW_DBUG_THSICE
          IF (dBug(i,j,bi,bj)) WRITE(6,1010)
     &          'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
     &           iceFrac(i,j), iceHeight(i,j,bi,bj),
     &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
#endif

        ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
          icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
          icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
        ELSE
          icFlxAtm(i,j,bi,bj) = 0. _d 0
          icFrwAtm(i,j,bi,bj) = 0. _d 0
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_SALT_PLUME
      IF ( useSALT_PLUME ) THEN
        CALL THSICE_SALT_PLUME(
     I          iceMask(1-OLx,1-OLy,bi,bj),
     I          sOceMxL(1-OLx,1-OLy,bi,bj),
     I          frw2oc,
     I          iMin,iMax, jMin,jMax, bi, bj,
     I          0, myTime, myIter, myThid )
      ENDIF
#endif /* ALLOW_SALT_PLUME */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
C    part.3 : freezing of sea-water
C     over ice-free fraction and what is left from ice-covered fraction
C-------
      DO j = 1-OLy, sNy+OLy
       DO i = 1-OLx, sNx+OLx
         flx2oc(i,j) = 0. _d 0
         frw2oc(i,j) = 0. _d 0
         fsalt (i,j) = 0. _d 0
       ENDDO
      ENDDO
      CALL THSICE_EXTEND(
     I          bi, bj,
     I          iMin,iMax, jMin,jMax, dBugFlag,
     I          frzmltMxL, tFrzOce,
     I          tOceMxL(siLo,sjLo,bi,bj),
     U          iceFrac, iceHeight(siLo,sjLo,bi,bj),
     U          snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
     U          Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
     U          Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
     O          flx2oc, frw2oc, fsalt,
     I          myTime, myIter, myThid )

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE snowHeight(:,:,bi,bj) =
CADJ &     comlev1_bibj,key=ticekey,byte=isbyte
#endif
      DO j = jMin, jMax
       DO i = iMin, iMax
#ifdef THSICE_OLD_STEP_FWD
        IF (frzmltMxL(i,j).GT.0. _d 0) THEN
#endif
C--    Net fluxes :
#ifdef ALLOW_ATM2D
          pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
          pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
          sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
#else
          Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
          EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)
          saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
#endif

#ifdef ALLOW_DBUG_THSICE
          IF (dBug(i,j,bi,bj)) WRITE(6,1010)
     &         'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
     &           iceFrac(i,j), iceHeight(i,j,bi,bj),
     &           snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
#endif
#ifdef THSICE_OLD_STEP_FWD
        ENDIF
#endif

        IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
     &    isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
     &                  .AND.   iceFrac(i,j) .LE.0. _d 0
        IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
          iceMask(i,j,bi,bj)=iceFrac(i,j)
          IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
     &     snowAge(i,j,bi,bj) = 0. _d 0
        ELSE
          iceMask(i,j,bi,bj)  = 0. _d 0
          iceHeight(i,j,bi,bj)= 0. _d 0
          snowHeight(i,j,bi,bj)=0. _d 0
          snowAge(i,j,bi,bj)  = 0. _d 0
          Tsrf(i,j,bi,bj)     = tOceMxL(i,j,bi,bj)
          Tice1(i,j,bi,bj)    = 0. _d 0
          Tice2(i,j,bi,bj)    = 0. _d 0
          Qice1(i,j,bi,bj)    = Lfresh
          Qice2(i,j,bi,bj)    = Lfresh
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_SALT_PLUME
      IF ( useSALT_PLUME ) THEN
        CALL THSICE_SALT_PLUME(
     I          frzmltMxL,
     I          sOceMxL(1-OLx,1-OLy,bi,bj),
     I          frw2oc,
     I          iMin,iMax, jMin,jMax, bi, bj,
     I          1, myTime, myIter, myThid )
      ENDIF
#endif /* ALLOW_SALT_PLUME */

# ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE snowHeight(:,:,bi,bj) =
CADJ &     comlev1_bibj,key=ticekey,byte=isbyte
# endif
      DO j = jMin, jMax
       DO i = iMin, iMax
C--     Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
        sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
     &                         + iceHeight(i,j,bi,bj)*rhoi
     &                        )*iceMask(i,j,bi,bj)
#ifdef ALLOW_ATM2D
        pass_sIceLoad(i,j)=sIceLoad(i,j,bi,bj)
#endif
       ENDDO
      ENDDO

      IF ( thSIceAdvScheme.GT.0 ) THEN
C--   note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
       DO j = jMin, jMax
        DO i = iMin, iMax
         IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
          Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
          EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
          saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
         ENDIF
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_BULK_FORCE
      IF ( useBulkForce ) THEN
        CALL BULKF_FLUX_ADJUST(
     I                          bi, bj, iMin, iMax, jMin, jMax,
     I                          isIceFree, myTime, myIter, myThid )
      ENDIF
#endif /* ALLOW_BULK_FORCE */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#endif /* ALLOW_THSICE */

      RETURN
      END
1	jmc	1.34	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.33 2012/02/04 17:21:23 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	jscott	1.23	#ifdef ALLOW_ATM2D
6			# include "ctrparam.h"
7			#endif
8	jmc	1.34	#define THSICE_OLD_STEP_FWD
9	jmc	1.3
10			CBOP
11	jmc	1.1	C !ROUTINE: THSICE_STEP_FWD
12			C !INTERFACE:
13	jmc	1.13	SUBROUTINE THSICE_STEP_FWD(
14	jmc	1.1	I bi, bj, iMin, iMax, jMin, jMax,
15	jmc	1.3	I prcAtm,
16	jmc	1.1	I myTime, myIter, myThid )
17	jmc	1.3	C !DESCRIPTION: \bv
18	jmc	1.1	C ==========================================================
19	jmc	1.13	C \| S/R THSICE_STEP_FWD
20	jmc	1.1	C \| o Step Forward Therm-SeaIce model.
21			C ==========================================================
22	jmc	1.3	C \ev
23	jmc	1.1
24			C !USES:
25			IMPLICIT NONE
26	jmc	1.3
27	jmc	1.1	C === Global variables ===
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "FFIELDS.h"
32	jscott	1.23	#ifdef ALLOW_ATM2D
33			# include "ATMSIZE.h"
34			# include "ATM2D_VARS.h"
35			#endif
36	jmc	1.1	#include "THSICE_SIZE.h"
37			#include "THSICE_PARAMS.h"
38	jmc	1.3	#include "THSICE_VARS.h"
39			#include "THSICE_TAVE.h"
40	heimbach	1.20	#ifdef ALLOW_AUTODIFF_TAMC
41			# include "tamc.h"
42			# include "tamc_keys.h"
43			#endif
44
45	jmc	1.17	INTEGER siLo, siHi, sjLo, sjHi
46			PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
47			PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
48	jmc	1.13
49	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
50			C === Routine arguments ===
51	jmc	1.17	C- input:
52	jmc	1.3	C bi,bj :: tile indices
53			C iMin,iMax :: computation domain: 1rst index range
54			C jMin,jMax :: computation domain: 2nd index range
55	jmc	1.17	C prcAtm :: total precip from the atmosphere [kg/m2/s]
56			C myTime :: current Time of simulation [s]
57			C myIter :: current Iteration number in simulation
58			C myThid :: my Thread Id number
59			C-- Use fluxes hold in commom blocks
60	jmc	1.3	C- input:
61	jmc	1.17	C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
62			C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
63			C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
64	jmc	1.3	C- output
65	jmc	1.17	C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
66			C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
67	jmc	1.1	INTEGER bi,bj
68			INTEGER iMin, iMax
69			INTEGER jMin, jMax
70	jmc	1.3	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71	jmc	1.1	_RL myTime
72			INTEGER myIter
73			INTEGER myThid
74	jmc	1.3	CEOP
75	jmc	1.1
76			#ifdef ALLOW_THSICE
77			C !LOCAL VARIABLES:
78			C === Local variables ===
79	jmc	1.17	C iceFrac :: fraction of grid area covered in ice
80	jmc	1.3	C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
81			C frw2oc :: fresh-water flux from the ice to the ocean
82			C fsalt :: mass salt flux to the ocean
83			C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
84	jmc	1.17	C tFrzOce :: sea-water freezing temperature [oC] (function of S)
85	jmc	1.14	C isIceFree :: true for ice-free grid-cell that remains ice-free
86	jmc	1.17	C ageFac :: snow aging factor [1]
87			C snowFac :: snowing refreshing-age factor [units of 1/snowPr]
88			LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89			_RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90			_RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91			_RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92			_RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93			_RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
94			_RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95			_RL ageFac
96			_RL snowFac
97			_RL cphm
98	jmc	1.3	_RL opFrac, icFrac
99	jmc	1.9	#ifdef ALLOW_DIAGNOSTICS
100	jmc	1.10	_RL tmpFac
101	jmc	1.9	#endif
102	jmc	1.17	INTEGER i,j
103			LOGICAL dBugFlag
104	jmc	1.1
105	jmc	1.17	C- define grid-point location where to print debugging values
106			#include "THSICE_DEBUG.h"
107	jmc	1.1
108	jmc	1.13	1010 FORMAT(A,1P4E14.6)
109	jmc	1.17
110			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
111
112	heimbach	1.20	#ifdef ALLOW_AUTODIFF_TAMC
113			act1 = bi - myBxLo(myThid)
114			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
115			act2 = bj - myByLo(myThid)
116			max2 = myByHi(myThid) - myByLo(myThid) + 1
117			act3 = myThid - 1
118			max3 = nTx*nTy
119			act4 = ikey_dynamics - 1
120	gforget	1.30	ticekey = (act1 + 1) + act2*max1
121	heimbach	1.20	& + act3max1max2
122			& + act4max1max2*max3
123			#endif /* ALLOW_AUTODIFF_TAMC */
124
125	jmc	1.17	C- Initialise
126	jmc	1.31	dBugFlag = debugLevel.GE.debLevC
127	jmc	1.17	DO j = 1-OLy, sNy+OLy
128	jmc	1.34	DO i = 1-OLx, sNx+OLx
129	jmc	1.14	isIceFree(i,j) = .FALSE.
130	jscott	1.23	#ifdef ALLOW_ATM2D
131			sFluxFromIce(i,j) = 0. _d 0
132			#else
133	jmc	1.17	saltFlux(i,j,bi,bj) = 0. _d 0
134	jscott	1.23	#endif
135	heimbach	1.16	#ifdef ALLOW_AUTODIFF_TAMC
136	jmc	1.17	iceFrac(i,j) = 0.
137	jmc	1.34	flx2oc(i,j) = 0. _d 0
138			frw2oc(i,j) = 0. _d 0
139			fsalt (i,j) = 0. _d 0
140	heimbach	1.16	#endif
141	jmc	1.34	ENDDO
142	jmc	1.17	ENDDO
143	jmc	1.1
144	jmc	1.17	ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
145			snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
146	heimbach	1.21
147			#ifdef ALLOW_AUTODIFF_TAMC
148	gforget	1.30	CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
149			CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
150			CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
151			CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
152			CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
153			CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
154	heimbach	1.21	#endif
155	jmc	1.17	DO j = jMin, jMax
156			DO i = iMin, iMax
157			IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
158			C-- Snow aging :
159			snowAge(i,j,bi,bj) = thSIce_deltaT
160			& + snowAge(i,j,bi,bj)*ageFac
161			IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
162			& snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
163			& * EXP( - snowFac*snowPrc(i,j,bi,bj) )
164	jmc	1.3	C-------
165	jmc	1.17	C note: Any flux of mass (here fresh water) that enter or leave the system
166			C with a non zero energy HAS TO be counted: add snow precip.
167			icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
168			& - Lfresh*snowPrc(i,j,bi,bj)
169			C--
170			ENDIF
171	jmc	1.3	ENDDO
172	jmc	1.17	ENDDO
173	jmc	1.3
174	jmc	1.9	#ifdef ALLOW_DIAGNOSTICS
175			IF ( useDiagnostics ) THEN
176	jmc	1.10	tmpFac = 1. _d 0
177	jmc	1.22	CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
178	jmc	1.10	CALL DIAGNOSTICS_FRACT_FILL(
179	jmc	1.11	I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
180	jmc	1.10	I 0,1,1,bi,bj,myThid)
181	jmc	1.12	CALL DIAGNOSTICS_FRACT_FILL(
182			I siceAlb, iceMask,tmpFac,1,'SIalbedo',
183			I 0,1,1,bi,bj,myThid)
184	jmc	1.9	ENDIF
185			#endif /* ALLOW_DIAGNOSTICS */
186	jmc	1.12	DO j = jMin, jMax
187			DO i = iMin, iMax
188			siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
189			ENDDO
190			ENDDO
191	jmc	1.9
192	jmc	1.3	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
193	jmc	1.13	C part.2 : ice-covered fraction ;
194			C change in ice/snow thickness and ice-fraction
195	jmc	1.3	C note: can only reduce the ice-fraction but not increase it.
196			C-------
197	jmc	1.1	DO j = jMin, jMax
198			DO i = iMin, iMax
199	jmc	1.3
200	jmc	1.17	tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
201	jmc	1.3	cphm = cpwaterrhoswhOceMxL(i,j,bi,bj)
202	jmc	1.17	frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
203			& * cphm/ocean_deltaT
204			iceFrac(i,j) = iceMask(i,j,bi,bj)
205			flx2oc(i,j) = icFlxSW(i,j,bi,bj)
206	jmc	1.3	C-------
207	jmc	1.17	#ifdef ALLOW_DBUG_THSICE
208			IF ( dBug(i,j,bi,bj) ) THEN
209			IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
210	jmc	1.5	WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
211			WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
212	jmc	1.17	& iceFrac(i,j), iceHeight(i,j,bi,bj),
213	jmc	1.7	& snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
214	jmc	1.17	WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
215			& tOceMxL(i,j,bi,bj), tFrzOce(i,j),
216			& frzmltMxL(i,j), Qnet(i,j,bi,bj)
217			ENDIF
218			IF (iceFrac(i,j).GT.0.)
219			& WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
220			& iceFrac(i,j), icFlxAtm(i,j,bi,bj),
221			& icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
222	jmc	1.1	ENDIF
223	jmc	1.17	#endif
224			ENDDO
225			ENDDO
226	jmc	1.1
227	heimbach	1.20	#ifdef ALLOW_AUTODIFF_TAMC
228	gforget	1.30	CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
229	heimbach	1.20	#endif
230
231	jmc	1.17	CALL THSICE_CALC_THICKN(
232	heimbach	1.28	I bi, bj,
233	jmc	1.17	I iMin,iMax, jMin,jMax, dBugFlag,
234			I iceMask(siLo,sjLo,bi,bj), tFrzOce,
235			I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
236			I snowPrc(siLo,sjLo,bi,bj), prcAtm,
237			I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
238			U iceFrac, iceHeight(siLo,sjLo,bi,bj),
239			U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
240			U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
241			U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
242			O frw2oc, fsalt,
243			I myTime, myIter, myThid )
244	jmc	1.1
245	jmc	1.29	#ifdef ALLOW_AUTODIFF_TAMC
246	gforget	1.30	CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj,key=ticekey,byte=isbyte
247			CADJ STORE fsalt(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
248			CADJ STORE flx2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
249			CADJ STORE frw2oc(:,:) = comlev1_bibj,key=ticekey,byte=isbyte
250	jmc	1.29	#endif
251	jmc	1.1	C-- Net fluxes :
252	jmc	1.17	DO j = jMin, jMax
253			DO i = iMin, iMax
254	jmc	1.29	c#ifdef ALLOW_AUTODIFF_TAMC
255			c ikey_1 = i
256			c & + sNx*(j-1)
257			c & + sNxsNyact1
258			c & + sNxsNymax1*act2
259			c & + sNxsNymax1max2act3
260			c & + sNxsNymax1max2max3*act4
261			c#endif /* ALLOW_AUTODIFF_TAMC */
262			c#ifdef ALLOW_AUTODIFF_TAMC
263			cCADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
264			c#endif
265	jmc	1.17	IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
266	jmc	1.3	C- weighted average net fluxes:
267	jmc	1.29	c#ifdef ALLOW_AUTODIFF_TAMC
268			cCADJ STORE fsalt(i,j) = comlev1_thsice_1, key=ikey_1
269			cCADJ STORE flx2oc(i,j) = comlev1_thsice_1, key=ikey_1
270			cCADJ STORE frw2oc(i,j) = comlev1_thsice_1, key=ikey_1
271			cCADJ STORE icemask(i,j,bi,bj) = comlev1_thsice_1, key=ikey_1
272			c#endif
273	jmc	1.3	icFrac = iceMask(i,j,bi,bj)
274			opFrac= 1. _d 0-icFrac
275	jscott	1.23	#ifdef ALLOW_ATM2D
276			pass_qnet(i,j) = pass_qnet(i,j) - icFrac*flx2oc(i,j)
277			pass_evap(i,j) = pass_evap(i,j) - icFrac*frw2oc(i,j)/rhofw
278			sFluxFromIce(i,j) = -icFrac*fsalt(i,j)
279			#else
280	jmc	1.17	icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
281			& - opFrac*Qnet(i,j,bi,bj)
282			icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
283	jmc	1.26	& + opFrac*EmPmR(i,j,bi,bj)
284	jmc	1.17	Qnet(i,j,bi,bj) = -icFracflx2oc(i,j) + opFracQnet(i,j,bi,bj)
285	jmc	1.26	EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)
286	jmc	1.17	& + opFrac*EmPmR(i,j,bi,bj)
287			saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
288	jscott	1.23	#endif
289	jmc	1.17
290			#ifdef ALLOW_DBUG_THSICE
291			IF (dBug(i,j,bi,bj)) WRITE(6,1010)
292			& 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
293			& iceFrac(i,j), iceHeight(i,j,bi,bj),
294			& snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
295			#endif
296	jmc	1.1
297	jmc	1.26	ELSEIF (hOceMxL(i,j,bi,bj).GT.0. _d 0) THEN
298	jmc	1.17	icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
299	jmc	1.26	icFrwAtm(i,j,bi,bj) = EmPmR(i,j,bi,bj)
300	jmc	1.3	ELSE
301	jmc	1.17	icFlxAtm(i,j,bi,bj) = 0. _d 0
302			icFrwAtm(i,j,bi,bj) = 0. _d 0
303	jmc	1.1	ENDIF
304	jmc	1.17	ENDDO
305			ENDDO
306	jmc	1.1
307	jmc	1.33	#ifdef ALLOW_SALT_PLUME
308			IF ( useSALT_PLUME ) THEN
309			CALL THSICE_SALT_PLUME(
310	jmc	1.34	I iceMask(1-OLx,1-OLy,bi,bj),
311			I sOceMxL(1-OLx,1-OLy,bi,bj),
312	jmc	1.33	I frw2oc,
313	jmc	1.34	I iMin,iMax, jMin,jMax, bi, bj,
314	jmc	1.33	I 0, myTime, myIter, myThid )
315			ENDIF
316			#endif /* ALLOW_SALT_PLUME */
317
318	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
319	jmc	1.3	C part.3 : freezing of sea-water
320	jmc	1.1	C over ice-free fraction and what is left from ice-covered fraction
321			C-------
322	jmc	1.34	DO j = 1-OLy, sNy+OLy
323			DO i = 1-OLx, sNx+OLx
324			flx2oc(i,j) = 0. _d 0
325			frw2oc(i,j) = 0. _d 0
326			fsalt (i,j) = 0. _d 0
327			ENDDO
328			ENDDO
329	jmc	1.17	CALL THSICE_EXTEND(
330	heimbach	1.28	I bi, bj,
331	jmc	1.17	I iMin,iMax, jMin,jMax, dBugFlag,
332			I frzmltMxL, tFrzOce,
333			I tOceMxL(siLo,sjLo,bi,bj),
334			U iceFrac, iceHeight(siLo,sjLo,bi,bj),
335			U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
336			U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
337			U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
338			O flx2oc, frw2oc, fsalt,
339			I myTime, myIter, myThid )
340
341	heimbach	1.21	#ifdef ALLOW_AUTODIFF_TAMC
342	jmc	1.29	CADJ STORE snowHeight(:,:,bi,bj) =
343	gforget	1.30	CADJ & comlev1_bibj,key=ticekey,byte=isbyte
344	heimbach	1.21	#endif
345	jmc	1.17	DO j = jMin, jMax
346			DO i = iMin, iMax
347	jmc	1.34	#ifdef THSICE_OLD_STEP_FWD
348	jmc	1.17	IF (frzmltMxL(i,j).GT.0. _d 0) THEN
349	jmc	1.34	#endif
350	jmc	1.1	C-- Net fluxes :
351	jscott	1.23	#ifdef ALLOW_ATM2D
352			pass_qnet(i,j) = pass_qnet(i,j) - flx2oc(i,j)
353			pass_evap(i,j) = pass_evap(i,j) - frw2oc(i,j)/rhofw
354			sFluxFromIce(i,j)= sFluxFromIce(i,j) - fsalt(i,j)
355			#else
356	jmc	1.17	Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
357	jmc	1.26	EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)
358	jmc	1.17	saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
359	jscott	1.23	#endif
360	jmc	1.17
361			#ifdef ALLOW_DBUG_THSICE
362			IF (dBug(i,j,bi,bj)) WRITE(6,1010)
363			& 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
364			& iceFrac(i,j), iceHeight(i,j,bi,bj),
365			& snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
366			#endif
367	jmc	1.34	#ifdef THSICE_OLD_STEP_FWD
368	jmc	1.1	ENDIF
369	jmc	1.34	#endif
370	jmc	1.1
371	jmc	1.14	IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
372			& isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
373	jmc	1.18	& .AND. iceFrac(i,j) .LE.0. _d 0
374	jmc	1.17	IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
375			iceMask(i,j,bi,bj)=iceFrac(i,j)
376			IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
377			& snowAge(i,j,bi,bj) = 0. _d 0
378	jmc	1.1	ELSE
379			iceMask(i,j,bi,bj) = 0. _d 0
380			iceHeight(i,j,bi,bj)= 0. _d 0
381			snowHeight(i,j,bi,bj)=0. _d 0
382	jmc	1.3	snowAge(i,j,bi,bj) = 0. _d 0
383	jmc	1.17	Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
384	jmc	1.1	Tice1(i,j,bi,bj) = 0. _d 0
385			Tice2(i,j,bi,bj) = 0. _d 0
386	jmc	1.19	Qice1(i,j,bi,bj) = Lfresh
387			Qice2(i,j,bi,bj) = Lfresh
388	jmc	1.1	ENDIF
389	heimbach	1.21	ENDDO
390			ENDDO
391	jmc	1.1
392	jmc	1.33	#ifdef ALLOW_SALT_PLUME
393			IF ( useSALT_PLUME ) THEN
394			CALL THSICE_SALT_PLUME(
395	jmc	1.34	I frzmltMxL,
396			I sOceMxL(1-OLx,1-OLy,bi,bj),
397	jmc	1.33	I frw2oc,
398	jmc	1.34	I iMin,iMax, jMin,jMax, bi, bj,
399	jmc	1.33	I 1, myTime, myIter, myThid )
400			ENDIF
401			#endif /* ALLOW_SALT_PLUME */
402
403	heimbach	1.21	# ifdef ALLOW_AUTODIFF_TAMC
404	jmc	1.29	CADJ STORE snowHeight(:,:,bi,bj) =
405	gforget	1.30	CADJ & comlev1_bibj,key=ticekey,byte=isbyte
406	heimbach	1.21	# endif
407			DO j = jMin, jMax
408			DO i = iMin, iMax
409	jmc	1.6	C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
410			sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
411			& + iceHeight(i,j,bi,bj)*rhoi
412			& )*iceMask(i,j,bi,bj)
413	jscott	1.32	#ifdef ALLOW_ATM2D
414			pass_sIceLoad(i,j)=sIceLoad(i,j,bi,bj)
415			#endif
416	jmc	1.1	ENDDO
417			ENDDO
418
419	jmc	1.18	IF ( thSIceAdvScheme.GT.0 ) THEN
420			C-- note: those fluxes should to be added directly to Qnet, EmPmR & saltFlux
421			DO j = jMin, jMax
422			DO i = iMin, iMax
423			IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 ) THEN
424			Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - oceQnet(i,j,bi,bj)
425	jmc	1.26	EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- oceFWfx(i,j,bi,bj)
426	jmc	1.18	saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - oceSflx(i,j,bi,bj)
427			ENDIF
428			ENDDO
429			ENDDO
430			ENDIF
431
432	jmc	1.14	#ifdef ALLOW_BULK_FORCE
433			IF ( useBulkForce ) THEN
434			CALL BULKF_FLUX_ADJUST(
435			I bi, bj, iMin, iMax, jMin, jMax,
436			I isIceFree, myTime, myIter, myThid )
437			ENDIF
438			#endif /* ALLOW_BULK_FORCE */
439
440	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
441			#endif /* ALLOW_THSICE */
442
443			RETURN
444			END