pkg/thsice/thsice_step_fwd.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.16 2006/04/09 17:35:30 heimbach Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"

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"
#include "THSICE_SIZE.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#include "THSICE_TAVE.h"
      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-|--+----|

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

      ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
      snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
      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_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

      CALL THSICE_CALC_THICKN(
     I          bi, bj, siLo, siHi, sjLo, sjHi,
     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 )

C--    Net fluxes :
      DO j = jMin, jMax
       DO i = iMin, iMax
        IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
C-     weighted average net fluxes:
          icFrac = iceMask(i,j,bi,bj)
          opFrac= 1. _d 0-icFrac
          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*rhofw*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)/rhofw
     &                    +  opFrac*EmPmR(i,j,bi,bj)
          saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)

#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) = rhofw*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, siLo, siHi, sjLo, sjHi,
     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 )

      DO j = jMin, jMax
       DO i = iMin, iMax
        IF (frzmltMxL(i,j).GT.0. _d 0) THEN
C--    Net fluxes :
          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)/rhofw
          saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)

#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)    = 0. _d 0
          Qice2(i,j,bi,bj)    = 0. _d 0
        ENDIF

#ifdef ATMOSPHERIC_LOADING
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)
#endif

       ENDDO
      ENDDO

#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.17	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_step_fwd.F,v 1.16 2006/04/09 17:35:30 heimbach Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	jmc	1.3
6			CBOP
7	jmc	1.1	C !ROUTINE: THSICE_STEP_FWD
8			C !INTERFACE:
9	jmc	1.13	SUBROUTINE THSICE_STEP_FWD(
10	jmc	1.1	I bi, bj, iMin, iMax, jMin, jMax,
11	jmc	1.3	I prcAtm,
12	jmc	1.1	I myTime, myIter, myThid )
13	jmc	1.3	C !DESCRIPTION: \bv
14	jmc	1.1	C ==========================================================
15	jmc	1.13	C \| S/R THSICE_STEP_FWD
16	jmc	1.1	C \| o Step Forward Therm-SeaIce model.
17			C ==========================================================
18	jmc	1.3	C \ev
19	jmc	1.1
20			C !USES:
21			IMPLICIT NONE
22	jmc	1.3
23	jmc	1.1	C === Global variables ===
24			#include "SIZE.h"
25			#include "EEPARAMS.h"
26			#include "PARAMS.h"
27			#include "FFIELDS.h"
28			#include "THSICE_SIZE.h"
29			#include "THSICE_PARAMS.h"
30	jmc	1.3	#include "THSICE_VARS.h"
31			#include "THSICE_TAVE.h"
32	jmc	1.17	INTEGER siLo, siHi, sjLo, sjHi
33			PARAMETER ( siLo = 1-OLx , siHi = sNx+OLx )
34			PARAMETER ( sjLo = 1-OLy , sjHi = sNy+OLy )
35	jmc	1.13
36	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
37			C === Routine arguments ===
38	jmc	1.17	C- input:
39	jmc	1.3	C bi,bj :: tile indices
40			C iMin,iMax :: computation domain: 1rst index range
41			C jMin,jMax :: computation domain: 2nd index range
42	jmc	1.17	C prcAtm :: total precip from the atmosphere [kg/m2/s]
43			C myTime :: current Time of simulation [s]
44			C myIter :: current Iteration number in simulation
45			C myThid :: my Thread Id number
46			C-- Use fluxes hold in commom blocks
47	jmc	1.3	C- input:
48	jmc	1.17	C icFlxSW :: net short-wave heat flux (+=down) below sea-ice, into ocean
49			C icFlxAtm :: net Atmospheric surf. heat flux over sea-ice [W/m2], (+=down)
50			C icFrwAtm :: evaporation over sea-ice to the atmosphere [kg/m2/s] (+=up)
51	jmc	1.3	C- output
52	jmc	1.17	C icFlxAtm :: net Atmospheric surf. heat flux over ice+ocean [W/m2], (+=down)
53			C icFrwAtm :: net fresh-water flux (E-P) from the atmosphere [m/s] (+=up)
54	jmc	1.1	INTEGER bi,bj
55			INTEGER iMin, iMax
56			INTEGER jMin, jMax
57	jmc	1.3	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58	jmc	1.1	_RL myTime
59			INTEGER myIter
60			INTEGER myThid
61	jmc	1.3	CEOP
62	jmc	1.1
63			#ifdef ALLOW_THSICE
64			C !LOCAL VARIABLES:
65			C === Local variables ===
66	jmc	1.17	C iceFrac :: fraction of grid area covered in ice
67	jmc	1.3	C flx2oc :: net heat flux from the ice to the ocean (+=down) [W/m2]
68			C frw2oc :: fresh-water flux from the ice to the ocean
69			C fsalt :: mass salt flux to the ocean
70			C frzmltMxL :: ocean mixed-layer freezing/melting potential [W/m2]
71	jmc	1.17	C tFrzOce :: sea-water freezing temperature [oC] (function of S)
72	jmc	1.14	C isIceFree :: true for ice-free grid-cell that remains ice-free
73	jmc	1.17	C ageFac :: snow aging factor [1]
74			C snowFac :: snowing refreshing-age factor [units of 1/snowPr]
75			LOGICAL isIceFree(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76			_RL iceFrac (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77			_RL flx2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78			_RL frw2oc (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79			_RL fsalt (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80			_RL tFrzOce (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81			_RL frzmltMxL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82			_RL ageFac
83			_RL snowFac
84			_RL cphm
85	jmc	1.3	_RL opFrac, icFrac
86	jmc	1.9	#ifdef ALLOW_DIAGNOSTICS
87	jmc	1.10	_RL tmpFac
88	jmc	1.9	#endif
89	jmc	1.17	INTEGER i,j
90			LOGICAL dBugFlag
91	jmc	1.1
92	jmc	1.17	C- define grid-point location where to print debugging values
93			#include "THSICE_DEBUG.h"
94	jmc	1.1
95	jmc	1.13	1010 FORMAT(A,1P4E14.6)
96	jmc	1.17
97			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
98
99			C- Initialise
100			dBugFlag = debugLevel.GE.debLevB
101			DO j = 1-OLy, sNy+OLy
102			DO i = 1-OLx, sNx+OLx
103	jmc	1.14	isIceFree(i,j) = .FALSE.
104	jmc	1.17	saltFlux(i,j,bi,bj) = 0. _d 0
105	heimbach	1.16	#ifdef ALLOW_AUTODIFF_TAMC
106	jmc	1.17	iceFrac(i,j) = 0.
107	heimbach	1.16	#endif
108	jmc	1.4	ENDDO
109	jmc	1.17	ENDDO
110	jmc	1.1
111	jmc	1.17	ageFac = 1. _d 0 - thSIce_deltaT/snowAgTime
112			snowFac = thSIce_deltaT/(rhos*hNewSnowAge)
113			DO j = jMin, jMax
114			DO i = iMin, iMax
115			IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
116			C-- Snow aging :
117			snowAge(i,j,bi,bj) = thSIce_deltaT
118			& + snowAge(i,j,bi,bj)*ageFac
119			IF ( snowPrc(i,j,bi,bj).GT.0. _d 0 )
120			& snowAge(i,j,bi,bj) = snowAge(i,j,bi,bj)
121			& * EXP( - snowFac*snowPrc(i,j,bi,bj) )
122			c & * EXP( -(thSIce_deltaT*snowPrc(i,j,bi,bj)/rhos)
123			c & /hNewSnowAge )
124	jmc	1.3	C-------
125	jmc	1.17	C note: Any flux of mass (here fresh water) that enter or leave the system
126			C with a non zero energy HAS TO be counted: add snow precip.
127			icFlxAtm(i,j,bi,bj) = icFlxAtm(i,j,bi,bj)
128			& - Lfresh*snowPrc(i,j,bi,bj)
129			C--
130			ENDIF
131	jmc	1.3	ENDDO
132	jmc	1.17	ENDDO
133	jmc	1.3
134	jmc	1.9	#ifdef ALLOW_DIAGNOSTICS
135			IF ( useDiagnostics ) THEN
136	jmc	1.10	tmpFac = 1. _d 0
137			CALL DIAGNOSTICS_FRACT_FILL(
138	jmc	1.11	I snowPrc, iceMask,tmpFac,1,'SIsnwPrc',
139	jmc	1.10	I 0,1,1,bi,bj,myThid)
140	jmc	1.12	CALL DIAGNOSTICS_FRACT_FILL(
141			I siceAlb, iceMask,tmpFac,1,'SIalbedo',
142			I 0,1,1,bi,bj,myThid)
143	jmc	1.9	ENDIF
144			#endif /* ALLOW_DIAGNOSTICS */
145	jmc	1.12	DO j = jMin, jMax
146			DO i = iMin, iMax
147			siceAlb(i,j,bi,bj) = iceMask(i,j,bi,bj)*siceAlb(i,j,bi,bj)
148			ENDDO
149			ENDDO
150	jmc	1.9
151	jmc	1.3	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
152	jmc	1.13	C part.2 : ice-covered fraction ;
153			C change in ice/snow thickness and ice-fraction
154	jmc	1.3	C note: can only reduce the ice-fraction but not increase it.
155			C-------
156	jmc	1.1	DO j = jMin, jMax
157			DO i = iMin, iMax
158	jmc	1.3
159	jmc	1.17	tFrzOce(i,j) = -mu_Tf*sOceMxL(i,j,bi,bj)
160	jmc	1.3	cphm = cpwaterrhoswhOceMxL(i,j,bi,bj)
161	jmc	1.17	frzmltMxL(i,j) = ( tFrzOce(i,j)-tOceMxL(i,j,bi,bj) )
162			& * cphm/ocean_deltaT
163			iceFrac(i,j) = iceMask(i,j,bi,bj)
164			flx2oc(i,j) = icFlxSW(i,j,bi,bj)
165	jmc	1.3	C-------
166	jmc	1.17	#ifdef ALLOW_DBUG_THSICE
167			IF ( dBug(i,j,bi,bj) ) THEN
168			IF (frzmltMxL(i,j).GT.0. .OR. iceFrac(i,j).GT.0.) THEN
169	jmc	1.5	WRITE(6,'(A,2I4,2I2)') 'ThSI_FWD: i,j=',i,j,bi,bj
170			WRITE(6,1010) 'ThSI_FWD:-1- iceMask, hIc, hSn, Tsf =',
171	jmc	1.17	& iceFrac(i,j), iceHeight(i,j,bi,bj),
172	jmc	1.7	& snowHeight(i,j,bi,bj), Tsrf(i,j,bi,bj)
173	jmc	1.17	WRITE(6,1010) 'ThSI_FWD: ocTs,tFrzOce,frzmltMxL,Qnet=',
174			& tOceMxL(i,j,bi,bj), tFrzOce(i,j),
175			& frzmltMxL(i,j), Qnet(i,j,bi,bj)
176			ENDIF
177			IF (iceFrac(i,j).GT.0.)
178			& WRITE(6,1010) 'ThSI_FWD: icFrac,flxAtm,evpAtm,flxSnw=',
179			& iceFrac(i,j), icFlxAtm(i,j,bi,bj),
180			& icFrwAtm(i,j,bi,bj),-Lfresh*snowPrc(i,j,bi,bj)
181	jmc	1.1	ENDIF
182	jmc	1.17	#endif
183			ENDDO
184			ENDDO
185	jmc	1.1
186	jmc	1.17	CALL THSICE_CALC_THICKN(
187			I bi, bj, siLo, siHi, sjLo, sjHi,
188			I iMin,iMax, jMin,jMax, dBugFlag,
189			I iceMask(siLo,sjLo,bi,bj), tFrzOce,
190			I tOceMxL(siLo,sjLo,bi,bj), v2ocMxL(siLo,sjLo,bi,bj),
191			I snowPrc(siLo,sjLo,bi,bj), prcAtm,
192			I sHeating(siLo,sjLo,bi,bj), flxCndBt(siLo,sjLo,bi,bj),
193			U iceFrac, iceHeight(siLo,sjLo,bi,bj),
194			U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
195			U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
196			U icFrwAtm(siLo,sjLo,bi,bj), frzmltMxL, flx2oc,
197			O frw2oc, fsalt,
198			I myTime, myIter, myThid )
199	jmc	1.1
200			C-- Net fluxes :
201	jmc	1.17	DO j = jMin, jMax
202			DO i = iMin, iMax
203			IF (iceMask(i,j,bi,bj).GT.0. _d 0) THEN
204	jmc	1.3	C- weighted average net fluxes:
205			icFrac = iceMask(i,j,bi,bj)
206			opFrac= 1. _d 0-icFrac
207	jmc	1.17	icFlxAtm(i,j,bi,bj) = icFrac*icFlxAtm(i,j,bi,bj)
208			& - opFrac*Qnet(i,j,bi,bj)
209			icFrwAtm(i,j,bi,bj) = icFrac*icFrwAtm(i,j,bi,bj)
210			& + opFracrhofwEmPmR(i,j,bi,bj)
211			Qnet(i,j,bi,bj) = -icFracflx2oc(i,j) + opFracQnet(i,j,bi,bj)
212			EmPmR(i,j,bi,bj)= -icFrac*frw2oc(i,j)/rhofw
213			& + opFrac*EmPmR(i,j,bi,bj)
214			saltFlux(i,j,bi,bj) = -icFrac*fsalt(i,j)
215
216			#ifdef ALLOW_DBUG_THSICE
217			IF (dBug(i,j,bi,bj)) WRITE(6,1010)
218			& 'ThSI_FWD:-3- iceFrac, hIc, hSn, Qnet =',
219			& iceFrac(i,j), iceHeight(i,j,bi,bj),
220			& snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
221			#endif
222	jmc	1.1
223	jmc	1.3	ELSEIF (hOceMxL(i,j,bi,bj).gt.0. _d 0) THEN
224	jmc	1.17	icFlxAtm(i,j,bi,bj) = -Qnet(i,j,bi,bj)
225			icFrwAtm(i,j,bi,bj) = rhofw*EmPmR(i,j,bi,bj)
226	jmc	1.3	ELSE
227	jmc	1.17	icFlxAtm(i,j,bi,bj) = 0. _d 0
228			icFrwAtm(i,j,bi,bj) = 0. _d 0
229	jmc	1.1	ENDIF
230	jmc	1.17	ENDDO
231			ENDDO
232	jmc	1.1
233			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
234	jmc	1.3	C part.3 : freezing of sea-water
235	jmc	1.1	C over ice-free fraction and what is left from ice-covered fraction
236			C-------
237	jmc	1.17	CALL THSICE_EXTEND(
238			I bi, bj, siLo, siHi, sjLo, sjHi,
239			I iMin,iMax, jMin,jMax, dBugFlag,
240			I frzmltMxL, tFrzOce,
241			I tOceMxL(siLo,sjLo,bi,bj),
242			U iceFrac, iceHeight(siLo,sjLo,bi,bj),
243			U snowHeight(siLo,sjLo,bi,bj), Tsrf(siLo,sjLo,bi,bj),
244			U Tice1(siLo,sjLo,bi,bj), Tice2(siLo,sjLo,bi,bj),
245			U Qice1(siLo,sjLo,bi,bj), Qice2(siLo,sjLo,bi,bj),
246			O flx2oc, frw2oc, fsalt,
247			I myTime, myIter, myThid )
248
249			DO j = jMin, jMax
250			DO i = iMin, iMax
251			IF (frzmltMxL(i,j).GT.0. _d 0) THEN
252	jmc	1.1	C-- Net fluxes :
253	jmc	1.17	Qnet(i,j,bi,bj) = Qnet(i,j,bi,bj) - flx2oc(i,j)
254			EmPmR(i,j,bi,bj)= EmPmR(i,j,bi,bj)- frw2oc(i,j)/rhofw
255			saltFlux(i,j,bi,bj)=saltFlux(i,j,bi,bj) - fsalt(i,j)
256
257			#ifdef ALLOW_DBUG_THSICE
258			IF (dBug(i,j,bi,bj)) WRITE(6,1010)
259			& 'ThSI_FWD:-4- iceFrac, hIc, hSn, Qnet =',
260			& iceFrac(i,j), iceHeight(i,j,bi,bj),
261			& snowHeight(i,j,bi,bj), Qnet(i,j,bi,bj)
262			#endif
263	jmc	1.1	ENDIF
264
265	jmc	1.14	IF ( hOceMxL(i,j,bi,bj).GT.0. _d 0 )
266			& isIceFree(i,j) = iceMask(i,j,bi,bj).LE.0. _d 0
267	jmc	1.17	& .AND. iceFrac(i,j) .LE.0. _d 0
268			IF ( iceFrac(i,j) .GT. 0. _d 0 ) THEN
269			iceMask(i,j,bi,bj)=iceFrac(i,j)
270			IF ( snowHeight(i,j,bi,bj).EQ.0. _d 0 )
271			& snowAge(i,j,bi,bj) = 0. _d 0
272	jmc	1.1	ELSE
273			iceMask(i,j,bi,bj) = 0. _d 0
274			iceHeight(i,j,bi,bj)= 0. _d 0
275			snowHeight(i,j,bi,bj)=0. _d 0
276	jmc	1.3	snowAge(i,j,bi,bj) = 0. _d 0
277	jmc	1.17	Tsrf(i,j,bi,bj) = tOceMxL(i,j,bi,bj)
278	jmc	1.1	Tice1(i,j,bi,bj) = 0. _d 0
279			Tice2(i,j,bi,bj) = 0. _d 0
280			Qice1(i,j,bi,bj) = 0. _d 0
281			Qice2(i,j,bi,bj) = 0. _d 0
282			ENDIF
283
284	jmc	1.6	#ifdef ATMOSPHERIC_LOADING
285			C-- Compute Sea-Ice Loading (= mass of sea-ice + snow / area unit)
286			sIceLoad(i,j,bi,bj) = ( snowHeight(i,j,bi,bj)*rhos
287			& + iceHeight(i,j,bi,bj)*rhoi
288			& )*iceMask(i,j,bi,bj)
289			#endif
290
291	jmc	1.1	ENDDO
292			ENDDO
293
294	jmc	1.14	#ifdef ALLOW_BULK_FORCE
295			IF ( useBulkForce ) THEN
296			CALL BULKF_FLUX_ADJUST(
297			I bi, bj, iMin, iMax, jMin, jMax,
298			I isIceFree, myTime, myIter, myThid )
299			ENDIF
300			#endif /* ALLOW_BULK_FORCE */
301
302	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
303			#endif /* ALLOW_THSICE */
304
305			RETURN
306			END