pkg/seaice/seaice_advdiff.F

C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advdiff.F,v 1.17 2007/04/30 22:48:52 mlosch Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"

CBOP
C !ROUTINE: SEAICE_ADVDIFF

C !INTERFACE: ==========================================================
      SUBROUTINE SEAICE_ADVDIFF(
     I     myTime, myIter, myThid )

C !DESCRIPTION: \bv
C     *===========================================================*
C     | SUBROUTINE SEAICE_ADVDIFF
C     | o driver for different advection routines
C     |   calls an adaption of gad_advection to call different
C     |   advection routines of pkg/generic_advdiff
C     *===========================================================*
C \ev

C !USES: ===============================================================
      IMPLICIT NONE

C     === Global variables ===
C     UICE/VICE :: ice velocity
C     HEFF      :: scalar field to be advected
C     HEFFM     :: mask for scalar field
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "GAD.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"

#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif

C !INPUT PARAMETERS: ===================================================
C     === Routine arguments ===
C     myTime    :: current time
C     myIter    :: iteration number
C     myThid    :: Thread no. that called this routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid
CEndOfInterface

#ifdef ALLOW_SEAICE
C !LOCAL VARIABLES: ====================================================
C     === Local variables ===
C     i,j,bi,bj :: Loop counters
C     ks        :: surface level index
C     uc/vc     :: current ice velocity on C-grid
C     uTrans    :: volume transport, x direction
C     vTrans    :: volume transport, y direction
C     iceFld    :: copy of seaice field
C     afx       :: horizontal advective flux, x direction
C     afy       :: horizontal advective flux, y direction
C     gFld      :: tendency of seaice field
C     xA,yA     :: "areas" of X and Y face of tracer cells
C     msgBuf    :: Informational/error meesage buffer
      INTEGER i, j, bi, bj
      INTEGER ks
      LOGICAL SEAICEmultiDimAdvection
      CHARACTER*(MAX_LEN_MBUF) msgBuf

      _RL uc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL uTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL iceFld    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afx       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afy       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL gFld      (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL fld3d     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy)
      _RS xA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

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

      ks = 1

C--   Get rid of the time dimension for velocities and interpolate
C--   to C-points if necessary
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifdef SEAICE_CGRID
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          uc(i,j,bi,bj)=UICE(i,j,1,bi,bj)
          vc(i,j,bi,bj)=VICE(i,j,1,bi,bj)
         ENDDO
        ENDDO
#else /* not SEAICE_CGRID = BGRID */
C average seaice velocity to c-grid
        DO j=1-Oly,sNy+Oly-1
         DO i=1-Olx,sNx+Olx-1
          uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,1,bi,bj)+UICE(i,j+1,1,bi,bj))
          vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,1,bi,bj)+VICE(I+1,J,1,bi,bj))
         ENDDO
        ENDDO
#endif /* SEAICE_CGRID */
       ENDDO
      ENDDO

#ifndef SEAICE_CGRID
C     Do we need this? I am afraid so.
      CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid)
#endif /* not SEAICE_CGRID */

      SEAICEmultidimadvection = .TRUE.
      IF ( SEAICEadvScheme.EQ.ENUM_CENTERED_2ND
     & .OR.SEAICEadvScheme.EQ.ENUM_UPWIND_3RD
     & .OR.SEAICEadvScheme.EQ.ENUM_CENTERED_4TH ) THEN
       SEAICEmultiDimAdvection = .FALSE.
      ENDIF


      IF ( SEAICEmultiDimAdvection ) THEN
#ifndef ALLOW_AUTODIFF_TAMC
C     This has to be done to comply with the time stepping in advect.F:
C     Making sure that the following routines see the different
C     time levels correctly
C     At the end of the routine ADVECT,
C     timelevel 1 is updated with advection contribution
C                 and diffusion contribution
C                 (which was computed in DIFFUS on timelevel 3)
C     timelevel 2 is the previous timelevel 1
C     timelevel 3 is the total diffusion tendency * deltaT
C                 (empty if no diffusion)

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uc     = comlev1, key = ikey_dynamics
CADJ STORE vc     = comlev1, key = ikey_dynamics
#endif /* ALLOW_AUTODIFF_TAMC */

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
C---   loops on tile indices bi,bj

