pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.13 2007/04/16 22:38:24 heimbach Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"

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 "GRID.h"
#include "SURFACE.h"
#include "DYNVARS.h"
#include "FFIELDS.h"
#include "THSICE_PARAMS.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
C--
# include "THSICE_2DYN.h"
C--
#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 ===
      INTEGER i,j
      INTEGER bi,bj
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      _RL prcAtm(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-|--+----|

      IF ( stressReduction.GT. 0. _d 0 ) THEN
C-     needs new Ice Fraction in halo region to apply wind-stress reduction
       iMin = 1-OLx
       iMax = sNx+OLx-1
       jMin = 1-OLy
       jMax = sNy+OLy-1
#ifdef ATMOSPHERIC_LOADING
      ELSEIF ( useRealFreshWaterFlux .AND. .NOT.useSEAICE ) THEN
C-     needs sea-ice loading in part of the halo regions for grad.Phi0surf
C      to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
       iMin = 0
       iMax = sNx+1
       jMin = 0
       jMax = sNy+1
#endif /* ATMOSPHERIC_LOADING */
      ELSE
       iMin = 1
       iMax = sNx
       jMin = 1
       jMax = sNy
      ENDIF

      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
          iicekey = (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=iicekey, byte=isbyte
CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
# ifdef ALLOW_EXF
CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
# endif
#endif

cph(
          print *, 'ph-thsice-1 in thsice_main'
cph)
C--     Mixed layer thickness: take the 1rst layer
#ifdef NONLIN_FRSURF
        IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
         IF ( select_rStar.GT.0 ) THEN
          DO j = jMin, jMax
           DO i = iMin, iMax
             hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
     &                                  *rStarFacC(i,j,bi,bj)
           ENDDO
          ENDDO
         ELSE
          DO j = jMin, jMax
           DO i = iMin, iMax
            IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
             hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
            ELSE
             hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
            ENDIF
           ENDDO
          ENDDO
         ENDIF
        ELSE
#else /* ndef NONLIN_FRSURF */
        IF (.TRUE.) THEN
#endif /* NONLIN_FRSURF */
          DO j = jMin, jMax
           DO i = iMin, iMax
             hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
           ENDDO
          ENDDO
        ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
#endif

        DO j = jMin, jMax
         DO i = iMin, iMax
          tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
          sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
          v2ocMxL(i,j,bi,bj) =
     &              ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
     &              + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
     &              + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
     &              + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
     &              )*0.5 _d 0
          prcAtm(i,j) = 0.
          icFrwAtm(i,j,bi,bj) = 0. _d 0
          icFlxAtm(i,j,bi,bj) = 0. _d 0
          icFlxSW (i,j,bi,bj) = 0. _d 0
          snowPrc(i,j,bi,bj) = 0. _d 0
          siceAlb(i,j,bi,bj) = 0. _d 0
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceMask = comlev1, key = iicekey
CADJ STORE iceHeight  = comlev1, key = iicekey
CADJ STORE snowHeight = comlev1, key = iicekey
CADJ STORE Tsrf    = comlev1, key = iicekey
CADJ STORE Qice1   = comlev1, key = iicekey
CADJ STORE Qice2   = comlev1, key = iicekey
CADJ STORE snowAge = comlev1, key = iicekey
CADJ STORE snowPrc  = comlev1, key = iicekey

CADJ STORE hOceMxL = comlev1, key = iicekey
CADJ STORE tOceMxL = comlev1, key = iicekey
CADJ STORE sOceMxL = comlev1, key = iicekey
CADJ STORE v2ocMxL = comlev1, key = iicekey

CADJ STORE empmr   = comlev1, key = iicekey
CADJ STORE qnet    = comlev1, key = iicekey
#endif

cph(
          print *, 'ph-thsice-2 in thsice_main'
cph)
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 )

