/[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.3 - (hide annotations) (download)
Mon Dec 10 22:19:49 2012 UTC (12 years, 7 months ago) by torge
Branch: MAIN
Changes since 1.2: +5 -1 lines
include updates from main branch

1 torge 1.3 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advdiff.F,v 1.62 2012/10/23 21:28:08 heimbach 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 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 torge 1.3 # ifdef SEAICE_VARIABLE_SALINITY
198     CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
199     CADJ & key = itmpkey, kind=isbyte
200     # endif
201 dimitri 1.1 #endif /* ALLOW_AUTODIFF_TAMC */
202    
203     DO j=1-OLy,sNy+OLy
204     DO i=1-OLx,sNx+OLx
205     #if ( defined (SEAICE_GROWTH_LEGACY) || defined (ALLOW_SITRACER) )
206     hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj)
207     areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj)
208     #endif
209     recip_heff(i,j) = 1. _d 0
210     ENDDO
211     ENDDO
212    
213     C- Calculate "volume transports" through tracer cell faces.
214     DO j=1-OLy,sNy+OLy
215     DO i=1-OLx,sNx+OLx
216     uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j,bi,bj)
217     vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j,bi,bj)
218     ENDDO
219     ENDDO
220    
221     C-- Effective Thickness (Volume)
222     IF ( SEAICEadvHeff ) THEN
223 dimitri 1.2 CToM<<<
224     #ifdef SEAICE_ITD
225     DO k=1,nITD
226     DO j=1-OLy,sNy+OLy
227     DO i=1-OLx,sNx+OLx
228     HEFF(i,j,bi,bj)=HEFFITD(i,j,k,bi,bj)
229     ENDDO
230     ENDDO
231     #endif
232     C>>>ToM
233 dimitri 1.1 CALL SEAICE_ADVECTION(
234     I GAD_HEFF, SEAICEadvSchHeff,
235     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
236     I uTrans, vTrans, HEFF(1-OLx,1-OLy,bi,bj), recip_heff,
237     O gFld, afx, afy,
238     I bi, bj, myTime, myIter, myThid )
239     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
240     C- Add tendency due to diffusion
241     CALL SEAICE_DIFFUSION(
242     I GAD_HEFF, SEAICEdiffKhHeff, ONE,
243     I HEFF(1-OLx,1-OLy,bi,bj), HEFFM,
244     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
245     U gFld,
246     I bi, bj, myTime, myIter, myThid )
247     ENDIF
248     C now do the "explicit" time step
249     DO j=1,sNy
250     DO i=1,sNx
251     HEFF(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
252     & HEFF(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
253     & )
254     ENDDO
255     ENDDO
256 dimitri 1.2 CToM<<<
257     #ifdef SEAICE_ITD
258     DO j=1-OLy,sNy+OLy
259     DO i=1-OLx,sNx+OLx
260     HEFFITD(i,j,k,bi,bj)=HEFF(i,j,bi,bj)
261     ENDDO
262     ENDDO
263     ENDDO
264     #endif
265     C>>>ToM
266 dimitri 1.1 ENDIF
267    
268     C-- Fractional area
269     IF ( SEAICEadvArea ) THEN
270 dimitri 1.2 CToM<<<
271     #ifdef SEAICE_ITD
272     DO k=1,nITD
273     DO j=1-OLy,sNy+OLy
274     DO i=1-OLx,sNx+OLx
275     AREA(i,j,bi,bj)=AREAITD(i,j,k,bi,bj)
276     ENDDO
277     ENDDO
278     #endif
279     C>>>ToM
280 dimitri 1.1 CALL SEAICE_ADVECTION(
281     I GAD_AREA, SEAICEadvSchArea,
282     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
283     I uTrans, vTrans, AREA(1-OLx,1-OLy,bi,bj), recip_heff,
284     O gFld, afx, afy,
285     I bi, bj, myTime, myIter, myThid )
286     IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
287     C- Add tendency due to diffusion
288     CALL SEAICE_DIFFUSION(
289     I GAD_AREA, SEAICEdiffKhArea, ONE,
290     I AREA(1-OLx,1-OLy,bi,bj), HEFFM,
291     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
292     U gFld,
293     I bi, bj, myTime, myIter, myThid )
294     ENDIF
295     C now do the "explicit" time step
296     DO j=1,sNy
297     DO i=1,sNx
298     AREA(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
299     & AREA(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
300     & )
301     ENDDO
302     ENDDO
303 dimitri 1.2 CToM<<<
304     #ifdef SEAICE_ITD
305     DO j=1-OLy,sNy+OLy
306     DO i=1-OLx,sNx+OLx
307     AREAITD(i,j,k,bi,bj)=AREA(i,j,bi,bj)
308     ENDDO
309     ENDDO
310     ENDDO
311     #endif
312     C>>>ToM
313 dimitri 1.1 ENDIF
314    
315     C-- Effective Snow Thickness (Volume)
316     IF ( SEAICEadvSnow ) THEN
317 dimitri 1.2 CToM<<<
318     #ifdef SEAICE_ITD
319     DO k=1,nITD
320     DO j=1-OLy,sNy+OLy
321     DO i=1-OLx,sNx+OLx
322     HSNOW(i,j,bi,bj)=HSNOWITD(i,j,k,bi,bj)
323     ENDDO
324     ENDDO
325     #endif
326     C>>>ToM
327 dimitri 1.1 CALL SEAICE_ADVECTION(
328     I GAD_SNOW, SEAICEadvSchSnow,
329     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
330     I uTrans, vTrans, HSNOW(1-OLx,1-OLy,bi,bj), recip_heff,
331     O gFld, afx, afy,
332     I bi, bj, myTime, myIter, myThid )
333     IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
334     C-- Add tendency due to diffusion
335     CALL SEAICE_DIFFUSION(
336     I GAD_SNOW, SEAICEdiffKhSnow, ONE,
337     I HSNOW(1-OLx,1-OLy,bi,bj), HEFFM,
338     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
339     U gFld,
340     I bi, bj, myTime, myIter, myThid )
341     ENDIF
342     C now do the "explicit" time step
343     DO j=1,sNy
344     DO i=1,sNx
345     HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
346     & HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
347     & )
348     ENDDO
349     ENDDO
350 dimitri 1.2 CToM<<<
351     #ifdef SEAICE_ITD
352     DO j=1-OLy,sNy+OLy
353     DO i=1-OLx,sNx+OLx
354     HSNOWITD(i,j,k,bi,bj)=HSNOW(i,j,bi,bj)
355     ENDDO
356     ENDDO
357     ENDDO
358     #endif
359     C>>>ToM
360 dimitri 1.1 ENDIF
361    
362     #ifdef SEAICE_VARIABLE_SALINITY
363     C-- Effective Sea Ice Salinity (Mass of salt)
364     IF ( SEAICEadvSalt ) THEN
365     CALL SEAICE_ADVECTION(
366     I GAD_SALT, SEAICEadvSchSalt,
367     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
368     I uTrans, vTrans, HSALT(1-OLx,1-OLy,bi,bj), recip_heff,
369     O gFld, afx, afy,
370     I bi, bj, myTime, myIter, myThid )
371     IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
372     C-- Add tendency due to diffusion
373     CALL SEAICE_DIFFUSION(
374     I GAD_SALT, SEAICEdiffKhSalt, ONE,
375     I HSALT(1-OLx,1-OLy,bi,bj), HEFFM,
376     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
377     U gFld,
378     I bi, bj, myTime, myIter, myThid )
379     ENDIF
380     C now do the "explicit" time step
381     DO j=1,sNy
382     DO i=1,sNx
383     HSALT(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
384     & HSALT(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
385     & )
386     ENDDO
387     ENDDO
388     ENDIF
389     #endif /* SEAICE_VARIABLE_SALINITY */
390    
391     #ifdef ALLOW_SITRACER
392     C-- Sea Ice Tracers
393     DO iTr = 1, SItrNumInUse
394     IF ( (SEAICEadvHEFF.AND.(SItrMate(iTr).EQ.'HEFF')).OR.
395     & (SEAICEadvAREA.AND.(SItrMate(iTr).EQ.'AREA')) ) THEN
396     C-- scale to effective value
397     IF (SItrMate(iTr).EQ.'HEFF') THEN
398     SEAICEadvSchSItr=SEAICEadvSchHEFF
399     SEAICEdiffKhSItr=SEAICEdiffKhHEFF
400     DO j=1-OLy,sNy+OLy
401     DO i=1-OLx,sNx+OLx
402     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
403     & SItracer(i,j,bi,bj,iTr) * hEffNm1(i,j,bi,bj)
404     ENDDO
405     ENDDO
406     c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
407     ELSE
408     SEAICEadvSchSItr=SEAICEadvSchAREA
409     SEAICEdiffKhSItr=SEAICEdiffKhAREA
410     DO j=1-OLy,sNy+OLy
411     DO i=1-OLx,sNx+OLx
412     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
413     & SItracer(i,j,bi,bj,iTr) * areaNm1(i,j,bi,bj)
414     ENDDO
415     ENDDO
416     ENDIF
417     C-- store a couple things
418     DO j=1-OLy,sNy+OLy
419     DO i=1-OLx,sNx+OLx
420     #ifdef ALLOW_SITRACER_ADVCAP
421     C-- store previous value for spurious maxima treament
422     SItrPrev(i,j,bi,bj)=SItracer(i,j,bi,bj,iTr)
423     #endif
424     #ifdef ALLOW_SITRACER_DEBUG_DIAG
425     diagArray(I,J,2+(iTr-1)*5) = SItrExt(i,j,bi,bj)
426     #endif
427     ENDDO
428     ENDDO
429     C-- compute advective tendency
430     CALL SEAICE_ADVECTION(
431     I GAD_SITR+iTr-1, SEAICEadvSchSItr,
432     I uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
433     I uTrans, vTrans, SItrExt(1-OLx,1-OLy,bi,bj),
434     I recip_heff,
435     O gFld, afx, afy,
436     I bi, bj, myTime, myIter, myThid )
437     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
438     C-- add diffusive tendency
439     CALL SEAICE_DIFFUSION(
440     I GAD_SITR+iTr-1, SEAICEdiffKhSItr, ONE,
441     I SItrExt(1-OLx,1-OLy,bi,bj), HEFFM,
442     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
443     U gFld,
444     I bi, bj, myTime, myIter, myThid )
445     ENDIF
446     C-- apply tendency
447     DO j=1,sNy
448     DO i=1,sNx
449     SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
450     & SItrExt(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) )
451     ENDDO
452     ENDDO
453     C-- scale back to actual value, or move effective value to ocean bucket
454     IF (SItrMate(iTr).EQ.'HEFF') THEN
455     DO j=1,sNy
456     DO i=1,sNx
457     if (HEFF(I,J,bi,bj).GE.siEps) then
458     SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/HEFF(I,J,bi,bj)
459     SItrBucket(i,j,bi,bj,iTr)=0. _d 0
460     else
461     SItracer(i,j,bi,bj,iTr)=0. _d 0
462     SItrBucket(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)
463     endif
464     #ifdef ALLOW_SITRACER_ADVCAP
465     C hack to try avoid 'spontaneous generation' of maxima, which supposedly would
466     C occur less frequently if we advected SItr with uXheff instead SItrXheff with u
467     tmpscal1=max(SItrPrev(i,j,bi,bj),
468     & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
469     & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
470     tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
471     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
472     SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
473     & +tmpscal2*HEFF(I,J,bi,bj)
474     #endif
475     C treat case of potential negative value
476     if (HEFF(I,J,bi,bj).GE.siEps) then
477     tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
478     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
479     SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
480     & +HEFF(I,J,bi,bj)*tmpscal1
481     endif
482     #ifdef ALLOW_SITRACER_DEBUG_DIAG
483     diagArray(I,J,1+(iTr-1)*5)= - SItrBucket(i,j,bi,bj,iTr)
484     & *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
485     tmpscal1= ( HEFF(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)
486     & + SItrBucket(i,j,bi,bj,iTr) )*HEFFM(I,J,bi,bj)
487     diagArray(I,J,2+(iTr-1)*5)= tmpscal1-diagArray(I,J,2+(iTr-1)*5)
488     diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
489     & SEAICE_deltaTtherm * gFld(i,j)
490     #endif
491     ENDDO
492     ENDDO
493     c TAF? ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
494     ELSE
495     DO j=1,sNy
496     DO i=1,sNx
497     if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
498     SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/AREA(I,J,bi,bj)
499     else
500     SItracer(i,j,bi,bj,iTr)=0. _d 0
501     endif
502     SItrBucket(i,j,bi,bj,iTr)=0. _d 0
503     #ifdef ALLOW_SITRACER_ADVCAP
504     tmpscal1=max(SItrPrev(i,j,bi,bj),
505     & SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
506     & SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
507     tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
508     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
509     #endif
510     C treat case of potential negative value
511     if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
512     tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
513     SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
514     endif
515     #ifdef ALLOW_SITRACER_DEBUG_DIAG
516     diagArray(I,J,1+(iTr-1)*5)= 0. _d 0
517     diagArray(I,J,2+(iTr-1)*5)= - diagArray(I,J,2+(iTr-1)*5)
518     & + AREA(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)*HEFFM(I,J,bi,bj)
519     diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
520     & SEAICE_deltaTtherm * gFld(i,j)
521     #endif
522     ENDDO
523     ENDDO
524     ENDIF
525     C--
526     ENDIF
527     ENDDO
528     #ifdef ALLOW_SITRACER_DEBUG_DIAG
529     c CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG2 ',0,Nr,2,bi,bj,myThid)
530     #endif
531     #endif /* ALLOW_SITRACER */
532    
533     C--- end bi,bj loops
534     ENDDO
535     ENDDO
536    
537     ELSE
538     C-- if not multiDimAdvection
539    
540     Cold This has to be done to comply with the time stepping in advect.F:
541     Cold Making sure that the following routines see the different
542     Cold time levels correctly
543     Cold At the end of the routine ADVECT,
544     Cold timelevel 1 is updated with advection contribution
545     Cold and diffusion contribution
546     Cold (which was computed in DIFFUS on timelevel 3)
547     Cold timelevel 2 is the previous timelevel 1
548     Cold timelevel 3 is the total diffusion tendency * deltaT
549     Cold (empty if no diffusion)
550     C-- This is what remains from old 3-level storage of AREA & HEFF: still
551     C needed for SEAICE_GROWTH, Legacy branch. Left old comments here above.
552     #ifdef SEAICE_GROWTH_LEGACY
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     hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj)
558     areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj)
559     ENDDO
560     ENDDO
561     ENDDO
562     ENDDO
563     #endif
564    
565     IF ( SEAICEadvHEff ) THEN
566     #ifdef ALLOW_AUTODIFF_TAMC
567     CADJ STORE heff = comlev1, key = ikey_dynamics, kind=isbyte
568     #endif
569 dimitri 1.2 CToM<<<
570     #ifdef SEAICE_ITD
571     DO k=1,nITD
572     DO bj=myByLo(myThid),myByHi(myThid)
573     DO bi=myBxLo(myThid),myBxHi(myThid)
574     DO j=1-OLy,sNy+OLy
575     DO i=1-OLx,sNx+OLx
576     HEFF(i,j,bi,bj)=HEFFITD(i,j,k,bi,bj)
577     ENDDO
578     ENDDO
579     ENDDO
580     ENDDO
581     #endif
582     C>>>ToM
583 dimitri 1.1 CALL ADVECT( uc, vc, hEff, fldNm1, HEFFM, myThid )
584     IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
585     C- Add tendency due to diffusion
586     DO bj=myByLo(myThid),myByHi(myThid)
587     DO bi=myBxLo(myThid),myBxHi(myThid)
588     CALL SEAICE_DIFFUSION(
589     I GAD_HEFF, SEAICEdiffKhHeff, SEAICE_deltaTtherm,
590     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
591     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
592     U HEFF(1-OLx,1-OLy,bi,bj),
593     I bi, bj, myTime, myIter, myThid )
594     ENDDO
595     ENDDO
596     ENDIF
597 dimitri 1.2 CToM<<<
598     #ifdef SEAICE_ITD
599     DO bj=myByLo(myThid),myByHi(myThid)
600     DO bi=myBxLo(myThid),myBxHi(myThid)
601     DO j=1-OLy,sNy+OLy
602     DO i=1-OLx,sNx+OLx
603     HEFFITD(i,j,k,bi,bj)=HEFF(i,j,bi,bj)
604     ENDDO
605     ENDDO
606     ENDDO
607     ENDDO
608     ENDDO
609     #endif
610     C>>>ToM
611 dimitri 1.1 ENDIF
612     IF ( SEAICEadvArea ) THEN
613     #ifdef ALLOW_AUTODIFF_TAMC
614     CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte
615     #endif
616 dimitri 1.2 CToM<<<
617     #ifdef SEAICE_ITD
618     DO k=1,nITD
619     DO bj=myByLo(myThid),myByHi(myThid)
620     DO bi=myBxLo(myThid),myBxHi(myThid)
621     DO j=1-OLy,sNy+OLy
622     DO i=1-OLx,sNx+OLx
623     AREA(i,j,bi,bj)=AREAITD(i,j,k,bi,bj)
624     ENDDO
625     ENDDO
626     ENDDO
627     ENDDO
628     #endif
629     C>>>ToM
630 dimitri 1.1 CALL ADVECT( uc, vc, area, fldNm1, HEFFM, myThid )
631     IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
632     C- Add tendency due to diffusion
633     DO bj=myByLo(myThid),myByHi(myThid)
634     DO bi=myBxLo(myThid),myBxHi(myThid)
635     CALL SEAICE_DIFFUSION(
636     I GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm,
637     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
638     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
639     U Area(1-OLx,1-OLy,bi,bj),
640     I bi, bj, myTime, myIter, myThid )
641     ENDDO
642     ENDDO
643     ENDIF
644 dimitri 1.2 CToM<<<
645     #ifdef SEAICE_ITD
646     DO bj=myByLo(myThid),myByHi(myThid)
647     DO bi=myBxLo(myThid),myBxHi(myThid)
648     DO j=1-OLy,sNy+OLy
649     DO i=1-OLx,sNx+OLx
650     AREAITD(i,j,k,bi,bj)=AREA(i,j,bi,bj)
651     ENDDO
652     ENDDO
653     ENDDO
654     ENDDO
655     ENDDO
656     #endif
657     C>>>ToM
658 dimitri 1.1 ENDIF
659     IF ( SEAICEadvSnow ) THEN
660     #ifdef ALLOW_AUTODIFF_TAMC
661     CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte
662     #endif
663 dimitri 1.2 CToM<<<
664     #ifdef SEAICE_ITD
665     DO k=1,nITD
666     DO bj=myByLo(myThid),myByHi(myThid)
667     DO bi=myBxLo(myThid),myBxHi(myThid)
668     DO j=1-OLy,sNy+OLy
669     DO i=1-OLx,sNx+OLx
670     HSNOW(i,j,bi,bj)=HSNOWITD(i,j,k,bi,bj)
671     ENDDO
672     ENDDO
673     ENDDO
674     ENDDO
675     #endif
676     C>>>ToM
677 dimitri 1.1 CALL ADVECT( uc, vc, HSNOW, fldNm1, HEFFM, myThid )
678     IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
679     C- Add tendency due to diffusion
680     DO bj=myByLo(myThid),myByHi(myThid)
681     DO bi=myBxLo(myThid),myBxHi(myThid)
682     CALL SEAICE_DIFFUSION(
683     I GAD_SNOW, SEAICEdiffKhSnow, SEAICE_deltaTtherm,
684     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
685     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
686     U HSNOW(1-OLx,1-OLy,bi,bj),
687     I bi, bj, myTime, myIter, myThid )
688     ENDDO
689     ENDDO
690     ENDIF
691 dimitri 1.2 CToM<<<
692     #ifdef SEAICE_ITD
693     DO bj=myByLo(myThid),myByHi(myThid)
694     DO bi=myBxLo(myThid),myBxHi(myThid)
695     DO j=1-OLy,sNy+OLy
696     DO i=1-OLx,sNx+OLx
697     HSNOWITD(i,j,k,bi,bj)=HSNOW(i,j,bi,bj)
698     ENDDO
699     ENDDO
700     ENDDO
701     ENDDO
702     ENDDO
703     #endif
704     C>>>ToM
705 dimitri 1.1 ENDIF
706    
707     #ifdef SEAICE_VARIABLE_SALINITY
708     IF ( SEAICEadvSalt ) THEN
709     #ifdef ALLOW_AUTODIFF_TAMC
710     CADJ STORE hsalt = comlev1, key = ikey_dynamics, kind=isbyte
711     #endif
712     CALL ADVECT( uc, vc, HSALT, fldNm1, HEFFM, myThid )
713     IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
714     C- Add tendency due to diffusion
715     DO bj=myByLo(myThid),myByHi(myThid)
716     DO bi=myBxLo(myThid),myBxHi(myThid)
717     CALL SEAICE_DIFFUSION(
718     I GAD_SALT, SEAICEdiffKhSalt, SEAICE_deltaTtherm,
719     I fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
720     I xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
721     U HSALT(1-OLx,1-OLy,bi,bj),
722     I bi, bj, myTime, myIter, myThid )
723     ENDDO
724     ENDDO
725     ENDIF
726     ENDIF
727     #endif /* SEAICE_VARIABLE_SALINITY */
728    
729     C-- end if multiDimAdvection
730     ENDIF
731    
732     #ifdef ALLOW_AUTODIFF_TAMC
733     CADJ STORE AREA = comlev1, key = ikey_dynamics, kind=isbyte
734     #endif
735     IF ( .NOT. usePW79thermodynamics ) THEN
736     C Hiblers "ridging function": Do it now if not in seaice_growth
737     C in principle we should add a "real" ridging function here (or
738     C somewhere after doing the advection)
739     DO bj=myByLo(myThid),myByHi(myThid)
740     DO bi=myBxLo(myThid),myBxHi(myThid)
741     DO j=1-OLy,sNy+OLy
742     DO i=1-OLx,sNx+OLx
743     AREA(I,J,bi,bj) = MIN(ONE,AREA(I,J,bi,bj))
744     ENDDO
745     ENDDO
746     ENDDO
747     ENDDO
748     ENDIF
749    
750     RETURN
751     END

  ViewVC Help
Powered by ViewVC 1.1.22