/[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.2 - (hide annotations) (download)
Fri Apr 27 22:25:23 2012 UTC (13 years, 3 months ago) by dimitri
Branch: MAIN
Changes since 1.1: +154 -0 lines
first check in of itd code

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

  ViewVC Help
Powered by ViewVC 1.1.22