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

Contents 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 - (show 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 C $Header: /u/gcmpack/MITgcm_contrib/torge/itd/code/seaice_advdiff.F,v 1.4 2013/03/27 18:59:52 torge 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 #ifdef SEAICE_ITD
57 C it :: Loop counter for ice thickness categories
58 #endif
59 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 #ifdef SEAICE_ITD
69 INTEGER it
70 #endif
71 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 # ifdef SEAICE_VARIABLE_SALINITY
192 CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,
193 CADJ & key = itmpkey, kind=isbyte
194 # endif
195 #endif /* ALLOW_AUTODIFF_TAMC */
196
197 DO j=1-OLy,sNy+OLy
198 DO i=1-OLx,sNx+OLx
199 #ifdef ALLOW_SITRACER
200 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 #ifdef SEAICE_ITD
218 DO it=1,nITD
219 DO j=1-OLy,sNy+OLy
220 DO i=1-OLx,sNx+OLx
221 HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
222 ENDDO
223 ENDDO
224 #endif
225 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 #ifdef SEAICE_ITD
249 DO j=1-OLy,sNy+OLy
250 DO i=1-OLx,sNx+OLx
251 HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
252 ENDDO
253 ENDDO
254 ENDDO
255 #endif
256 ENDIF
257
258 C-- Fractional area
259 IF ( SEAICEadvArea ) THEN
260 #ifdef SEAICE_ITD
261 DO it=1,nITD
262 DO j=1-OLy,sNy+OLy
263 DO i=1-OLx,sNx+OLx
264 AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
265 ENDDO
266 ENDDO
267 #endif
268 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 #ifdef SEAICE_ITD
292 DO j=1-OLy,sNy+OLy
293 DO i=1-OLx,sNx+OLx
294 AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
295 ENDDO
296 ENDDO
297 ENDDO
298 #endif
299 ENDIF
300
301 C-- Effective Snow Thickness (Volume)
302 IF ( SEAICEadvSnow ) THEN
303 #ifdef SEAICE_ITD
304 DO it=1,nITD
305 DO j=1-OLy,sNy+OLy
306 DO i=1-OLx,sNx+OLx
307 HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
308 ENDDO
309 ENDDO
310 #endif
311 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 #ifdef SEAICE_ITD
335 DO j=1-OLy,sNy+OLy
336 DO i=1-OLx,sNx+OLx
337 HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
338 ENDDO
339 ENDDO
340 ENDDO
341 #endif
342 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 #ifdef SEAICE_ITD
527 DO it=1,nITD
528 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 HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
533 ENDDO
534 ENDDO
535 ENDDO
536 ENDDO
537 #endif
538 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 #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 HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
558 ENDDO
559 ENDDO
560 ENDDO
561 ENDDO
562 ENDDO
563 #endif
564 ENDIF
565 IF ( SEAICEadvArea ) THEN
566 #ifdef ALLOW_AUTODIFF_TAMC
567 CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte
568 #endif
569 #ifdef SEAICE_ITD
570 DO it=1,nITD
571 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 AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
576 ENDDO
577 ENDDO
578 ENDDO
579 ENDDO
580 #endif
581 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 #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 AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
601 ENDDO
602 ENDDO
603 ENDDO
604 ENDDO
605 ENDDO
606 #endif
607 ENDIF
608 IF ( SEAICEadvSnow ) THEN
609 #ifdef ALLOW_AUTODIFF_TAMC
610 CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte
611 #endif
612 #ifdef SEAICE_ITD
613 DO it=1,nITD
614 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 HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
619 ENDDO
620 ENDDO
621 ENDDO
622 ENDDO
623 #endif
624 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 #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 HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
644 ENDDO
645 ENDDO
646 ENDDO
647 ENDDO
648 ENDDO
649 #endif
650 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