pkg/thsice/thsice_main.F

C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.26 2012/11/21 01:53:34 jmc 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_SIZE.h"
#include "THSICE_VARS.h"
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#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,nSx,nSy)
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 ( useEXF .OR. useSEAICE ) THEN
C-    EXF does not provide valid fields in overlap
       iMin = 1
       iMax = sNx
       jMin = 1
       jMax = sNy
      ELSEIF ( 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 ) 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
          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
# ifdef ALLOW_EXF
CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
# endif
#endif

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=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
     &                                  *rStarFacC(i,j,bi,bj)
           ENDDO
          ENDDO
         ELSE
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            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=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
             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=ticekey, byte=isbyte
CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, 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,bi,bj) = 0. _d 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(:,:,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 snowPrc(:,:,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
#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 */

#ifdef ALLOW_BULK_FORCE
        IF ( useBulkforce ) THEN
         CALL THSICE_GET_PRECIP(
     I                  iceMask,
     O                  prcAtm(1-OLx,1-OLy,bi,bj),
     O                  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(1-OLx,1-OLy,bi,bj),
     O                  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_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.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
#endif
          CALL THSICE_STEP_FWD(
     I                     bi, bj, iMin, iMax, jMin, jMax,
     I                     prcAtm(1-OLx,1-OLy,bi,bj),
     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

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        CALL THSICE_AVE(
     I                     bi,bj, myTime, myIter, myThid )
       ENDDO
      ENDDO

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. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
        CALL THSICE_DO_EXCH( myThid )
      ENDIF
#ifdef OLD_THSICE_CALL_SEQUENCE
#ifdef ATMOSPHERIC_LOADING
      IF ( useRealFreshWaterFlux .AND.
     &    ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
     &  _EXCH_XY_RS( sIceLoad, myThid )
#endif
#else /* OLD_THSICE_CALL_SEQUENCE */
#ifdef ATMOSPHERIC_LOADING
      IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
     &                           .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--   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 = 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

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

      RETURN
      END
1	jmc	1.27	C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_main.F,v 1.26 2012/11/21 01:53:34 jmc 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	gforget	1.23	#include "THSICE_SIZE.h"
32	jmc	1.2	#include "THSICE_VARS.h"
33	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
34			# include "tamc.h"
35			# include "tamc_keys.h"
36			#endif
37	jmc	1.8
38	jmc	1.1	C !INPUT/OUTPUT PARAMETERS:
39			C === Routine arguments ===
40	jmc	1.12	C myTime :: Current time in simulation (s)
41			C myIter :: Current iteration number
42			C myThid :: My Thread Id. number
43			_RL myTime
44	jmc	1.1	INTEGER myIter
45			INTEGER myThid
46	jmc	1.2	CEOP
47	jmc	1.1
48			#ifdef ALLOW_THSICE
49			C !LOCAL VARIABLES:
50			C === Local variables ===
51			INTEGER i,j
52			INTEGER bi,bj
53			INTEGER iMin, iMax
54			INTEGER jMin, jMax
55	jmc	1.25	_RL prcAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
56	jmc	1.8	c _RL evpAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57			c _RL flxAtm(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58			c _RL flxSW (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59	jmc	1.1
60			_RL tauFac
61
62			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
63
64	jmc	1.24	IF ( useEXF .OR. useSEAICE ) THEN
65			C- EXF does not provide valid fields in overlap
66			iMin = 1
67			iMax = sNx
68			jMin = 1
69			jMax = sNy
70			ELSEIF ( stressReduction.GT. 0. _d 0 ) THEN
71	jmc	1.5	C- needs new Ice Fraction in halo region to apply wind-stress reduction
72	jmc	1.8	iMin = 1-OLx
73			iMax = sNx+OLx-1
74			jMin = 1-OLy
75			jMax = sNy+OLy-1
76	jmc	1.5	#ifdef ATMOSPHERIC_LOADING
77	jmc	1.24	ELSEIF ( useRealFreshWaterFlux ) THEN
78	jmc	1.5	C- needs sea-ice loading in part of the halo regions for grad.Phi0surf
79			C to be valid at the boundaries ( d/dx 1:sNx+1 ; d/dy 1:sNy+1 )
80			iMin = 0
81			iMax = sNx+1
82			jMin = 0
83			jMax = sNy+1
84			#endif /* ATMOSPHERIC_LOADING */
85	jmc	1.1	ELSE
86			iMin = 1
87			iMax = sNx
88			jMin = 1
89			jMax = sNy
90			ENDIF
91
92			DO bj=myByLo(myThid),myByHi(myThid)
93			DO bi=myBxLo(myThid),myBxHi(myThid)
94	jmc	1.2
95	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
96			act1 = bi - myBxLo(myThid)
97			max1 = myBxHi(myThid) - myBxLo(myThid) + 1
98			act2 = bj - myByLo(myThid)
99			max2 = myByHi(myThid) - myByLo(myThid) + 1
100			act3 = myThid - 1
101			max3 = nTx*nTy
102			act4 = ikey_dynamics - 1
103	gforget	1.23	ticekey = (act1 + 1) + act2*max1
104	heimbach	1.7	& + act3max1max2
105			& + act4max1max2*max3
106			#endif /* ALLOW_AUTODIFF_TAMC */
107
108	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
109	gforget	1.23	CADJ STORE ocefwfx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
110			CADJ STORE oceqnet(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
111			CADJ STORE ocesflx(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
112	heimbach	1.14	# ifdef ALLOW_EXF
113	gforget	1.23	CADJ STORE qsw(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
114	heimbach	1.14	# endif
115	heimbach	1.13	#endif
116
117	jmc	1.3	C-- Mixed layer thickness: take the 1rst layer
118			#ifdef NONLIN_FRSURF
119			IF ( staggerTimeStep .AND. nonlinFreeSurf.GT.0 ) THEN
120			IF ( select_rStar.GT.0 ) THEN
121	jmc	1.20	DO j=1-OLy,sNy+OLy
122			DO i=1-OLx,sNx+OLx
123	jmc	1.3	hOceMxL(i,j,bi,bj) = drF(1)*h0FacC(i,j,1,bi,bj)
124	jmc	1.8	& *rStarFacC(i,j,bi,bj)
125	jmc	1.3	ENDDO
126			ENDDO
127			ELSE
128	jmc	1.20	DO j=1-OLy,sNy+OLy
129			DO i=1-OLx,sNx+OLx
130	jmc	1.25	IF ( kSurfC(i,j,bi,bj).EQ.1 ) THEN
131	jmc	1.3	hOceMxL(i,j,bi,bj) = drF(1)*hFac_surfC(i,j,bi,bj)
132			ELSE
133	jmc	1.12	hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
134	jmc	1.3	ENDIF
135			ENDDO
136			ENDDO
137			ENDIF
138			ELSE
139			#else /* ndef NONLIN_FRSURF */
140			IF (.TRUE.) THEN
141			#endif /* NONLIN_FRSURF */
142	jmc	1.20	DO j=1-OLy,sNy+OLy
143			DO i=1-OLx,sNx+OLx
144	jmc	1.12	hOceMxL(i,j,bi,bj) = drF(1)*hFacC(i,j,1,bi,bj)
145	jmc	1.3	ENDDO
146			ENDDO
147			ENDIF
148
149	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
150	gforget	1.23	CADJ STORE uvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
151			CADJ STORE vvel (:,:,1,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
152	heimbach	1.7	#endif
153
154	mlosch	1.10	DO j = jMin, jMax
155			DO i = iMin, iMax
156			tOceMxL(i,j,bi,bj) = theta(i,j,1,bi,bj)
157			sOceMxL(i,j,bi,bj) = salt (i,j,1,bi,bj)
158			v2ocMxL(i,j,bi,bj) =
159	jmc	1.2	& ( uvel(i,j,1,bi,bj)*uvel(i,j,1,bi,bj)
160			& + uvel(i+1,j,1,bi,bj)*uvel(i+1,j,1,bi,bj)
161			& + vvel(i,j+1,1,bi,bj)*vvel(i,j+1,1,bi,bj)
162			& + vvel(i,j,1,bi,bj)*vvel(i,j,1,bi,bj)
163			& )*0.5 _d 0
164	jmc	1.25	prcAtm (i,j,bi,bj) = 0. _d 0
165	mlosch	1.10	icFrwAtm(i,j,bi,bj) = 0. _d 0
166			icFlxAtm(i,j,bi,bj) = 0. _d 0
167			icFlxSW (i,j,bi,bj) = 0. _d 0
168	jmc	1.25	snowPrc (i,j,bi,bj) = 0. _d 0
169			siceAlb (i,j,bi,bj) = 0. _d 0
170	jmc	1.6	ENDDO
171	mlosch	1.10	ENDDO
172	heimbach	1.7
173			#ifdef ALLOW_AUTODIFF_TAMC
174	gforget	1.23	CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
175			CADJ STORE iceHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
176			CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
177			CADJ STORE Tsrf(:,:,bi,bj) = comlev1_bibj, key = ticekey
178			CADJ STORE Qice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
179			CADJ STORE Qice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
180			CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
181			CADJ STORE snowPrc(:,:,bi,bj) = comlev1_bibj, key = ticekey
182
183			CADJ STORE hOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
184			CADJ STORE tOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
185			CADJ STORE sOceMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
186			CADJ STORE v2ocMxL(:,:,bi,bj) = comlev1_bibj, key = ticekey
187	heimbach	1.7	#endif
188
189	jmc	1.28	#ifdef OLD_THSICE_CALL_SEQUENCE
190	jmc	1.12	C- do sea-ice advection before getting surface fluxes
191			C Note: will inline this S/R once thSIce in Atmos. set-up is settled
192			IF ( thSIceAdvScheme.GT.0 )
193			& CALL THSICE_DO_ADVECT(
194			I bi,bj, myTime, myIter, myThid )
195	jmc	1.28	#endif /* OLD_THSICE_CALL_SEQUENCE */
196	jmc	1.12
197	jmc	1.2	#ifdef ALLOW_BULK_FORCE
198	mlosch	1.10	IF ( useBulkforce ) THEN
199			CALL THSICE_GET_PRECIP(
200	jmc	1.6	I iceMask,
201	jmc	1.25	O prcAtm(1-OLx,1-OLy,bi,bj),
202			O snowPrc(1-OLx,1-OLy,bi,bj),
203	jmc	1.8	O icFlxSW(1-OLx,1-OLy,bi,bj),
204	jmc	1.6	I iMin,iMax,jMin,jMax, bi,bj, myThid )
205	mlosch	1.10	ENDIF
206	jmc	1.2	#endif
207	mlosch	1.9	#ifdef ALLOW_EXF
208	mlosch	1.10	IF ( useEXF ) THEN
209			CALL THSICE_MAP_EXF(
210	mlosch	1.9	I iceMask,
211	jmc	1.25	O prcAtm(1-OLx,1-OLy,bi,bj),
212			O snowPrc(1-OLx,1-OLy,bi,bj),
213	mlosch	1.9	O icFlxSW(1-OLx,1-OLy,bi,bj),
214			I iMin,iMax,jMin,jMax, bi,bj, myThid )
215	mlosch	1.10	ENDIF
216	mlosch	1.9	#endif
217	jmc	1.2
218	heimbach	1.21	#ifdef ALLOW_AUTODIFF_TAMC
219	jmc	1.26	CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
220	gforget	1.23	CADJ STORE tice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
221			CADJ STORE tice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
222	jmc	1.26	#else
223			IF ( .NOT.thSIce_skipThermo ) THEN
224	heimbach	1.21	#endif
225	jmc	1.26	CALL THSICE_STEP_TEMP(
226	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
227			I myTime, myIter, myThid )
228
229	heimbach	1.13	#ifdef ALLOW_AUTODIFF_TAMC
230	gforget	1.23	CADJ STORE empmr(:,:,bi,bj) = comlev1_bibj, key = ticekey
231			CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = ticekey
232			CADJ STORE iceMask(:,:,bi,bj) = comlev1_bibj, key = ticekey
233			CADJ STORE iceHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
234			CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj, key = ticekey
235			cphCADJ STORE Tsrf(:,:,bi,bj) = comlev1_bibj, key = ticekey
236			CADJ STORE Qice1(:,:,bi,bj) = comlev1_bibj, key = ticekey
237			CADJ STORE Qice2(:,:,bi,bj) = comlev1_bibj, key = ticekey
238			CADJ STORE snowAge(:,:,bi,bj) = comlev1_bibj, key = ticekey
239	jmc	1.26	CADJ STORE sHeating(:,:,bi,bj) = comlev1_bibj, key = ticekey
240	heimbach	1.13	#endif
241	jmc	1.26	CALL THSICE_STEP_FWD(
242	jmc	1.8	I bi, bj, iMin, iMax, jMin, jMax,
243	jmc	1.25	I prcAtm(1-OLx,1-OLy,bi,bj),
244	jmc	1.1	I myTime, myIter, myThid )
245	jmc	1.26	#ifndef ALLOW_AUTODIFF_TAMC
246			ENDIF
247			#endif
248	jmc	1.2
249	jmc	1.25	C-- end bi,bj loop
250			ENDDO
251			ENDDO
252
253			#ifdef ALLOW_BALANCE_FLUXES
254			C-- Balance net Fresh-Water flux from Atm+Land
255			IF ( thSIceBalanceAtmFW.NE.0 ) THEN
256	jmc	1.26	CALL THSICE_BALANCE_FRW(
257	jmc	1.25	I iMin, iMax, jMin, jMax,
258			I prcAtm, myTime, myIter, myThid )
259			ENDIF
260			#endif
261
262			DO bj=myByLo(myThid),myByHi(myThid)
263			DO bi=myBxLo(myThid),myBxHi(myThid)
264	mlosch	1.10	CALL THSICE_AVE(
265	jmc	1.8	I bi,bj, myTime, myIter, myThid )
266	jmc	1.24	ENDDO
267			ENDDO
268	jmc	1.1
269	jmc	1.24	C add a small piece of code to check AddFluid implementation:
270			c#include "thsice_test_addfluid.h"
271
272			C-- Exchange fields that are advected by seaice dynamics
273	jmc	1.28	IF ( useSEAICE .OR. thSIceAdvScheme.GT.0
274			& .OR. ( useEXF .AND. stressReduction.GT.zeroRL ) ) THEN
275			CALL THSICE_DO_EXCH( myThid )
276	jmc	1.24	ENDIF
277	jmc	1.28	#ifdef OLD_THSICE_CALL_SEQUENCE
278	jmc	1.24	#ifdef ATMOSPHERIC_LOADING
279	jmc	1.27	IF ( useRealFreshWaterFlux .AND.
280			& ( useEXF .OR. useSEAICE .OR. thSIceAdvScheme.GT.0 ) )
281	jmc	1.24	& _EXCH_XY_RS( sIceLoad, myThid )
282			#endif
283	jmc	1.28	#else /* OLD_THSICE_CALL_SEQUENCE */
284			#ifdef ATMOSPHERIC_LOADING
285			IF ( useRealFreshWaterFlux .AND. (useEXF.OR.useSEAICE )
286			& .AND. thSIceAdvScheme.LE.0 )
287			& _EXCH_XY_RS( sIceLoad, myThid )
288			#endif
289
290			C- when useSEAICE=.true., this S/R is called from SEAICE_MODEL;
291			C otherwise, call it from here, after thsice-thermodynamics is done
292			IF ( thSIceAdvScheme.GT.0 .AND. .NOT.useSEAICE ) THEN
293			CALL THSICE_DO_ADVECT(
294			I 0, 0, myTime, myIter, myThid )
295			ENDIF
296			#endif /* OLD_THSICE_CALL_SEQUENCE */
297	jmc	1.24
298			DO bj=myByLo(myThid),myByHi(myThid)
299			DO bi=myBxLo(myThid),myBxHi(myThid)
300	jmc	1.11	C-- note: If useSEAICE=.true., the stress is computed in seaice_model,
301			C-- and stressReduction is always set to zero
302	heimbach	1.7	#ifdef ALLOW_AUTODIFF_TAMC
303	gforget	1.23	CADJ STORE fu(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
304			CADJ STORE fv(:,:,bi,bj) = comlev1_bibj, key=ticekey, byte=isbyte
305	heimbach	1.7	#endif
306	jmc	1.11	IF ( stressReduction.GT. 0. _d 0 ) THEN
307	mlosch	1.10	DO j = jMin, jMax
308			DO i = iMin+1,iMax
309	jmc	1.1	tauFac = stressReduction
310			& (iceMask(i-1,j,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
311			fu(i,j,bi,bj) = (1. _d 0 - tauFac)*fu(i,j,bi,bj)
312	mlosch	1.10	ENDDO
313	jmc	1.1	ENDDO
314	mlosch	1.10	DO j = jMin+1, jMax
315			DO i = iMin, iMax
316	jmc	1.1	tauFac = stressReduction
317			& (iceMask(i,j-1,bi,bj)+iceMask(i,j,bi,bj))0.5 _d 0
318			fv(i,j,bi,bj) = (1. _d 0 - tauFac)*fv(i,j,bi,bj)
319	mlosch	1.10	ENDDO
320	jmc	1.1	ENDDO
321			ENDIF
322
323			C-- end bi,bj loop
324			ENDDO
325			ENDDO
326
327			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
328			#endif /ALLOW_THSICE/
329
330			RETURN
331			END