pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.37 2013/05/02 20:11:06 jmc Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# ifdef ALLOW_EXF
#  include "EXF_OPTIONS.h"
# endif
#endif

CBOP
C     !ROUTINE: THSICE_MAIN
C     !INTERFACE:
      SUBROUTINE THSICE_MAIN(
     I                        myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | S/R  THSICE_MAIN
C     | o Therm_SeaIce main routine.
C     |   step forward Thermodynamic_SeaIce variables and modify
C     |    ocean surface forcing accordingly.
C     *==========================================================*

C     !USES:
      IMPLICIT NONE

C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "THSICE_PARAMS.h"
#include "THSICE_SIZE.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "THSICE_TAVE.h"
# include "THSICE_COST.h"
# include "DYNVARS.h"
# include "tamc.h"
# include "tamc_keys.h"
# ifdef ALLOW_EXF
#  include "EXF_FIELDS.h"
#  include "EXF_PARAM.h"
#  include "EXF_CONSTANTS.h"
# endif /* ALLOW_EXF */
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime    :: Current time in simulation (s)
C     myIter    :: Current iteration number
C     myThid    :: My Thread Id. number
      _RL     myTime
      INTEGER myIter
      INTEGER myThid
CEOP

#ifdef ALLOW_THSICE
C     !LOCAL VARIABLES:
C     === Local variables ===
C     prcAtm    :: total precip from the atmosphere [kg/m2/s]
C     snowPr    :: snow precipitation               [kg/m2/s]
C     qPrcRn    :: Energy content of Precip+RunOff (+=down) [W/m2]
      INTEGER i,j
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      _RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL snowPr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
c     _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tauFac

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

C-    Only compute/update seaice fields over the interior
C     (excluding overlap) and apply exchanges when needed
      iMin = 1
      iMax = sNx
      jMin = 1
      jMax = sNy

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#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 */

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE qsw(:,:,bi,bj)     = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
#endif

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          prcAtm  (i,j,bi,bj) = 0. _d 0
          snowPr  (i,j) = 0. _d 0
          qPrcRn  (i,j) = 0. _d 0
         ENDDO
        ENDDO

#ifndef ALLOW_AUTODIFF_TAMC
        IF ( .NOT.useCheapAML ) THEN
#endif
         CALL THSICE_GET_OCEAN(
     I                          bi, bj, myTime, myIter, myThid )

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey

CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
#else
C-   end if not useCheapAML
        ENDIF
#endif

#ifdef OLD_THSICE_CALL_SEQUENCE
C-      do sea-ice advection before getting surface fluxes
C Note: will inline this S/R once thSIce in Atmos. set-up is settled
        IF ( thSIceAdvScheme.GT.0 )
     &   CALL THSICE_DO_ADVECT(
     I                   bi,bj, myTime, myIter, myThid )
#endif /* OLD_THSICE_CALL_SEQUENCE */

#ifndef ALLOW_AUTODIFF_TAMC
        IF ( useBulkforce .OR. useCheapAML ) THEN
         CALL THSICE_GET_PRECIP(
     I                  iceMask, tOceMxL,
     O                  prcAtm(1-OLx,1-OLy,bi,bj),
     O                  snowPr, qPrcRn,
     O                  icFlxSW(1-OLx,1-OLy,bi,bj),
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
        ENDIF
#endif
        IF ( useEXF ) THEN
         CALL THSICE_MAP_EXF(
     I                  iceMask, tOceMxL,
     O                  prcAtm(1-OLx,1-OLy,bi,bj),
     O                  snowPr, qPrcRn,
     O                  icFlxSW(1-OLx,1-OLy,bi,bj),
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
        ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
#else
        IF ( .NOT.( useCheapAML .OR. thSIce_skipThermo ) ) THEN
#endif
         CALL THSICE_STEP_TEMP(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     myTime, myIter, myThid )

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE iceHeight(:,:,bi,bj)  = comlev1_bibj, key = ticekey
CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
cphCADJ STORE Tsrf(:,:,bi,bj)    = comlev1_bibj, key = ticekey
CADJ STORE Qice1(:,:,bi,bj)   = comlev1_bibj, key = ticekey
CADJ STORE Qice2(:,:,bi,bj)   = comlev1_bibj, key = ticekey
CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
#else
C-   end if not skipThermo / useCheapAML
        ENDIF
        IF ( .NOT.thSIce_skipThermo ) THEN
#endif
         CALL THSICE_STEP_FWD(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     prcAtm(1-OLx,1-OLy,bi,bj),
     I                     snowPr, qPrcRn,
     I                     myTime, myIter, myThid )
#ifndef ALLOW_AUTODIFF_TAMC
        ENDIF
#endif

C--  end bi,bj loop
       ENDDO
      ENDDO

#ifdef ALLOW_BALANCE_FLUXES
C--   Balance net Fresh-Water flux from Atm+Land
      IF ( thSIceBalanceAtmFW.NE.0 ) THEN
        CALL THSICE_BALANCE_FRW(
     I                      iMin, iMax, jMin, jMax,
     I                      prcAtm, myTime, myIter, myThid )
      ENDIF
#endif

C     add a small piece of code to check AddFluid implementation:
c#include "thsice_test_addfluid.h"

C--   Exchange fields that are advected by seaice dynamics
      IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
     &               .OR. stressReduction.GT.zeroRL ) THEN
        CALL THSICE_DO_EXCH( myThid )
      ENDIF
#ifdef OLD_THSICE_CALL_SEQUENCE
#ifdef ATMOSPHERIC_LOADING
      IF ( useRealFreshWaterFlux )
     &  _EXCH_XY_RS( sIceLoad, myThid )
#endif
#else /* OLD_THSICE_CALL_SEQUENCE */
#ifdef ATMOSPHERIC_LOADING
      IF ( useRealFreshWaterFlux .AND. thSIceAdvScheme.LE.0 )
     &  _EXCH_XY_RS( sIceLoad, myThid )
#endif

C-    when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
C     otherwise, call it from here, after thsice-thermodynamics is done
      IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
         CALL THSICE_DO_ADVECT(
     I                          0, 0, myTime, myIter, myThid )
      ENDIF
#endif /* OLD_THSICE_CALL_SEQUENCE */

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
C--   Cumulate time-averaged fields and also fill-up flux diagnostics
C     (if not done in THSICE_DO_ADVECT call)
#ifdef OLD_THSICE_CALL_SEQUENCE
        IF ( .TRUE. ) THEN
#else /* OLD_THSICE_CALL_SEQUENCE */
        IF ( thSIceAdvScheme.LE.0 ) THEN
#endif /* OLD_THSICE_CALL_SEQUENCE */
         CALL THSICE_AVE(
     I                     bi,bj, myTime, myIter, myThid )
        ENDIF
C--   note: If useSEAICE=.true., the stress is computed in seaice_model,
C--   and stressReduction is always set to zero
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
#endif
        IF ( stressReduction.GT. 0. _d 0 ) THEN
          DO j = 1-OLy,sNy+OLy-1
           DO i = 2-OLx,sNx+OLx-1
            tauFac = stressReduction
     &             *(iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
            fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
           ENDDO
          ENDDO
          DO j = 2-OLy,sNy+OLy-1
           DO i = 1-OLx,sNx+OLx-1
            tauFac = stressReduction
     &             *(iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))*0.5 _d 0
            fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
           ENDDO
          ENDDO
        ENDIF

C--  end bi,bj loop
       ENDDO
      ENDDO

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

      RETURN
      END
1	jmc	1.38	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.37 2013/05/02 20:11:06 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	heimbach	1.32	#ifdef ALLOW_AUTODIFF_TAMC
6			# ifdef ALLOW_EXF
7			# include "EXF_OPTIONS.h"
8			# endif
9			#endif
10	jmc	1.8
11	jmc	1.2	CBOP
12	jmc	1.1	C !ROUTINE: THSICE_MAIN
13			C !INTERFACE:
14	jmc	1.8	SUBROUTINE THSICE_MAIN(
15	jmc	1.1	I myTime, myIter, myThid )
16	jmc	1.2	C !DESCRIPTION: \bv
17	jmc	1.1	C ==========================================================
18	jmc	1.8	C \| S/R THSICE_MAIN
19			C \| o Therm_SeaIce main routine.
20	jmc	1.1	C \| step forward Thermodynamic_SeaIce variables and modify
21			C \| ocean surface forcing accordingly.
22			C ==========================================================
23
24			C !USES:
25			IMPLICIT NONE
26	jmc	1.2
27	jmc	1.1	C === Global variables ===
28			#include "SIZE.h"
29			#include "EEPARAMS.h"
30			#include "PARAMS.h"
31			#include "FFIELDS.h"
32			#include "THSICE_PARAMS.h"
33	gforget	1.23	#include "THSICE_SIZE.h"
34	jmc	1.2	#include "THSICE_VARS.h"
35	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
36	jmc	1.33	# include "THSICE_TAVE.h"
37			# include "THSICE_COST.h"
38	jmc	1.29	# include "DYNVARS.h"
39	heimbach	1.7	# include "tamc.h"
40			# include "tamc_keys.h"
41	heimbach	1.32	# ifdef ALLOW_EXF
42			# include "EXF_FIELDS.h"
43			# include "EXF_PARAM.h"
44			# include "EXF_CONSTANTS.h"
45			# endif /* ALLOW_EXF */
46	heimbach	1.7	#endif
47	jmc	1.8
48	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
49			C === Routine arguments ===
50	jmc	1.12	C myTime :: Current time in simulation (s)
51			C myIter :: Current iteration number
52			C myThid :: My Thread Id. number
53			_RL myTime
54	jmc	1.1	INTEGER myIter
55			INTEGER myThid
56	jmc	1.2	CEOP
57	jmc	1.1
58			#ifdef ALLOW_THSICE
59			C !LOCAL VARIABLES:
60			C === Local variables ===
61	jmc	1.37	C prcAtm :: total precip from the atmosphere [kg/m2/s]
62			C snowPr :: snow precipitation [kg/m2/s]
63	jmc	1.31	C qPrcRn :: Energy content of Precip+RunOff (+=down) [W/m2]
64	jmc	1.1	INTEGER i,j
65			INTEGER bi,bj
66			INTEGER iMin, iMax
67			INTEGER jMin, jMax
68	jmc	1.25	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
69	jmc	1.36	_RL snowPr(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
70	jmc	1.31	_RL qPrcRn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
71	jmc	1.8	c _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72			c _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73			c _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74	jmc	1.1	_RL tauFac
75
76			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
77
78	jmc	1.34	C- Only compute/update seaice fields over the interior
79			C (excluding overlap) and apply exchanges when needed
80			iMin = 1
81			iMax = sNx
82			jMin = 1
83			jMax = sNy
84	jmc	1.1
85			DO bj=myByLo(myThid),myByHi(myThid)
86			DO bi=myBxLo(myThid),myBxHi(myThid)
87	jmc	1.2
88	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
89			act1 = bi - myBxLo(myThid)
90			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
91			act2 = bj - myByLo(myThid)
92			max2 = myByHi(myThid) - myByLo(myThid) + 1
93			act3 = myThid - 1
94			max3 = nTx*nTy
95			act4 = ikey_dynamics - 1
96	gforget	1.23	ticekey = (act1 + 1) + act2*max1
97	heimbach	1.7	& + act3max1max2
98			& + act4max1max2*max3
99			#endif /* ALLOW_AUTODIFF_TAMC */
100
101	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
102	gforget	1.23	CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
103			CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
104			CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
105	jmc	1.38	CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
106	gforget	1.23	CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
107			CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
108	heimbach	1.7	#endif
109
110	jmc	1.29	DO j=1-OLy,sNy+OLy
111			DO i=1-OLx,sNx+OLx
112	jmc	1.25	prcAtm (i,j,bi,bj) = 0. _d 0
113	jmc	1.36	snowPr (i,j) = 0. _d 0
114	jmc	1.31	qPrcRn (i,j) = 0. _d 0
115	jmc	1.6	ENDDO
116	mlosch	1.10	ENDDO
117	heimbach	1.7
118	jmc	1.38	#ifndef ALLOW_AUTODIFF_TAMC
119	jmc	1.33	IF ( .NOT.useCheapAML ) THEN
120			#endif
121	jmc	1.38	CALL THSICE_GET_OCEAN(
122			I bi, bj, myTime, myIter, myThid )
123	jmc	1.29
124	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
125	gforget	1.23	CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
126			CADJ STORE iceHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
127			CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
128			CADJ STORE Tsrf(:,:,bi,bj) = comlev1_bibj, key = ticekey
129			CADJ STORE Qice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
130			CADJ STORE Qice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
131			CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
132
133			CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
134			CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
135			CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
136			CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
137	jmc	1.38	#else
138			C- end if not useCheapAML
139			ENDIF
140	heimbach	1.7	#endif
141
142	jmc	1.28	#ifdef OLD_THSICE_CALL_SEQUENCE
143	jmc	1.12	C- do sea-ice advection before getting surface fluxes
144			C Note: will inline this S/R once thSIce in Atmos. set-up is settled
145			IF ( thSIceAdvScheme.GT.0 )
146			& CALL THSICE_DO_ADVECT(
147			I bi,bj, myTime, myIter, myThid )
148	jmc	1.28	#endif /* OLD_THSICE_CALL_SEQUENCE */
149	jmc	1.12
150	jmc	1.38	#ifndef ALLOW_AUTODIFF_TAMC
151			IF ( useBulkforce .OR. useCheapAML ) THEN
152	mlosch	1.10	CALL THSICE_GET_PRECIP(
153	jmc	1.35	I iceMask, tOceMxL,
154	jmc	1.25	O prcAtm(1-OLx,1-OLy,bi,bj),
155	jmc	1.36	O snowPr, qPrcRn,
156	jmc	1.8	O icFlxSW(1-OLx,1-OLy,bi,bj),
157	jmc	1.6	I iMin,iMax,jMin,jMax, bi,bj, myThid )
158	mlosch	1.10	ENDIF
159	jmc	1.2	#endif
160	mlosch	1.10	IF ( useEXF ) THEN
161			CALL THSICE_MAP_EXF(
162	jmc	1.31	I iceMask, tOceMxL,
163	jmc	1.25	O prcAtm(1-OLx,1-OLy,bi,bj),
164	jmc	1.36	O snowPr, qPrcRn,
165	mlosch	1.9	O icFlxSW(1-OLx,1-OLy,bi,bj),
166			I iMin,iMax,jMin,jMax, bi,bj, myThid )
167	mlosch	1.10	ENDIF
168	jmc	1.2
169	heimbach	1.21	#ifdef ALLOW_AUTODIFF_TAMC
170	jmc	1.26	CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
171	gforget	1.23	CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
172			CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
173	jmc	1.26	#else
174	jmc	1.38	IF ( .NOT.( useCheapAML .OR. thSIce_skipThermo ) ) THEN
175	heimbach	1.21	#endif
176	jmc	1.38	CALL THSICE_STEP_TEMP(
177	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
178			I myTime, myIter, myThid )
179
180	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
181	gforget	1.23	CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
182			CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
183			CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
184			CADJ STORE iceHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
185			CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
186			cphCADJ STORE Tsrf(:,:,bi,bj) = comlev1_bibj, key = ticekey
187			CADJ STORE Qice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
188			CADJ STORE Qice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
189			CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
190	jmc	1.26	CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
191	jmc	1.29	#else
192	jmc	1.38	C- end if not skipThermo / useCheapAML
193	jmc	1.29	ENDIF
194			IF ( .NOT.thSIce_skipThermo ) THEN
195	heimbach	1.13	#endif
196	jmc	1.38	CALL THSICE_STEP_FWD(
197	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
198	jmc	1.36	I prcAtm(1-OLx,1-OLy,bi,bj),
199			I snowPr, qPrcRn,
200	jmc	1.1	I myTime, myIter, myThid )
201	jmc	1.26	#ifndef ALLOW_AUTODIFF_TAMC
202			ENDIF
203			#endif
204	jmc	1.2
205	jmc	1.25	C-- end bi,bj loop
206			ENDDO
207			ENDDO
208
209			#ifdef ALLOW_BALANCE_FLUXES
210			C-- Balance net Fresh-Water flux from Atm+Land
211			IF ( thSIceBalanceAtmFW.NE.0 ) THEN
212	jmc	1.26	CALL THSICE_BALANCE_FRW(
213	jmc	1.25	I iMin, iMax, jMin, jMax,
214			I prcAtm, myTime, myIter, myThid )
215			ENDIF
216			#endif
217
218	jmc	1.24	C add a small piece of code to check AddFluid implementation:
219			c#include "thsice_test_addfluid.h"
220
221			C-- Exchange fields that are advected by seaice dynamics
222	jmc	1.28	IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
223	jmc	1.34	& .OR. stressReduction.GT.zeroRL ) THEN
224	jmc	1.28	CALL THSICE_DO_EXCH( myThid )
225	jmc	1.24	ENDIF
226	jmc	1.28	#ifdef OLD_THSICE_CALL_SEQUENCE
227	jmc	1.24	#ifdef ATMOSPHERIC_LOADING
228	jmc	1.34	IF ( useRealFreshWaterFlux )
229	jmc	1.24	& _EXCH_XY_RS( sIceLoad, myThid )
230			#endif
231	jmc	1.28	#else /* OLD_THSICE_CALL_SEQUENCE */
232			#ifdef ATMOSPHERIC_LOADING
233	jmc	1.34	IF ( useRealFreshWaterFlux .AND. thSIceAdvScheme.LE.0 )
234	jmc	1.28	& _EXCH_XY_RS( sIceLoad, myThid )
235			#endif
236
237			C- when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
238			C otherwise, call it from here, after thsice-thermodynamics is done
239			IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
240			CALL THSICE_DO_ADVECT(
241			I 0, 0, myTime, myIter, myThid )
242			ENDIF
243			#endif /* OLD_THSICE_CALL_SEQUENCE */
244	jmc	1.24
245			DO bj=myByLo(myThid),myByHi(myThid)
246			DO bi=myBxLo(myThid),myBxHi(myThid)
247	jmc	1.30	C-- Cumulate time-averaged fields and also fill-up flux diagnostics
248			C (if not done in THSICE_DO_ADVECT call)
249			#ifdef OLD_THSICE_CALL_SEQUENCE
250			IF ( .TRUE. ) THEN
251			#else /* OLD_THSICE_CALL_SEQUENCE */
252			IF ( thSIceAdvScheme.LE.0 ) THEN
253			#endif /* OLD_THSICE_CALL_SEQUENCE */
254			CALL THSICE_AVE(
255			I bi,bj, myTime, myIter, myThid )
256			ENDIF
257	jmc	1.11	C-- note: If useSEAICE=.true., the stress is computed in seaice_model,
258			C-- and stressReduction is always set to zero
259	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
260	gforget	1.23	CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
261			CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
262	heimbach	1.7	#endif
263	jmc	1.11	IF ( stressReduction.GT. 0. _d 0 ) THEN
264	jmc	1.34	DO j = 1-OLy,sNy+OLy-1
265			DO i = 2-OLx,sNx+OLx-1
266	jmc	1.1	tauFac = stressReduction
267			& (iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
268			fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
269	mlosch	1.10	ENDDO
270	jmc	1.1	ENDDO
271	jmc	1.34	DO j = 2-OLy,sNy+OLy-1
272			DO i = 1-OLx,sNx+OLx-1
273	jmc	1.1	tauFac = stressReduction
274			& (iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
275			fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
276	mlosch	1.10	ENDDO
277	jmc	1.1	ENDDO
278			ENDIF
279
280			C-- end bi,bj loop
281			ENDDO
282			ENDDO
283
284			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
285			#endif /ALLOW_THSICE/
286
287			RETURN
288			END