/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_advdiff.F
ViewVC logotype

Annotation of /MITgcm_contrib/torge/itd/code/seaice_advdiff.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.5 - (hide annotations) (download)
Fri May 3 18:58:30 2013 UTC (12 years, 3 months ago) by torge
Branch: MAIN
CVS Tags: HEAD
Changes since 1.4: +21 -49 lines 
- changing counter variable looping through ITD space from "k" to "it"
  ("it" is used in seaice_growth)
- removing all "ToM" comments
1 torge 1.5 C $Header: /u/gcmpack/MITgcm_contrib/torge/itd/code/seaice_advdiff.F,v 1.4 2013/03/27 18:59:52 torge Exp $
2 dimitri 1.1 C $Name: $
3    
4     #include "SEAICE_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: SEAICE_ADVDIFF
8    
9     C !INTERFACE: ==========================================================
10     SUBROUTINE SEAICE_ADVDIFF(
11     I myTime, myIter, myThid )
12    
13     C !DESCRIPTION: \bv
14     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     IMPLICIT NONE
24    
25     C === Global variables ===
26     C UICE/VICE :: ice velocity
27     C HEFF :: scalar field to be advected
28     C HEFFM :: mask for scalar field
29     #include "SIZE.h"
30     #include "EEPARAMS.h"
31     #include "PARAMS.h"
32     #include "GRID.h"
33     #include "GAD.h"
34     #include "SEAICE_SIZE.h"
35     #include "SEAICE_PARAMS.h"
36     #include "SEAICE.h"
37     #include "SEAICE_TRACER.h"
38    
39     #ifdef ALLOW_AUTODIFF_TAMC
40     # include "tamc.h"
41     #endif
42    
43     C !INPUT PARAMETERS: ===================================================
44     C === Routine arguments ===
45     C myTime :: current time
46     C myIter :: iteration number
47     C myThid :: Thread no. that called this routine.
48     _RL myTime
49     INTEGER myIter
50     INTEGER myThid
51     CEndOfInterface
52    
53     C !LOCAL VARIABLES: ====================================================
54     C === Local variables ===
55     C i,j,bi,bj :: Loop counters
56 dimitri 1.2 #ifdef SEAICE_ITD
57 torge 1.5 C it :: Loop counter for ice thickness categories
58 dimitri 1.2 #endif
59 dimitri 1.1 C ks :: surface level index
60     C uc/vc :: current ice velocity on C-grid
61     C uTrans :: volume transport, x direction
62     C vTrans :: volume transport, y direction
63     C afx :: horizontal advective flux, x direction
64     C afy :: horizontal advective flux, y direction
65     C gFld :: tendency of seaice field
66     C xA,yA :: "areas" of X and Y face of tracer cells
67     INTEGER i, j, bi, bj
68 dimitri 1.2 #ifdef SEAICE_ITD
69 torge 1.5 INTEGER it
70 dimitri 1.2 #endif
71 dimitri 1.1 INTEGER ks
72     LOGICAL SEAICEmultiDimAdvection
73     #ifdef ALLOW_AUTODIFF_TAMC
74     INTEGER itmpkey
75     #endif /* ALLOW_AUTODIFF_TAMC */
76     #ifdef ALLOW_SITRACER
77     _RL hEffNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78     _RL areaNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79     INTEGER iTr, SEAICEadvSchSItr
80     _RL SEAICEdiffKhSItr
81     _RL SItrExt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82     _RL tmpscal1, tmpscal2
83     # ifdef ALLOW_SITRACER_ADVCAP
84     _RL SItrPrev (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
85     # endif
86     # ifdef ALLOW_SITRACER_DEBUG_DIAG
87     _RL DIAGarray (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
88     # endif
89     #endif /* ALLOW_SITRACER */
90    
91     _RL uc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
92     _RL vc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
93     _RL fldNm1 (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
94     _RL uTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
95     _RL vTrans (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
96     _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
97     _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
98     _RL gFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
99     _RS xA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
100     _RS yA (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
101     _RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
102     CEOP
103    
104     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
105    
106     ks = 1
107    
108     C-- make a local copy of the velocities for compatibility with B-grid
109     C-- alternatively interpolate to C-points if necessary
110     DO bj=myByLo(myThid),myByHi(myThid)
111     DO bi=myBxLo(myThid),myBxHi(myThid)
112     #ifdef SEAICE_CGRID
113     DO j=1-OLy,sNy+OLy
114     DO i=1-OLx,sNx+OLx
115     uc(i,j,bi,bj)=UICE(i,j,bi,bj)
116     vc(i,j,bi,bj)=VICE(i,j,bi,bj)
117     ENDDO
118     ENDDO
119     #else /* not SEAICE_CGRID = BGRID */
120     C average seaice velocity to C-grid
121     DO j=1-OLy,sNy+OLy-1
122     DO i=1-OLx,sNx+OLx-1
123     uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj))
124     vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj))
125     ENDDO
126     ENDDO
127     #endif /* SEAICE_CGRID */
128     C- compute cell areas used by all tracers
129     DO j=1-OLy,sNy+OLy
130     DO i=1-OLx,sNx+OLx
131     xA(i,j,bi,bj) = _dyG(i,j,bi,bj)*_maskW(i,j,ks,bi,bj)
132     yA(i,j,bi,bj) = _dxG(i,j,bi,bj)*_maskS(i,j,ks,bi,bj)
133     ENDDO
134     ENDDO
135     ENDDO
136     ENDDO
137    
138     #ifndef SEAICE_CGRID
139     C Do we need this? I am afraid so.
140     CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid)
141     #endif /* not SEAICE_CGRID */
142    
143     #ifdef ALLOW_AUTODIFF_TAMC
144     CADJ STORE uc = comlev1, key = ikey_dynamics, kind=isbyte
145     CADJ STORE vc = comlev1, key = ikey_dynamics, kind=isbyte
146     #endif /* ALLOW_AUTODIFF_TAMC */
147    
148     SEAICEmultidimadvection = .TRUE.
149     IF ( SEAICEadvScheme.EQ.ENUM_CENTERED_2ND
150     & .OR.SEAICEadvScheme.EQ.ENUM_UPWIND_3RD
151     & .OR.SEAICEadvScheme.EQ.ENUM_CENTERED_4TH ) THEN
152     SEAICEmultiDimAdvection = .FALSE.
153     ENDIF
154    
155     #ifdef ALLOW_AUTODIFF_TAMC
156     CADJ STORE heffm = comlev1, key = ikey_dynamics, kind=isbyte
157     #endif /* ALLOW_AUTODIFF_TAMC */
158     IF ( SEAICEmultiDimAdvection ) THEN
159    
160     DO bj=myByLo(myThid),myByHi(myThid)
161     DO bi=myBxLo(myThid),myBxHi(myThid)
162     C--- loops on tile indices bi,bj
163    
164     #ifdef ALLOW_AUTODIFF_TAMC
165     C Initialise for TAF
166     DO j=1-OLy,sNy+OLy
167     DO i=1-OLx,sNx+OLx
168     gFld(i,j) = 0. _d 0
169     ENDDO
170     ENDDO
171     C
172     act1 = bi - myBxLo(myThid)
173     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
174     act2 = bj - myByLo(myThid)
175     max2 = myByHi(myThid) - myByLo(myThid) + 1
176     act3 = myThid - 1
177     max3 = nTx*nTy
178     act4 = ikey_dynamics - 1
179     itmpkey = (act1 + 1) + act2*max1
180     & + act3*max1*max2
181     & + act4*max1*max2*max3
182     C
183     CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
184     CADJ & key = itmpkey, kind=isbyte
185     CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
186     CADJ & key = itmpkey, kind=isbyte
187     CADJ STORE heffm(:,:,bi,bj) = comlev1_bibj,
188     CADJ & key = itmpkey, kind=isbyte
189     CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
190     CADJ & key = itmpkey, kind=isbyte
191 torge 1.3 # ifdef SEAICE_VARIABLE_SALINITY
192     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
193     CADJ & key = itmpkey, kind=isbyte
194     # endif
195 dimitri 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
196    
197     DO j=1-OLy,sNy+OLy
198     DO i=1-OLx,sNx+OLx
199 torge 1.4 #ifdef ALLOW_SITRACER
200 dimitri 1.1 hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj)
201     areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj)
202     #endif
203     recip_heff(i,j) = 1. _d 0
204     ENDDO
205     ENDDO
206    
207     C- Calculate "volume transports" through tracer cell faces.
208     DO j=1-OLy,sNy+OLy
209     DO i=1-OLx,sNx+OLx
210     uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j,bi,bj)
211     vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j,bi,bj)
212     ENDDO
213     ENDDO
214    
215     C-- Effective Thickness (Volume)
216     IF ( SEAICEadvHeff ) THEN
217 dimitri 1.2 #ifdef SEAICE_ITD
218 torge 1.5 DO it=1,nITD
219 dimitri 1.2 DO j=1-OLy,sNy+OLy
220     DO i=1-OLx,sNx+OLx
221 torge 1.5 HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
222 dimitri 1.2 ENDDO
223     ENDDO
224     #endif
225 dimitri 1.1 CALL SEAICE_ADVECTION(
226     I GAD_HEFF, SEAICEadvSchHeff,
227     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
228     I uTrans, vTrans, HEFF(1-OLx,1-OLy,bi,bj), recip_heff,
229     O gFld, afx, afy,
230     I bi, bj, myTime, myIter, myThid )
231     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
232     C- Add tendency due to diffusion
233     CALL SEAICE_DIFFUSION(
234     I GAD_HEFF, SEAICEdiffKhHeff, ONE,
235     I HEFF(1-OLx,1-OLy,bi,bj), HEFFM,
236     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
237     U gFld,
238     I bi, bj, myTime, myIter, myThid )
239     ENDIF
240     C now do the "explicit" time step
241     DO j=1,sNy
242     DO i=1,sNx
243     HEFF(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
244     & HEFF(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
245     & )
246     ENDDO
247     ENDDO
248 dimitri 1.2 #ifdef SEAICE_ITD
249     DO j=1-OLy,sNy+OLy
250     DO i=1-OLx,sNx+OLx
251 torge 1.5 HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
252 dimitri 1.2 ENDDO
253     ENDDO
254     ENDDO
255     #endif
256 dimitri 1.1 ENDIF
257    
258     C-- Fractional area
259     IF ( SEAICEadvArea ) THEN
260 dimitri 1.2 #ifdef SEAICE_ITD
261 torge 1.5 DO it=1,nITD
262 dimitri 1.2 DO j=1-OLy,sNy+OLy
263     DO i=1-OLx,sNx+OLx
264 torge 1.5 AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
265 dimitri 1.2 ENDDO
266     ENDDO
267     #endif
268 dimitri 1.1 CALL SEAICE_ADVECTION(
269     I GAD_AREA, SEAICEadvSchArea,
270     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
271     I uTrans, vTrans, AREA(1-OLx,1-OLy,bi,bj), recip_heff,
272     O gFld, afx, afy,
273     I bi, bj, myTime, myIter, myThid )
274     IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
275     C- Add tendency due to diffusion
276     CALL SEAICE_DIFFUSION(
277     I GAD_AREA, SEAICEdiffKhArea, ONE,
278     I AREA(1-OLx,1-OLy,bi,bj), HEFFM,
279     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
280     U gFld,
281     I bi, bj, myTime, myIter, myThid )
282     ENDIF
283     C now do the "explicit" time step
284     DO j=1,sNy
285     DO i=1,sNx
286     AREA(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
287     & AREA(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
288     & )
289     ENDDO
290     ENDDO
291 dimitri 1.2 #ifdef SEAICE_ITD
292     DO j=1-OLy,sNy+OLy
293     DO i=1-OLx,sNx+OLx
294 torge 1.5 AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
295 dimitri 1.2 ENDDO
296     ENDDO
297     ENDDO
298     #endif
299 dimitri 1.1 ENDIF
300    
301     C-- Effective Snow Thickness (Volume)
302     IF ( SEAICEadvSnow ) THEN
303 dimitri 1.2 #ifdef SEAICE_ITD
304 torge 1.5 DO it=1,nITD
305 dimitri 1.2 DO j=1-OLy,sNy+OLy
306     DO i=1-OLx,sNx+OLx
307 torge 1.5 HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
308 dimitri 1.2 ENDDO
309     ENDDO
310     #endif
311 dimitri 1.1 CALL SEAICE_ADVECTION(
312     I GAD_SNOW, SEAICEadvSchSnow,
313     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
314     I uTrans, vTrans, HSNOW(1-OLx,1-OLy,bi,bj), recip_heff,
315     O gFld, afx, afy,
316     I bi, bj, myTime, myIter, myThid )
317     IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
318     C-- Add tendency due to diffusion
319     CALL SEAICE_DIFFUSION(
320     I GAD_SNOW, SEAICEdiffKhSnow, ONE,
321     I HSNOW(1-OLx,1-OLy,bi,bj), HEFFM,
322     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
323     U gFld,
324     I bi, bj, myTime, myIter, myThid )
325     ENDIF
326     C now do the "explicit" time step
327     DO j=1,sNy
328     DO i=1,sNx
329     HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
330     & HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
331     & )
332     ENDDO
333     ENDDO
334 dimitri 1.2 #ifdef SEAICE_ITD
335     DO j=1-OLy,sNy+OLy
336     DO i=1-OLx,sNx+OLx
337 torge 1.5 HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
338 dimitri 1.2 ENDDO
339     ENDDO
340     ENDDO
341     #endif
342 dimitri 1.1 ENDIF
343    
344     #ifdef SEAICE_VARIABLE_SALINITY
345     C-- Effective Sea Ice Salinity (Mass of salt)
346     IF ( SEAICEadvSalt ) THEN
347     CALL SEAICE_ADVECTION(
348     I GAD_SALT, SEAICEadvSchSalt,
349     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
350     I uTrans, vTrans, HSALT(1-OLx,1-OLy,bi,bj), recip_heff,
351     O gFld, afx, afy,
352     I bi, bj, myTime, myIter, myThid )
353     IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
354     C-- Add tendency due to diffusion
355     CALL SEAICE_DIFFUSION(
356     I GAD_SALT, SEAICEdiffKhSalt, ONE,
357     I HSALT(1-OLx,1-OLy,bi,bj), HEFFM,
358     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
359     U gFld,
360     I bi, bj, myTime, myIter, myThid )
361     ENDIF
362     C now do the "explicit" time step
363     DO j=1,sNy
364     DO i=1,sNx
365     HSALT(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
366     & HSALT(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
367     & )
368     ENDDO
369     ENDDO
370     ENDIF
371     #endif /* SEAICE_VARIABLE_SALINITY */
372    
373     #ifdef ALLOW_SITRACER
374     C-- Sea Ice Tracers
375     DO iTr = 1, SItrNumInUse
376     IF ( (SEAICEadvHEFF.AND.(SItrMate(iTr).EQ.'HEFF')).OR.
377     & (SEAICEadvAREA.AND.(SItrMate(iTr).EQ.'AREA')) ) THEN
378     C-- scale to effective value
379     IF (SItrMate(iTr).EQ.'HEFF') THEN
380     SEAICEadvSchSItr=SEAICEadvSchHEFF
381     SEAICEdiffKhSItr=SEAICEdiffKhHEFF
382     DO j=1-OLy,sNy+OLy
383     DO i=1-OLx,sNx+OLx
384     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
385     & SItracer(i,j,bi,bj,iTr) * hEffNm1(i,j,bi,bj)
386     ENDDO
387     ENDDO
388     c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
389     ELSE
390     SEAICEadvSchSItr=SEAICEadvSchAREA
391     SEAICEdiffKhSItr=SEAICEdiffKhAREA
392     DO j=1-OLy,sNy+OLy
393     DO i=1-OLx,sNx+OLx
394     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
395     & SItracer(i,j,bi,bj,iTr) * areaNm1(i,j,bi,bj)
396     ENDDO
397     ENDDO
398     ENDIF
399     C-- store a couple things
400     DO j=1-OLy,sNy+OLy
401     DO i=1-OLx,sNx+OLx
402     #ifdef ALLOW_SITRACER_ADVCAP
403     C-- store previous value for spurious maxima treament
404     SItrPrev(i,j,bi,bj)=SItracer(i,j,bi,bj,iTr)
405     #endif
406     #ifdef ALLOW_SITRACER_DEBUG_DIAG
407     diagArray(I,J,2+(iTr-1)*5) = SItrExt(i,j,bi,bj)
408     #endif
409     ENDDO
410     ENDDO
411     C-- compute advective tendency
412     CALL SEAICE_ADVECTION(
413     I GAD_SITR+iTr-1, SEAICEadvSchSItr,
414     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
415     I uTrans, vTrans, SItrExt(1-OLx,1-OLy,bi,bj),
416     I recip_heff,
417     O gFld, afx, afy,
418     I bi, bj, myTime, myIter, myThid )
419     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
420     C-- add diffusive tendency
421     CALL SEAICE_DIFFUSION(
422     I GAD_SITR+iTr-1, SEAICEdiffKhSItr, ONE,
423     I SItrExt(1-OLx,1-OLy,bi,bj), HEFFM,
424     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
425     U gFld,
426     I bi, bj, myTime, myIter, myThid )
427     ENDIF
428     C-- apply tendency
429     DO j=1,sNy
430     DO i=1,sNx
431     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
432     & SItrExt(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) )
433     ENDDO
434     ENDDO
435     C-- scale back to actual value, or move effective value to ocean bucket
436     IF (SItrMate(iTr).EQ.'HEFF') THEN
437     DO j=1,sNy
438     DO i=1,sNx
439     if (HEFF(I,J,bi,bj).GE.siEps) then
440     SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/HEFF(I,J,bi,bj)
441     SItrBucket(i,j,bi,bj,iTr)=0. _d 0
442     else
443     SItracer(i,j,bi,bj,iTr)=0. _d 0
444     SItrBucket(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)
445     endif
446     #ifdef ALLOW_SITRACER_ADVCAP
447     C hack to try avoid 'spontaneous generation' of maxima, which supposedly would
448     C occur less frequently if we advected SItr with uXheff instead SItrXheff with u
449     tmpscal1=max(SItrPrev(i,j,bi,bj),
450     & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
451     & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
452     tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
453     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
454     SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
455     & +tmpscal2*HEFF(I,J,bi,bj)
456     #endif
457     C treat case of potential negative value
458     if (HEFF(I,J,bi,bj).GE.siEps) then
459     tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
460     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
461     SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
462     & +HEFF(I,J,bi,bj)*tmpscal1
463     endif
464     #ifdef ALLOW_SITRACER_DEBUG_DIAG
465     diagArray(I,J,1+(iTr-1)*5)= - SItrBucket(i,j,bi,bj,iTr)
466     & *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
467     tmpscal1= ( HEFF(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)
468     & + SItrBucket(i,j,bi,bj,iTr) )*HEFFM(I,J,bi,bj)
469     diagArray(I,J,2+(iTr-1)*5)= tmpscal1-diagArray(I,J,2+(iTr-1)*5)
470     diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
471     & SEAICE_deltaTtherm * gFld(i,j)
472     #endif
473     ENDDO
474     ENDDO
475     c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
476     ELSE
477     DO j=1,sNy
478     DO i=1,sNx
479     if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
480     SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/AREA(I,J,bi,bj)
481     else
482     SItracer(i,j,bi,bj,iTr)=0. _d 0
483     endif
484     SItrBucket(i,j,bi,bj,iTr)=0. _d 0
485     #ifdef ALLOW_SITRACER_ADVCAP
486     tmpscal1=max(SItrPrev(i,j,bi,bj),
487     & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
488     & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
489     tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
490     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
491     #endif
492     C treat case of potential negative value
493     if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
494     tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
495     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
496     endif
497     #ifdef ALLOW_SITRACER_DEBUG_DIAG
498     diagArray(I,J,1+(iTr-1)*5)= 0. _d 0
499     diagArray(I,J,2+(iTr-1)*5)= - diagArray(I,J,2+(iTr-1)*5)
500     & + AREA(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)*HEFFM(I,J,bi,bj)
501     diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
502     & SEAICE_deltaTtherm * gFld(i,j)
503     #endif
504     ENDDO
505     ENDDO
506     ENDIF
507     C--
508     ENDIF
509     ENDDO
510     #ifdef ALLOW_SITRACER_DEBUG_DIAG
511     c CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG2 ',0,Nr,2,bi,bj,myThid)
512     #endif
513     #endif /* ALLOW_SITRACER */
514    
515     C--- end bi,bj loops
516     ENDDO
517     ENDDO
518    
519     ELSE
520     C-- if not multiDimAdvection
521    
522     IF ( SEAICEadvHEff ) THEN
523     #ifdef ALLOW_AUTODIFF_TAMC
524     CADJ STORE heff = comlev1, key = ikey_dynamics, kind=isbyte
525     #endif
526 dimitri 1.2 #ifdef SEAICE_ITD
527 torge 1.5 DO it=1,nITD
528 dimitri 1.2 DO bj=myByLo(myThid),myByHi(myThid)
529     DO bi=myBxLo(myThid),myBxHi(myThid)
530     DO j=1-OLy,sNy+OLy
531     DO i=1-OLx,sNx+OLx
532 torge 1.5 HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
533 dimitri 1.2 ENDDO
534     ENDDO
535     ENDDO
536     ENDDO
537     #endif
538 dimitri 1.1 CALL ADVECT( uc, vc, hEff, fldNm1, HEFFM, myThid )
539     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
540     C- Add tendency due to diffusion
541     DO bj=myByLo(myThid),myByHi(myThid)
542     DO bi=myBxLo(myThid),myBxHi(myThid)
543     CALL SEAICE_DIFFUSION(
544     I GAD_HEFF, SEAICEdiffKhHeff, SEAICE_deltaTtherm,
545     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
546     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
547     U HEFF(1-OLx,1-OLy,bi,bj),
548     I bi, bj, myTime, myIter, myThid )
549     ENDDO
550     ENDDO
551     ENDIF
552 dimitri 1.2 #ifdef SEAICE_ITD
553     DO bj=myByLo(myThid),myByHi(myThid)
554     DO bi=myBxLo(myThid),myBxHi(myThid)
555     DO j=1-OLy,sNy+OLy
556     DO i=1-OLx,sNx+OLx
557 torge 1.5 HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
558 dimitri 1.2 ENDDO
559     ENDDO
560     ENDDO
561     ENDDO
562     ENDDO
563     #endif
564 dimitri 1.1 ENDIF
565     IF ( SEAICEadvArea ) THEN
566     #ifdef ALLOW_AUTODIFF_TAMC
567     CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte
568     #endif
569 dimitri 1.2 #ifdef SEAICE_ITD
570 torge 1.5 DO it=1,nITD
571 dimitri 1.2 DO bj=myByLo(myThid),myByHi(myThid)
572     DO bi=myBxLo(myThid),myBxHi(myThid)
573     DO j=1-OLy,sNy+OLy
574     DO i=1-OLx,sNx+OLx
575 torge 1.5 AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
576 dimitri 1.2 ENDDO
577     ENDDO
578     ENDDO
579     ENDDO
580     #endif
581 dimitri 1.1 CALL ADVECT( uc, vc, area, fldNm1, HEFFM, myThid )
582     IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
583     C- Add tendency due to diffusion
584     DO bj=myByLo(myThid),myByHi(myThid)
585     DO bi=myBxLo(myThid),myBxHi(myThid)
586     CALL SEAICE_DIFFUSION(
587     I GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm,
588     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
589     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
590     U Area(1-OLx,1-OLy,bi,bj),
591     I bi, bj, myTime, myIter, myThid )
592     ENDDO
593     ENDDO
594     ENDIF
595 dimitri 1.2 #ifdef SEAICE_ITD
596     DO bj=myByLo(myThid),myByHi(myThid)
597     DO bi=myBxLo(myThid),myBxHi(myThid)
598     DO j=1-OLy,sNy+OLy
599     DO i=1-OLx,sNx+OLx
600 torge 1.5 AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
601 dimitri 1.2 ENDDO
602     ENDDO
603     ENDDO
604     ENDDO
605     ENDDO
606     #endif
607 dimitri 1.1 ENDIF
608     IF ( SEAICEadvSnow ) THEN
609     #ifdef ALLOW_AUTODIFF_TAMC
610     CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte
611     #endif
612 dimitri 1.2 #ifdef SEAICE_ITD
613 torge 1.5 DO it=1,nITD
614 dimitri 1.2 DO bj=myByLo(myThid),myByHi(myThid)
615     DO bi=myBxLo(myThid),myBxHi(myThid)
616     DO j=1-OLy,sNy+OLy
617     DO i=1-OLx,sNx+OLx
618 torge 1.5 HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
619 dimitri 1.2 ENDDO
620     ENDDO
621     ENDDO
622     ENDDO
623     #endif
624 dimitri 1.1 CALL ADVECT( uc, vc, HSNOW, fldNm1, HEFFM, myThid )
625     IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
626     C- Add tendency due to diffusion
627     DO bj=myByLo(myThid),myByHi(myThid)
628     DO bi=myBxLo(myThid),myBxHi(myThid)
629     CALL SEAICE_DIFFUSION(
630     I GAD_SNOW, SEAICEdiffKhSnow, SEAICE_deltaTtherm,
631     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
632     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
633     U HSNOW(1-OLx,1-OLy,bi,bj),
634     I bi, bj, myTime, myIter, myThid )
635     ENDDO
636     ENDDO
637     ENDIF
638 dimitri 1.2 #ifdef SEAICE_ITD
639     DO bj=myByLo(myThid),myByHi(myThid)
640     DO bi=myBxLo(myThid),myBxHi(myThid)
641     DO j=1-OLy,sNy+OLy
642     DO i=1-OLx,sNx+OLx
643 torge 1.5 HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
644 dimitri 1.2 ENDDO
645     ENDDO
646     ENDDO
647     ENDDO
648     ENDDO
649     #endif
650 dimitri 1.1 ENDIF
651    
652     #ifdef SEAICE_VARIABLE_SALINITY
653     IF ( SEAICEadvSalt ) THEN
654     #ifdef ALLOW_AUTODIFF_TAMC
655     CADJ STORE hsalt = comlev1, key = ikey_dynamics, kind=isbyte
656     #endif
657     CALL ADVECT( uc, vc, HSALT, fldNm1, HEFFM, myThid )
658     IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
659     C- Add tendency due to diffusion
660     DO bj=myByLo(myThid),myByHi(myThid)
661     DO bi=myBxLo(myThid),myBxHi(myThid)
662     CALL SEAICE_DIFFUSION(
663     I GAD_SALT, SEAICEdiffKhSalt, SEAICE_deltaTtherm,
664     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
665     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
666     U HSALT(1-OLx,1-OLy,bi,bj),
667     I bi, bj, myTime, myIter, myThid )
668     ENDDO
669     ENDDO
670     ENDIF
671     ENDIF
672     #endif /* SEAICE_VARIABLE_SALINITY */
673    
674     C-- end if multiDimAdvection
675     ENDIF
676    
677     #ifdef ALLOW_AUTODIFF_TAMC
678     CADJ STORE AREA = comlev1, key = ikey_dynamics, kind=isbyte
679     #endif
680     IF ( .NOT. usePW79thermodynamics ) THEN
681     C Hiblers "ridging function": Do it now if not in seaice_growth
682     C in principle we should add a "real" ridging function here (or
683     C somewhere after doing the advection)
684     DO bj=myByLo(myThid),myByHi(myThid)
685     DO bi=myBxLo(myThid),myBxHi(myThid)
686     DO j=1-OLy,sNy+OLy
687     DO i=1-OLx,sNx+OLx
688     AREA(I,J,bi,bj) = MIN(ONE,AREA(I,J,bi,bj))
689     ENDDO
690     ENDDO
691     ENDDO
692     ENDDO
693     ENDIF
694    
695     RETURN
696     END
  ViewVC Help
Powered by ViewVC 1.1.22