#ifdef ALLOW_AUTODIFF_TAMC
C     Initialise for TAF
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          iceFld(i,j)     = 0. _d 0
          gFld(i,j)       = 0. _d 0
         ENDDO
        ENDDO
#endif /* ALLOW_AUTODIFF_TAMC */

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          HEFF(i,j,3,bi,bj) = 0. _d 0
          HEFF(i,j,2,bi,bj) = HEFF(i,j,1,bi,bj)
          AREA(i,j,3,bi,bj) = 0. _d 0
          AREA(i,j,2,bi,bj) = AREA(i,j,1,bi,bj)
          recip_heff(i,j)   = 1. _d 0
         ENDDO
        ENDDO

C-    first compute cell areas used by all tracers
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          xA(i,j) = _dyG(i,j,bi,bj)*_maskW(i,j,ks,bi,bj)
          yA(i,j) = _dxG(i,j,bi,bj)*_maskS(i,j,ks,bi,bj)
         ENDDO
        ENDDO
C-    Calculate "volume transports" through tracer cell faces.
        DO j=1-Oly,sNy+Oly
         DO i=1-Olx,sNx+Olx
          uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j)
          vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j)
         ENDDO
        ENDDO

C--   Effective Thickness (Volume)
        IF ( SEAICEadvHeff ) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           iceFld(i,j) = HEFF(i,j,1,bi,bj)
          ENDDO
         ENDDO
         CALL SEAICE_ADVECTION(
     I        GAD_HEFF, SEAICEadvSchHeff,
     I        uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I        uTrans, vTrans, iceFld, recip_heff,
     O        gFld, afx, afy,
     I        bi, bj, myTime, myIter, myThid )
         IF ( diff1 .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_HEFF,
     I         HEFF(1-OLx,1-OLy,1,bi,bj), HEFFM, xA, yA,
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           HEFF(i,j,1,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          HEFF(i,j,1,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
        ENDIF

C--   Fractional area
        IF ( SEAICEadvArea ) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           iceFld(i,j) = AREA(i,j,1,bi,bj)
          ENDDO
         ENDDO
         CALL SEAICE_ADVECTION(
     I        GAD_AREA, SEAICEadvSchArea,
     I        uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I        uTrans, vTrans, iceFld, recip_heff,
     O        gFld, afx, afy,
     I        bi, bj, myTime, myIter, myThid )
         IF ( diff1 .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_AREA,
     I         AREA(1-OLx,1-OLy,1,bi,bj), HEFFM, xA, yA,
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           AREA(i,j,1,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          AREA(i,j,1,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
        ENDIF

C--   Effective Snow Thickness (Volume)
        IF ( SEAICEadvSnow ) THEN
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           iceFld(i,j) = HSNOW(i,j,bi,bj)
          ENDDO
         ENDDO
         CALL SEAICE_ADVECTION(
     I        GAD_SNOW, SEAICEadvSchSnow,
     I        uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I        uTrans, vTrans, iceFld, recip_heff,
     O        gFld, afx, afy,
     I        bi, bj, myTime, myIter, myThid )
         IF ( diff1 .GT. 0. _d 0 ) THEN
C--   Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_SNOW,
     I         HSNOW(1-OLx,1-OLy,bi,bj), HEFFM, xA, yA,
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
        ENDIF
C---   end bi,bj loops
       ENDDO
      ENDDO

#else
      STOP 'SEAICEmultiDimAdvection not yet implemented for adjoint'
#endif /* ndef ALLOW_AUTODIFF_TAMC */
      ELSE
C--   if not multiDimAdvection

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uc   = comlev1, key = ikey_dynamics
CADJ STORE vc   = comlev1, key = ikey_dynamics
#endif /* ALLOW_AUTODIFF_TAMC */

       IF ( SEAICEadvHEff ) CALL ADVECT( uc, vc, HEFF, HEFFM, myThid )
       IF ( SEAICEadvArea ) CALL ADVECT( uc, vc, AREA, HEFFM, myThid )
       IF ( SEAICEadvSnow ) THEN
C     another fudge: copy 2D field HSNOW to 3D field fld3d to be able to
C     use ADVECT. Works only with LAD=2 (Backward Euler scheme) and NOT
C     with LAD=1 (Leapfrog).
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            fld3d(I,J,3,bi,bj) = 0. _d 0 
            fld3d(I,J,2,bi,bj) = HSNOW(I,J,bi,bj)
            fld3d(I,J,1,bi,bj) = HSNOW(I,J,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        CALL ADVECT( uc, vc, fld3d, HEFFM, myThid )       
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            HSNOW(I,J,bi,bj) = fld3d(I,J,1,bi,bj) 
           ENDDO
          ENDDO
         ENDDO
        ENDDO
       ENDIF

C--   end if multiDimAdvection
      ENDIF
C--   end if SEAICEuseDynamics

      IF ( .NOT. usePW79thermodynamics ) THEN
C     Hiblers "ridging function": Do it now if not in seaice_growth
C     in principle we should add a "real" ridging function here (or 
C     somewhere after doing the advection)
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           AREA(I,J,1,bi,bj) = MIN(ONE, AREA(I,J,1,bi,bj))
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDIF
#endif /* ALLOW_SEAICE */

      RETURN
      END
1	mlosch	1.18	C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advdiff.F,v 1.17 2007/04/30 22:48:52 mlosch Exp $
2	mlosch	1.1	C $Name: $
3
4			#include "SEAICE_OPTIONS.h"
5	jmc	1.11
6	mlosch	1.1	CBOP
7			C !ROUTINE: SEAICE_ADVDIFF
8
9			C !INTERFACE: ==========================================================
10	jmc	1.11	SUBROUTINE SEAICE_ADVDIFF(
11	mlosch	1.1	I myTime, myIter, myThid )
12
13			C !DESCRIPTION: \bv
14	jmc	1.11	C ===========================================================
15			C \| SUBROUTINE SEAICE_ADVDIFF
16			C \| o driver for different advection routines
17			C \| calls an adaption of gad_advection to call different
18			C \| advection routines of pkg/generic_advdiff
19			C ===========================================================
20			C \ev
21
22			C !USES: ===============================================================
23	mlosch	1.1	IMPLICIT NONE
24	jmc	1.11
25	mlosch	1.1	C === Global variables ===
26	jmc	1.11	C UICE/VICE :: ice velocity
27			C HEFF :: scalar field to be advected
28			C HEFFM :: mask for scalar field
29	mlosch	1.1	#include "SIZE.h"
30			#include "EEPARAMS.h"
31			#include "PARAMS.h"
32			#include "GRID.h"
33			#include "GAD.h"
34			#include "SEAICE_PARAMS.h"
35			#include "SEAICE.h"
36
37			#ifdef ALLOW_AUTODIFF_TAMC
38			# include "tamc.h"
39			#endif
40
41			C !INPUT PARAMETERS: ===================================================
42			C === Routine arguments ===
43	jmc	1.11	C myTime :: current time
44			C myIter :: iteration number
45			C myThid :: Thread no. that called this routine.
46	mlosch	1.1	_RL myTime
47			INTEGER myIter
48			INTEGER myThid
49			CEndOfInterface
50
51			#ifdef ALLOW_SEAICE
52			C !LOCAL VARIABLES: ====================================================
53			C === Local variables ===
54	jmc	1.11	C i,j,bi,bj :: Loop counters
55			C ks :: surface level index
56			C uc/vc :: current ice velocity on C-grid
57			C uTrans :: volume transport, x direction
58			C vTrans :: volume transport, y direction
59			C iceFld :: copy of seaice field
60			C afx :: horizontal advective flux, x direction
61			C afy :: horizontal advective flux, y direction
62			C gFld :: tendency of seaice field
63			C xA,yA :: "areas" of X and Y face of tracer cells
64			C msgBuf :: Informational/error meesage buffer
65			INTEGER i, j, bi, bj
66			INTEGER ks
67	mlosch	1.1	LOGICAL SEAICEmultiDimAdvection
68	mlosch	1.6	CHARACTER*(MAX_LEN_MBUF) msgBuf
69	mlosch	1.1
70	mlosch	1.6	_RL uc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
71			_RL vc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
72	jmc	1.11	_RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
73			_RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
74			_RL iceFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
75			_RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
76			_RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
77			_RL gFld (1-Olx:sNx+Olx,1-Oly:sNy+Oly)
78	mlosch	1.12	_RL fld3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,3,nSx,nSy)
79	mlosch	1.6	_RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
80			_RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
81	mlosch	1.9	_RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
82	jmc	1.8	CEOP
83	jmc	1.11
84			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
85
86			ks = 1
87
88	mlosch	1.6	C-- Get rid of the time dimension for velocities and interpolate
89			C-- to C-points if necessary
90			DO bj=myByLo(myThid),myByHi(myThid)
91			DO bi=myBxLo(myThid),myBxHi(myThid)
92			#ifdef SEAICE_CGRID
93			DO j=1-Oly,sNy+Oly
94			DO i=1-Olx,sNx+Olx
95	jmc	1.11	uc(i,j,bi,bj)=UICE(i,j,1,bi,bj)
96			vc(i,j,bi,bj)=VICE(i,j,1,bi,bj)
97	mlosch	1.6	ENDDO
98			ENDDO
99			#else /* not SEAICE_CGRID = BGRID */
100			C average seaice velocity to c-grid
101			DO j=1-Oly,sNy+Oly-1
102			DO i=1-Olx,sNx+Olx-1
103	jmc	1.11	uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,1,bi,bj)+UICE(i,j+1,1,bi,bj))
104			vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,1,bi,bj)+VICE(I+1,J,1,bi,bj))
105	mlosch	1.6	ENDDO
106			ENDDO
107			#endif /* SEAICE_CGRID */
108			ENDDO
109			ENDDO
110
111			#ifndef SEAICE_CGRID
112			C Do we need this? I am afraid so.
113			CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid)
114			#endif /* not SEAICE_CGRID */
115
116	mlosch	1.1	SEAICEmultidimadvection = .TRUE.
117			IF ( SEAICEadvScheme.EQ.ENUM_CENTERED_2ND
118			& .OR.SEAICEadvScheme.EQ.ENUM_UPWIND_3RD
119			& .OR.SEAICEadvScheme.EQ.ENUM_CENTERED_4TH ) THEN
120			SEAICEmultiDimAdvection = .FALSE.
121			ENDIF
122
123
124			IF ( SEAICEmultiDimAdvection ) THEN
125	heimbach	1.7	#ifndef ALLOW_AUTODIFF_TAMC
126	mlosch	1.1	C This has to be done to comply with the time stepping in advect.F:
127	jmc	1.11	C Making sure that the following routines see the different
128	mlosch	1.1	C time levels correctly
129	jmc	1.11	C At the end of the routine ADVECT,
130			C timelevel 1 is updated with advection contribution
131			C and diffusion contribution
132	mlosch	1.1	C (which was computed in DIFFUS on timelevel 3)
133			C timelevel 2 is the previous timelevel 1
134			C timelevel 3 is the total diffusion tendency * deltaT
135			C (empty if no diffusion)
136
137	heimbach	1.5	#ifdef ALLOW_AUTODIFF_TAMC
138	mlosch	1.6	CADJ STORE uc = comlev1, key = ikey_dynamics
139			CADJ STORE vc = comlev1, key = ikey_dynamics
140	heimbach	1.5	#endif /* ALLOW_AUTODIFF_TAMC */
141
142	mlosch	1.1	DO bj=myByLo(myThid),myByHi(myThid)
143			DO bi=myBxLo(myThid),myBxHi(myThid)
144	jmc	1.11	C--- loops on tile indices bi,bj
145
146			#ifdef ALLOW_AUTODIFF_TAMC
147			C Initialise for TAF
148			DO j=1-Oly,sNy+Oly
149			DO i=1-Olx,sNx+Olx
150			iceFld(i,j) = 0. _d 0
151			gFld(i,j) = 0. _d 0
152			ENDDO
153			ENDDO
154			#endif /* ALLOW_AUTODIFF_TAMC */
155
156	mlosch	1.1	DO j=1-OLy,sNy+OLy
157			DO i=1-OLx,sNx+OLx
158	jmc	1.11	HEFF(i,j,3,bi,bj) = 0. _d 0
159			HEFF(i,j,2,bi,bj) = HEFF(i,j,1,bi,bj)
160			AREA(i,j,3,bi,bj) = 0. _d 0
161			AREA(i,j,2,bi,bj) = AREA(i,j,1,bi,bj)
162			recip_heff(i,j) = 1. _d 0
163	mlosch	1.1	ENDDO
164			ENDDO
165
166	jmc	1.11	C- first compute cell areas used by all tracers
167	mlosch	1.2	DO j=1-Oly,sNy+Oly
168			DO i=1-Olx,sNx+Olx
169	jmc	1.11	xA(i,j) = _dyG(i,j,bi,bj)*_maskW(i,j,ks,bi,bj)
170			yA(i,j) = _dxG(i,j,bi,bj)*_maskS(i,j,ks,bi,bj)
171	mlosch	1.2	ENDDO
172			ENDDO
173	jmc	1.11	C- Calculate "volume transports" through tracer cell faces.
174	mlosch	1.1	DO j=1-Oly,sNy+Oly
175			DO i=1-Olx,sNx+Olx
176	jmc	1.11	uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j)
177			vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j)
178	mlosch	1.1	ENDDO
179			ENDDO
180	jmc	1.11
181	jmc	1.8	C-- Effective Thickness (Volume)
182	mlosch	1.17	IF ( SEAICEadvHeff ) THEN
183			DO j=1-Oly,sNy+Oly
184			DO i=1-Olx,sNx+Olx
185			iceFld(i,j) = HEFF(i,j,1,bi,bj)
186			ENDDO
187	mlosch	1.1	ENDDO
188	mlosch	1.17	CALL SEAICE_ADVECTION(
189			I GAD_HEFF, SEAICEadvSchHeff,
190			I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
191			I uTrans, vTrans, iceFld, recip_heff,
192			O gFld, afx, afy,
193			I bi, bj, myTime, myIter, myThid )
194			IF ( diff1 .GT. 0. _d 0 ) THEN
195	jmc	1.11	C- Add tendency due to diffusion
196			CALL SEAICE_DIFFUSION(
197	mlosch	1.17	I GAD_HEFF,
198			I HEFF(1-OLx,1-OLy,1,bi,bj), HEFFM, xA, yA,
199			U gFld,
200			I bi, bj, myTime, myIter, myThid )
201			ENDIF
202	mlosch	1.1	C now do the "explicit" time step
203	mlosch	1.17	DO j=1,sNy
204			DO i=1,sNx
205			HEFF(i,j,1,bi,bj) = HEFFM(i,j,bi,bj) * (
206			& HEFF(i,j,1,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
207			& )
208			ENDDO
209	mlosch	1.1	ENDDO
210	mlosch	1.17	ENDIF
211	jmc	1.11
212			C-- Fractional area
213	mlosch	1.17	IF ( SEAICEadvArea ) THEN
214			DO j=1-Oly,sNy+Oly
215			DO i=1-Olx,sNx+Olx
216			iceFld(i,j) = AREA(i,j,1,bi,bj)
217			ENDDO
218	jmc	1.11	ENDDO
219	mlosch	1.17	CALL SEAICE_ADVECTION(
220			I GAD_AREA, SEAICEadvSchArea,
221			I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
222			I uTrans, vTrans, iceFld, recip_heff,
223			O gFld, afx, afy,
224			I bi, bj, myTime, myIter, myThid )
225			IF ( diff1 .GT. 0. _d 0 ) THEN
226	jmc	1.11	C- Add tendency due to diffusion
227			CALL SEAICE_DIFFUSION(
228	mlosch	1.17	I GAD_AREA,
229			I AREA(1-OLx,1-OLy,1,bi,bj), HEFFM, xA, yA,
230			U gFld,
231			I bi, bj, myTime, myIter, myThid )
232			ENDIF
233	jmc	1.11	C now do the "explicit" time step
234	mlosch	1.17	DO j=1,sNy
235			DO i=1,sNx
236			AREA(i,j,1,bi,bj) = HEFFM(i,j,bi,bj) * (
237			& AREA(i,j,1,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
238			& )
239			ENDDO
240	jmc	1.11	ENDDO
241	mlosch	1.17	ENDIF
242	jmc	1.11
243	mlosch	1.12	C-- Effective Snow Thickness (Volume)
244	mlosch	1.17	IF ( SEAICEadvSnow ) THEN
245			DO j=1-Oly,sNy+Oly
246			DO i=1-Olx,sNx+Olx
247			iceFld(i,j) = HSNOW(i,j,bi,bj)
248			ENDDO
249	mlosch	1.12	ENDDO
250	mlosch	1.17	CALL SEAICE_ADVECTION(
251			I GAD_SNOW, SEAICEadvSchSnow,
252			I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
253			I uTrans, vTrans, iceFld, recip_heff,
254			O gFld, afx, afy,
255			I bi, bj, myTime, myIter, myThid )
256			IF ( diff1 .GT. 0. _d 0 ) THEN
257			C-- Add tendency due to diffusion
258	mlosch	1.12	CALL SEAICE_DIFFUSION(
259	mlosch	1.17	I GAD_SNOW,
260			I HSNOW(1-OLx,1-OLy,bi,bj), HEFFM, xA, yA,
261			U gFld,
262			I bi, bj, myTime, myIter, myThid )
263			ENDIF
264	mlosch	1.12	C now do the "explicit" time step
265	mlosch	1.17	DO j=1,sNy
266			DO i=1,sNx
267			HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
268			& HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
269			& )
270			ENDDO
271	mlosch	1.12	ENDDO
272			ENDIF
273	jmc	1.11	C--- end bi,bj loops
274	mlosch	1.1	ENDDO
275			ENDDO
276
277	heimbach	1.7	#else
278			STOP 'SEAICEmultiDimAdvection not yet implemented for adjoint'
279			#endif /* ndef ALLOW_AUTODIFF_TAMC */
280	mlosch	1.1	ELSE
281	mlosch	1.6	C-- if not multiDimAdvection
282	heimbach	1.5
283			#ifdef ALLOW_AUTODIFF_TAMC
284	mlosch	1.6	CADJ STORE uc = comlev1, key = ikey_dynamics
285			CADJ STORE vc = comlev1, key = ikey_dynamics
286	heimbach	1.5	#endif /* ALLOW_AUTODIFF_TAMC */
287
288	mlosch	1.17	IF ( SEAICEadvHEff ) CALL ADVECT( uc, vc, HEFF, HEFFM, myThid )
289			IF ( SEAICEadvArea ) CALL ADVECT( uc, vc, AREA, HEFFM, myThid )
290			IF ( SEAICEadvSnow ) THEN
291	mlosch	1.12	C another fudge: copy 2D field HSNOW to 3D field fld3d to be able to
292			C use ADVECT. Works only with LAD=2 (Backward Euler scheme) and NOT
293			C with LAD=1 (Leapfrog).
294			DO bj=myByLo(myThid),myByHi(myThid)
295			DO bi=myBxLo(myThid),myBxHi(myThid)
296			DO j=1-OLy,sNy+OLy
297			DO i=1-OLx,sNx+OLx
298			fld3d(I,J,3,bi,bj) = 0. _d 0
299			fld3d(I,J,2,bi,bj) = HSNOW(I,J,bi,bj)
300			fld3d(I,J,1,bi,bj) = HSNOW(I,J,bi,bj)
301			ENDDO
302			ENDDO
303			ENDDO
304			ENDDO
305			CALL ADVECT( uc, vc, fld3d, HEFFM, myThid )
306			DO bj=myByLo(myThid),myByHi(myThid)
307			DO bi=myBxLo(myThid),myBxHi(myThid)
308			DO j=1-OLy,sNy+OLy
309			DO i=1-OLx,sNx+OLx
310			HSNOW(I,J,bi,bj) = fld3d(I,J,1,bi,bj)
311			ENDDO
312			ENDDO
313			ENDDO
314			ENDDO
315			ENDIF
316	heimbach	1.5
317	mlosch	1.6	C-- end if multiDimAdvection
318	mlosch	1.1	ENDIF
319	mlosch	1.6	C-- end if SEAICEuseDynamics
320	mlosch	1.18
321			IF ( .NOT. usePW79thermodynamics ) THEN
322			C Hiblers "ridging function": Do it now if not in seaice_growth
323			C in principle we should add a "real" ridging function here (or
324			C somewhere after doing the advection)
325			DO bj=myByLo(myThid),myByHi(myThid)
326			DO bi=myBxLo(myThid),myBxHi(myThid)
327			DO j=1-Oly,sNy+Oly
328			DO i=1-Olx,sNx+Olx
329			AREA(I,J,1,bi,bj) = MIN(ONE, AREA(I,J,1,bi,bj))
330			ENDDO
331			ENDDO
332			ENDDO
333			ENDDO
334			ENDIF
335	mlosch	1.1	#endif /* ALLOW_SEAICE */
336
337			RETURN
338			END