pkg/thsice/thsice_step_fwd.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.42 2013/06/10 20:05:15 jmc Exp $
C $Name:  $

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

CBOP
C     !ROUTINE: THSICE_STEP_FWD
C     !INTERFACE:
      SUBROUTINE THSICE_STEP_FWD(
     I             bi, bj, iMin, iMax, jMin, jMax,
     I             prcAtm, snowPrc, qPrcRnO,
     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     snowPrc :: snow precipitation               [kg/m2/s]
C     qPrcRnO :: Energy content of Precip+RunOff (+=down) [W/m2]
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 [kg/m2/s] (+=up)
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      _RL prcAtm (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL snowPrc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL qPrcRnO(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 (+=down)
C     fsalt     :: mass salt flux to the ocean                (+=down)
C     frzSeaWat :: seawater freezing rate (expressed as mass flux) [kg/m^2/s]
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     frzSeaWat(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
      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
          frzSeaWat(i,j) = 0. _d 0
#ifdef ALLOW_AUTODIFF_TAMC
          iceFrac(i,j) = 0.
C-   set these arrays everywhere: overlap are not set and not used,
C     but some arrays are stored and storage includes overlap.
          flx2oc(i,j) = 0. _d 0
          frw2oc(i,j) = 0. _d 0
          fsalt (i,j) = 0. _d 0
c         tFrzOce  (i,j) = 0. _d 0
c         frzmltMxL(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).GT.0. _d 0 )
     &      snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
     &          * EXP( - snowFac*snowPrc(i,j) )
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)
     &                        + qPrcRnO(i,j)
C--
        ENDIF
       ENDDO
      ENDDO

#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
# ifdef OLD_THSICE_CALL_SEQUENCE
        CALL DIAGNOSTICS_FILL(iceMask,'SI_FrcFx',0,1,1,bi,bj,myThid)
# endif
        CALL DIAGNOSTICS_FRACT_FILL( snowPrc,
     I                       iceMask(1-OLx,1-OLy,bi,bj), oneRL, 1,
     I                              'SIsnwPrc', 0,1,2,bi,bj,myThid )
        CALL DIAGNOSTICS_FRACT_FILL( siceAlb, iceMask, oneRL, 1,
     I                              'SIalbedo', 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) + qPrcRnO(i,j)
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)
        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), 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, frzSeaWat,
     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
C-    All seawater freezing (no reduction by surf. melting) from CALC_THICKN
c         frzSeaWat(i,j) = icFrac*frzSeaWat(i,j)
C-    Net seawater freezing (underestimated if there is surf. melting or rain)
          frzSeaWat(i,j) = MAX( -icFrac*frw2oc(i,j), 0. _d 0 )

#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

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
C--    Net fluxes : (only non-zero contribution where frzmltMxL > 0 )
#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
          frzSeaWat(i,j) = frzSeaWat(i,j) + MAX(-frw2oc(i,j), 0. _d 0 )

#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

        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          sOceMxL(1-OLx,1-OLy,bi,bj),
     I          frzSeaWat,
     I          iMin,iMax, jMin,jMax, bi, bj,
     I          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
#ifdef OLD_THSICE_CALL_SEQUENCE
      IF ( .TRUE. ) THEN
#else /* OLD_THSICE_CALL_SEQUENCE */
      IF ( thSIceAdvScheme.LE.0 ) THEN
C- note: 1) regarding sIceLoad in ocean-dynamics, in case thSIceAdvScheme > 0,
C          compute sIceLoad in THSICE_DO_ADVECT after seaice advection is done.
C        2) regarding sIceLoad in seaice-dynamics, probably better not to update
C          sIceLoad here, to keep the balance between sIceLoad and adjusted Eta.
C        3) not sure in the case of no advection (thSIceAdvScheme=0) but using
C          seaice dynamics (unlikely senario anyway).
#endif /* OLD_THSICE_CALL_SEQUENCE */
C--   Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
       DO j = jMin, jMax
        DO i = iMin, iMax
          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
      ENDIF

#ifdef OLD_THSICE_CALL_SEQUENCE
      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
#endif /* OLD_THSICE_CALL_SEQUENCE */

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