#ifdef ALLOW_BULK_FORCE
        IF ( useBulkforce ) THEN
         CALL THSICE_GET_PRECIP(
     I                  iceMask,
     O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
     O                  icFlxSW(1-OLx,1-OLy,bi,bj),
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
        ENDIF
#endif
#ifdef ALLOW_EXF
        IF ( useEXF ) THEN
         CALL THSICE_MAP_EXF(
     I                  iceMask,
     O                  prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
     O                  icFlxSW(1-OLx,1-OLy,bi,bj),
     I                  iMin,iMax,jMin,jMax, bi,bj, myThid )
        ENDIF
#endif


cph(
          print *, 'ph-thsice-3 in thsice_main'
cph)
        CALL THSICE_STEP_TEMP(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     myTime, myIter, myThid )

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE empmr, qnet = comlev1, key = iicekey
CADJ STORE iceMask = comlev1, key = iicekey
CADJ STORE iceHeight  = comlev1, key = iicekey
CADJ STORE snowHeight = comlev1, key = iicekey
CADJ STORE Tsrf    = comlev1, key = iicekey
CADJ STORE Qice1   = comlev1, key = iicekey
CADJ STORE Qice2   = comlev1, key = iicekey
CADJ STORE snowAge = comlev1, key = iicekey
#endif

cph(
          print *, 'ph-thsice-4 in thsice_main'
cph)
        CALL THSICE_STEP_FWD(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     prcAtm,
     I                     myTime, myIter, myThid )

cph(
          print *, 'ph-thsice-5 in thsice_main'
cph)
        CALL THSICE_AVE(
     I                     bi,bj, myTime, myIter, myThid )

c      ENDDO
c     ENDDO

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=iicekey, byte=isbyte
CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
#endif
        IF ( stressReduction.GT. 0. _d 0 ) THEN
          DO j = jMin, jMax
           DO i = iMin+1,iMax
            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 = jMin+1, jMax
           DO i = iMin, iMax
            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


cph(
          print *, 'ph-thsice-6 in thsice_main'
cph)
      IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN
C--   Exchange fields that are advected by seaice dynamics
        _EXCH_XY_R8( iceMask, myThid )
        _EXCH_XY_R8( iceHeight, myThid )
        _EXCH_XY_R8( snowHeight, myThid )
        _EXCH_XY_R8( Qice1, myThid )
        _EXCH_XY_R8( Qice2, myThid )
#ifdef ATMOSPHERIC_LOADING
        IF (useRealFreshWaterFlux)
     &  _EXCH_XY_RS( sIceLoad, myThid )
#endif
      ENDIF

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

      RETURN
      END
1	heimbach	1.14	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.13 2007/04/16 22:38:24 heimbach Exp $
2	jmc	1.1	C $Name: $
3
4			#include "THSICE_OPTIONS.h"
5	jmc	1.8
6	jmc	1.2	CBOP
7	jmc	1.1	C !ROUTINE: THSICE_MAIN
8			C !INTERFACE:
9	jmc	1.8	SUBROUTINE THSICE_MAIN(
10	jmc	1.1	I myTime, myIter, myThid )
11	jmc	1.2	C !DESCRIPTION: \bv
12	jmc	1.1	C ==========================================================
13	jmc	1.8	C \| S/R THSICE_MAIN
14			C \| o Therm_SeaIce main routine.
15	jmc	1.1	C \| step forward Thermodynamic_SeaIce variables and modify
16			C \| ocean surface forcing accordingly.
17			C ==========================================================
18
19			C !USES:
20			IMPLICIT NONE
21	jmc	1.2
22	jmc	1.1	C === Global variables ===
23			#include "SIZE.h"
24			#include "EEPARAMS.h"
25			#include "PARAMS.h"
26	jmc	1.2	#include "GRID.h"
27	jmc	1.3	#include "SURFACE.h"
28	jmc	1.2	#include "DYNVARS.h"
29	jmc	1.1	#include "FFIELDS.h"
30			#include "THSICE_PARAMS.h"
31	jmc	1.2	#include "THSICE_VARS.h"
32	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
33			# include "tamc.h"
34			# include "tamc_keys.h"
35	heimbach	1.13	C--
36			# include "THSICE_2DYN.h"
37			C--
38	heimbach	1.7	#endif
39	jmc	1.8
40	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
41			C === Routine arguments ===
42	jmc	1.12	C myTime :: Current time in simulation (s)
43			C myIter :: Current iteration number
44			C myThid :: My Thread Id. number
45			_RL myTime
46	jmc	1.1	INTEGER myIter
47			INTEGER myThid
48	jmc	1.2	CEOP
49	jmc	1.1
50			#ifdef ALLOW_THSICE
51			C !LOCAL VARIABLES:
52			C === Local variables ===
53			INTEGER i,j
54			INTEGER bi,bj
55			INTEGER iMin, iMax
56			INTEGER jMin, jMax
57	jmc	1.2	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58	jmc	1.8	c _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59			c _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
60			c _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
61	jmc	1.1
62			_RL tauFac
63
64			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
65
66			IF ( stressReduction.GT. 0. _d 0 ) THEN
67	jmc	1.5	C- needs new Ice Fraction in halo region to apply wind-stress reduction
68	jmc	1.8	iMin = 1-OLx
69			iMax = sNx+OLx-1
70			jMin = 1-OLy
71			jMax = sNy+OLy-1
72	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
73	jmc	1.11	ELSEIF ( useRealFreshWaterFlux .AND. .NOT.useSEAICE ) THEN
74	jmc	1.5	C- needs sea-ice loading in part of the halo regions for grad.Phi0surf
75			C to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
76			iMin = 0
77			iMax = sNx+1
78			jMin = 0
79			jMax = sNy+1
80			#endif /* ATMOSPHERIC_LOADING */
81	jmc	1.1	ELSE
82			iMin = 1
83			iMax = sNx
84			jMin = 1
85			jMax = sNy
86			ENDIF
87
88			DO bj=myByLo(myThid),myByHi(myThid)
89			DO bi=myBxLo(myThid),myBxHi(myThid)
90	jmc	1.2
91	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
92			act1 = bi - myBxLo(myThid)
93			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
94			act2 = bj - myByLo(myThid)
95			max2 = myByHi(myThid) - myByLo(myThid) + 1
96			act3 = myThid - 1
97			max3 = nTx*nTy
98			act4 = ikey_dynamics - 1
99			iicekey = (act1 + 1) + act2*max1
100			& + act3max1max2
101			& + act4max1max2*max3
102			#endif /* ALLOW_AUTODIFF_TAMC */
103
104	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
105			CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
106			CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
107			CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
108	heimbach	1.14	# ifdef ALLOW_EXF
109			CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
110			# endif
111	heimbach	1.13	#endif
112
113	heimbach	1.14	cph(
114			print *, 'ph-thsice-1 in thsice_main'
115			cph)
116	jmc	1.3	C-- Mixed layer thickness: take the 1rst layer
117			#ifdef NONLIN_FRSURF
118			IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
119			IF ( select_rStar.GT.0 ) THEN
120			DO j = jMin, jMax
121			DO i = iMin, iMax
122			hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
123	jmc	1.8	& *rStarFacC(i,j,bi,bj)
124	jmc	1.3	ENDDO
125			ENDDO
126			ELSE
127			DO j = jMin, jMax
128			DO i = iMin, iMax
129			IF ( ksurfC(i,j,bi,bj).EQ.1 ) THEN
130			hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
131			ELSE
132	jmc	1.12	hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
133	jmc	1.3	ENDIF
134			ENDDO
135			ENDDO
136			ENDIF
137			ELSE
138			#else /* ndef NONLIN_FRSURF */
139			IF (.TRUE.) THEN
140			#endif /* NONLIN_FRSURF */
141			DO j = jMin, jMax
142			DO i = iMin, iMax
143	jmc	1.12	hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
144	jmc	1.3	ENDDO
145			ENDDO
146			ENDIF
147
148	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
149			CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
150			CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
151			#endif
152
153	mlosch	1.10	DO j = jMin, jMax
154			DO i = iMin, iMax
155			tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
156			sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
157			v2ocMxL(i,j,bi,bj) =
158	jmc	1.2	& ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
159			& + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
160			& + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
161			& + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
162			& )*0.5 _d 0
163	mlosch	1.10	prcAtm(i,j) = 0.
164			icFrwAtm(i,j,bi,bj) = 0. _d 0
165			icFlxAtm(i,j,bi,bj) = 0. _d 0
166			icFlxSW (i,j,bi,bj) = 0. _d 0
167			snowPrc(i,j,bi,bj) = 0. _d 0
168			siceAlb(i,j,bi,bj) = 0. _d 0
169	jmc	1.6	ENDDO
170	mlosch	1.10	ENDDO
171	heimbach	1.7
172			#ifdef ALLOW_AUTODIFF_TAMC
173			CADJ STORE iceMask = comlev1, key = iicekey
174			CADJ STORE iceHeight = comlev1, key = iicekey
175			CADJ STORE snowHeight = comlev1, key = iicekey
176			CADJ STORE Tsrf = comlev1, key = iicekey
177			CADJ STORE Qice1 = comlev1, key = iicekey
178			CADJ STORE Qice2 = comlev1, key = iicekey
179			CADJ STORE snowAge = comlev1, key = iicekey
180			CADJ STORE snowPrc = comlev1, key = iicekey
181
182			CADJ STORE hOceMxL = comlev1, key = iicekey
183			CADJ STORE tOceMxL = comlev1, key = iicekey
184			CADJ STORE sOceMxL = comlev1, key = iicekey
185			CADJ STORE v2ocMxL = comlev1, key = iicekey
186
187			CADJ STORE empmr = comlev1, key = iicekey
188			CADJ STORE qnet = comlev1, key = iicekey
189			#endif
190
191	heimbach	1.14	cph(
192			print *, 'ph-thsice-2 in thsice_main'
193			cph)
194	jmc	1.12	C- do sea-ice advection before getting surface fluxes
195			C Note: will inline this S/R once thSIce in Atmos. set-up is settled
196			IF ( thSIceAdvScheme.GT.0 )
197			& CALL THSICE_DO_ADVECT(
198			I bi,bj, myTime, myIter, myThid )
199
200	jmc	1.2	#ifdef ALLOW_BULK_FORCE
201	mlosch	1.10	IF ( useBulkforce ) THEN
202			CALL THSICE_GET_PRECIP(
203	jmc	1.6	I iceMask,
204	jmc	1.8	O prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
205			O icFlxSW(1-OLx,1-OLy,bi,bj),
206	jmc	1.6	I iMin,iMax,jMin,jMax, bi,bj, myThid )
207	mlosch	1.10	ENDIF
208	jmc	1.2	#endif
209	mlosch	1.9	#ifdef ALLOW_EXF
210	mlosch	1.10	IF ( useEXF ) THEN
211			CALL THSICE_MAP_EXF(
212	mlosch	1.9	I iceMask,
213			O prcAtm, snowPrc(1-OLx,1-OLy,bi,bj),
214			O icFlxSW(1-OLx,1-OLy,bi,bj),
215			I iMin,iMax,jMin,jMax, bi,bj, myThid )
216	mlosch	1.10	ENDIF
217	mlosch	1.9	#endif
218	jmc	1.2
219	heimbach	1.7
220	heimbach	1.14	cph(
221			print *, 'ph-thsice-3 in thsice_main'
222			cph)
223	mlosch	1.10	CALL THSICE_STEP_TEMP(
224	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
225			I myTime, myIter, myThid )
226
227	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
228			CADJ STORE empmr, qnet = comlev1, key = iicekey
229			CADJ STORE iceMask = comlev1, key = iicekey
230			CADJ STORE iceHeight = comlev1, key = iicekey
231			CADJ STORE snowHeight = comlev1, key = iicekey
232			CADJ STORE Tsrf = comlev1, key = iicekey
233			CADJ STORE Qice1 = comlev1, key = iicekey
234			CADJ STORE Qice2 = comlev1, key = iicekey
235			CADJ STORE snowAge = comlev1, key = iicekey
236			#endif
237
238	heimbach	1.14	cph(
239			print *, 'ph-thsice-4 in thsice_main'
240			cph)
241	mlosch	1.10	CALL THSICE_STEP_FWD(
242	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
243			I prcAtm,
244	jmc	1.1	I myTime, myIter, myThid )
245	jmc	1.2
246	heimbach	1.14	cph(
247			print *, 'ph-thsice-5 in thsice_main'
248			cph)
249	mlosch	1.10	CALL THSICE_AVE(
250	jmc	1.8	I bi,bj, myTime, myIter, myThid )
251	jmc	1.2
252	jmc	1.1	c ENDDO
253			c ENDDO
254
255	jmc	1.11	C-- note: If useSEAICE=.true., the stress is computed in seaice_model,
256			C-- and stressReduction is always set to zero
257	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
258			CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
259			CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=iicekey, byte=isbyte
260			#endif
261	jmc	1.11	IF ( stressReduction.GT. 0. _d 0 ) THEN
262	mlosch	1.10	DO j = jMin, jMax
263			DO i = iMin+1,iMax
264	jmc	1.1	tauFac = stressReduction
265			& (iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
266			fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
267	mlosch	1.10	ENDDO
268	jmc	1.1	ENDDO
269	mlosch	1.10	DO j = jMin+1, jMax
270			DO i = iMin, iMax
271	jmc	1.1	tauFac = stressReduction
272			& (iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
273			fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
274	mlosch	1.10	ENDDO
275	jmc	1.1	ENDDO
276			ENDIF
277
278			C-- end bi,bj loop
279			ENDDO
280			ENDDO
281
282	jmc	1.5
283	heimbach	1.14	cph(
284			print *, 'ph-thsice-6 in thsice_main'
285			cph)
286	jmc	1.12	IF ( useSEAICE .OR. thSIceAdvScheme.GT.0 ) THEN
287	mlosch	1.10	C-- Exchange fields that are advected by seaice dynamics
288	jmc	1.11	_EXCH_XY_R8( iceMask, myThid )
289			_EXCH_XY_R8( iceHeight, myThid )
290			_EXCH_XY_R8( snowHeight, myThid )
291			_EXCH_XY_R8( Qice1, myThid )
292			_EXCH_XY_R8( Qice2, myThid )
293			#ifdef ATMOSPHERIC_LOADING
294			IF (useRealFreshWaterFlux)
295			& _EXCH_XY_RS( sIceLoad, myThid )
296			#endif
297	mlosch	1.10	ENDIF
298	jmc	1.11
299	jmc	1.1	C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
300			#endif /ALLOW_THSICE/
301
302			RETURN
303			END