pkg/thsice/thsice_step_fwd.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.28 2010/10/16 19:22:34 heimbach 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,
     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
      iicekey = (act1 + 1) + act2*max1
     &                     + act3*max1*max2
     &                     + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

C-    Initialise
      dBugFlag = debugLevel.GE.debLevB
      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.
#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=iicekey,byte=isbyte
CADJ STORE iceheight(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE icfrwatm(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE qice1(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE qice2(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE snowheight(:,:,bi,bj) = comlev1_bibj,key=iicekey,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    &          * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
c    &                  /hNewSnowAge )
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=iicekey,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=iicekey,byte=isbyte
CADJ STORE fsalt(:,:)  = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE flx2oc(:,:) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE frw2oc(:,:) = comlev1_bibj,key=iicekey,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

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