/[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.4 - (hide annotations) (download)
Wed Mar 27 18:59:52 2013 UTC (12 years, 4 months ago) by torge
Branch: MAIN
Changes since 1.3: +2 -29 lines
updating my MITgcm_contrib directory to include latest changes on main branch;
settings are to run a 1D test szenario with ITD code and 7 categories

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

  ViewVC Help
Powered by ViewVC 1.1.22