/[MITgcm]/MITgcm/pkg/seaice/seaice_advection.F
ViewVC logotype

Contents of /MITgcm/pkg/seaice/seaice_advection.F

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

Revision 1.33 - (show annotations) (download)
Mon Oct 20 03:20:57 2014 UTC (9 years, 6 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65g, HEAD
Changes since 1.32: +4 -1 lines 
- ECCO_OPTIONS.h is needed when including ecco_cost.h, ecco.h
- AUTODIFF_OPTIONS.h is needed when including tamc.h, tamc_keys.h
- CTRL_OPTIONS.h is needed when including ctrl.h, etc

- pkg/seaice/seaice_cost*.F : clean up CPP brackets
- SEAICE_SIZE.h : replace ALLOW_AUTODIFF_TAMC with ALLOW_AUTODIFF to
  avoid needing AUTODIFF_OPTIONS.h anytime SEAICE_SIZE.h is included
  (it seems that THSICE_SIZE.h, PTRACERS_SIZE.h have the same issue...)
1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advection.F,v 1.32 2012/11/09 22:15:18 heimbach Exp $
2 C $Name: $
3
4 #include "SEAICE_OPTIONS.h"
5 #ifdef ALLOW_GENERIC_ADVDIFF
6 # include "GAD_OPTIONS.h"
7 #endif
8 #ifdef ALLOW_AUTODIFF
9 # include "AUTODIFF_OPTIONS.h"
10 #endif
11
12 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
13 CBOP
14 C !ROUTINE: SEAICE_ADVECTION
15
16 C !INTERFACE: ==========================================================
17 SUBROUTINE SEAICE_ADVECTION(
18 I tracerIdentity,
19 I advectionScheme,
20 I uFld, vFld, uTrans, vTrans, iceFld, r_hFld,
21 O gFld, afx, afy,
22 I bi, bj, myTime, myIter, myThid)
23
24 C !DESCRIPTION:
25 C Calculates the tendency of a sea-ice field due to advection.
26 C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
27 C and can only be used for the non-linear advection schemes such as the
28 C direct-space-time method and flux-limiters.
29 C
30 C This routine is an adaption of the GAD_ADVECTION for 2D-fields.
31 C for Area, effective thickness or other "extensive" sea-ice field,
32 C the contribution iceFld*div(u) (that is present in gad_advection)
33 C is not included here.
34 C
35 C The algorithm is as follows:
36 C \begin{itemize}
37 C \item{$\theta^{(n+1/2)} = \theta^{(n)}
38 C - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
39 C \item{$\theta^{(n+2/2)} = \theta^{(n+1/2)}
40 C - \Delta t \partial_y (v\theta^{(n+1/2)}) + \theta^{(n)} \partial_y v$}
41 C \item{$G_\theta = ( \theta^{(n+2/2)} - \theta^{(n)} )/\Delta t$}
42 C \end{itemize}
43 C
44 C The tendency (output) is over-written by this routine.
45
46 C !USES: ===============================================================
47 IMPLICIT NONE
48 #include "SIZE.h"
49 #include "EEPARAMS.h"
50 #include "PARAMS.h"
51 #include "GRID.h"
52 #include "SEAICE_SIZE.h"
53 #include "SEAICE_PARAMS.h"
54 #ifdef ALLOW_GENERIC_ADVDIFF
55 # include "GAD.h"
56 #endif
57 #ifdef ALLOW_AUTODIFF_TAMC
58 # include "tamc.h"
59 # include "tamc_keys.h"
60 # ifdef ALLOW_PTRACERS
61 # include "PTRACERS_SIZE.h"
62 # endif
63 #endif
64 #ifdef ALLOW_EXCH2
65 #include "W2_EXCH2_SIZE.h"
66 #include "W2_EXCH2_TOPOLOGY.h"
67 #endif /* ALLOW_EXCH2 */
68 LOGICAL extensiveFld
69 PARAMETER ( extensiveFld = .TRUE. )
70
71 C !INPUT PARAMETERS: ===================================================
72 C tracerIdentity :: tracer identifier
73 C advectionScheme :: advection scheme to use (Horizontal plane)
74 C extensiveFld :: indicates to advect an "extensive" type of ice field
75 C uFld :: velocity, zonal component
76 C vFld :: velocity, meridional component
77 C uTrans,vTrans :: volume transports at U,V points
78 C iceFld :: sea-ice field
79 C r_hFld :: reciprocal of ice-thickness (only used for "intensive"
80 C type of sea-ice field)
81 C bi,bj :: tile indices
82 C myTime :: current time
83 C myIter :: iteration number
84 C myThid :: my Thread Id number
85 INTEGER tracerIdentity
86 INTEGER advectionScheme
87 _RL uFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 _RL vFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
89 _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
90 _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
91 _RL iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
92 _RL r_hFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
93 INTEGER bi,bj
94 _RL myTime
95 INTEGER myIter
96 INTEGER myThid
97
98 C !OUTPUT PARAMETERS: ==================================================
99 C gFld :: tendency array
100 C afx :: horizontal advective flux, x direction
101 C afy :: horizontal advective flux, y direction
102 _RL gFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
103 _RL afx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
104 _RL afy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
105
106 C !LOCAL VARIABLES: ====================================================
107 C maskLocW :: 2-D array for mask at West points
108 C maskLocS :: 2-D array for mask at South points
109 C iMin,iMax, :: loop range for called routines
110 C jMin,jMax :: loop range for called routines
111 C [iMin,iMax]Upd :: loop range to update sea-ice field
112 C [jMin,jMax]Upd :: loop range to update sea-ice field
113 C i,j,k :: loop indices
114 C af :: 2-D array for horizontal advective flux
115 C localTij :: 2-D array, temporary local copy of sea-ice fld
116 C calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
117 C calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
118 C interiorOnly :: only update the interior of myTile, but not the edges
119 C overlapOnly :: only update the edges of myTile, but not the interior
120 C nipass :: number of passes in multi-dimensional method
121 C ipass :: number of the current pass being made
122 C myTile :: variables used to determine which cube face
123 C nCFace :: owns a tile for cube grid runs using
124 C :: multi-dim advection.
125 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
126 C msgBuf :: Informational/error message buffer
127 _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
128 _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129 INTEGER iMin,iMax,jMin,jMax
130 INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
131 INTEGER i,j,k
132 _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
133 _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
134 LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
135 LOGICAL interiorOnly, overlapOnly
136 INTEGER nipass,ipass
137 INTEGER nCFace
138 LOGICAL N_edge, S_edge, E_edge, W_edge
139 CHARACTER*(MAX_LEN_MBUF) msgBuf
140 #ifdef ALLOW_EXCH2
141 INTEGER myTile
142 #endif
143 #ifdef ALLOW_DIAGNOSTICS
144 CHARACTER*8 diagName
145 CHARACTER*4 SEAICE_DIAG_SUFX, diagSufx
146 EXTERNAL SEAICE_DIAG_SUFX
147 #endif
148 LOGICAL dBug
149 INTEGER ioUnit
150 _RL tmpFac
151 CEOP
152
153 #ifdef ALLOW_AUTODIFF_TAMC
154 act0 = tracerIdentity - 1
155 max0 = maxpass
156 act1 = bi - myBxLo(myThid)
157 max1 = myBxHi(myThid) - myBxLo(myThid) + 1
158 act2 = bj - myByLo(myThid)
159 max2 = myByHi(myThid) - myByLo(myThid) + 1
160 act3 = myThid - 1
161 max3 = nTx*nTy
162 act4 = ikey_dynamics - 1
163 igadkey = (act0 + 1)
164 & + act1*max0
165 & + act2*max0*max1
166 & + act3*max0*max1*max2
167 & + act4*max0*max1*max2*max3
168 if (tracerIdentity.GT.maxpass) then
169 WRITE(msgBuf,'(A,2I3)')
170 & 'SEAICE_ADVECTION: maxpass < tracerIdentity ',
171 & maxpass, tracerIdentity
172 CALL PRINT_ERROR( msgBuf, myThid )
173 STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
174 endif
175 #endif /* ALLOW_AUTODIFF_TAMC */
176
177 #ifdef ALLOW_DIAGNOSTICS
178 C-- Set diagnostic suffix for the current tracer
179 IF ( useDiagnostics ) THEN
180 diagSufx = SEAICE_DIAG_SUFX( tracerIdentity, myThid )
181 ENDIF
182 #endif
183
184 ioUnit = standardMessageUnit
185 dBug = debugLevel.GE.debLevC
186 & .AND. myIter.EQ.nIter0
187 & .AND. ( tracerIdentity.EQ.GAD_HEFF .OR.
188 & tracerIdentity.EQ.GAD_QICE2 )
189 c & .AND. tracerIdentity.EQ.GAD_HEFF
190
191 C-- Set up work arrays with valid (i.e. not NaN) values
192 C These inital values do not alter the numerical results. They
193 C just ensure that all memory references are to valid floating
194 C point numbers. This prevents spurious hardware signals due to
195 C uninitialised but inert locations.
196 #ifdef ALLOW_AUTODIFF_TAMC
197 DO j=1-OLy,sNy+OLy
198 DO i=1-OLx,sNx+OLx
199 localTij(i,j) = 0. _d 0
200 ENDDO
201 ENDDO
202 #endif
203
204 C-- Set tile-specific parameters for horizontal fluxes
205 IF (useCubedSphereExchange) THEN
206 nipass=3
207 #ifdef ALLOW_AUTODIFF_TAMC
208 IF ( nipass.GT.maxcube ) THEN
209 WRITE(msgBuf,'(A)')
210 & 'SEAICE_ADVECTION: maxcube needs to be =3; check tamc.h '
211 CALL PRINT_ERROR( msgBuf, myThid )
212 STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
213 ENDIF
214 #endif
215 #ifdef ALLOW_EXCH2
216 myTile = W2_myTileList(bi,bj)
217 nCFace = exch2_myFace(myTile)
218 N_edge = exch2_isNedge(myTile).EQ.1
219 S_edge = exch2_isSedge(myTile).EQ.1
220 E_edge = exch2_isEedge(myTile).EQ.1
221 W_edge = exch2_isWedge(myTile).EQ.1
222 #else
223 nCFace = bi
224 N_edge = .TRUE.
225 S_edge = .TRUE.
226 E_edge = .TRUE.
227 W_edge = .TRUE.
228 #endif
229 ELSE
230 nipass=2
231 nCFace = bi
232 N_edge = .FALSE.
233 S_edge = .FALSE.
234 E_edge = .FALSE.
235 W_edge = .FALSE.
236 ENDIF
237
238 iMin = 1-OLx
239 iMax = sNx+OLx
240 jMin = 1-OLy
241 jMax = sNy+OLy
242
243 k = 1
244 C-- Start of k loop for horizontal fluxes
245 #ifdef ALLOW_AUTODIFF_TAMC
246 CADJ STORE iceFld =
247 CADJ & comlev1_bibj_k_gadice, key=igadkey, byte=isbyte
248 #endif /* ALLOW_AUTODIFF_TAMC */
249
250 C Content of CALC_COMMON_FACTORS, adapted for 2D fields
251 C-- Get temporary terms used by tendency routines
252
253 C-- Make local copy of sea-ice field and mask West & South
254 DO j=1-OLy,sNy+OLy
255 DO i=1-OLx,sNx+OLx
256 localTij(i,j)=iceFld(i,j)
257 #ifdef ALLOW_OBCS
258 maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
259 maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
260 #else /* ALLOW_OBCS */
261 maskLocW(i,j) = _maskW(i,j,k,bi,bj)
262 maskLocS(i,j) = _maskS(i,j,k,bi,bj)
263 #endif /* ALLOW_OBCS */
264 ENDDO
265 ENDDO
266
267 #ifdef ALLOW_AUTODIFF_TAMC
268 C- Initialise Advective flux in X & Y
269 DO j=1-OLy,sNy+OLy
270 DO i=1-OLx,sNx+OLx
271 afx(i,j) = 0.
272 afy(i,j) = 0.
273 ENDDO
274 ENDDO
275 #endif
276
277 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
278 IF (useCubedSphereExchange) THEN
279 withSigns = .FALSE.
280 CALL FILL_CS_CORNER_UV_RS(
281 & withSigns, maskLocW,maskLocS, bi,bj, myThid )
282 ENDIF
283 cph-exch2#endif
284
285 C-- Multiple passes for different directions on different tiles
286 C-- For cube need one pass for each of red, green and blue axes.
287 DO ipass=1,nipass
288 #ifdef ALLOW_AUTODIFF_TAMC
289 passkey = ipass + (igadkey-1)*maxpass
290 IF (nipass .GT. maxpass) THEN
291 WRITE(msgBuf,'(A,2I3)')
292 & 'SEAICE_ADVECTION: nipass > max[ass. check tamc.h ',
293 & nipass, maxpass
294 CALL PRINT_ERROR( msgBuf, myThid )
295 STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
296 ENDIF
297 #endif /* ALLOW_AUTODIFF_TAMC */
298
299 interiorOnly = .FALSE.
300 overlapOnly = .FALSE.
301 IF (useCubedSphereExchange) THEN
302 C-- CubedSphere : pass 3 times, with partial update of local seaice field
303 IF (ipass.EQ.1) THEN
304 overlapOnly = MOD(nCFace,3).EQ.0
305 interiorOnly = MOD(nCFace,3).NE.0
306 calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
307 calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
308 ELSEIF (ipass.EQ.2) THEN
309 overlapOnly = MOD(nCFace,3).EQ.2
310 calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
311 calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
312 ELSE
313 calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
314 calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
315 ENDIF
316 ELSE
317 C-- not CubedSphere
318 calc_fluxes_X = MOD(ipass,2).EQ.1
319 calc_fluxes_Y = .NOT.calc_fluxes_X
320 ENDIF
321 IF (dBug.AND.bi.EQ.3 ) WRITE(ioUnit,*)'ICE_adv:',tracerIdentity,
322 & ipass,calc_fluxes_X,calc_fluxes_Y,overlapOnly,interiorOnly
323
324 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
325 C-- X direction
326
327 #ifdef ALLOW_AUTODIFF_TAMC
328 CADJ STORE localTij(:,:) =
329 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
330 # ifndef DISABLE_MULTIDIM_ADVECTION
331 CADJ STORE af(:,:) =
332 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
333 # endif
334 #endif /* ALLOW_AUTODIFF_TAMC */
335 C
336 IF (calc_fluxes_X) THEN
337
338 C- Do not compute fluxes if
339 C a) needed in overlap only
340 C and b) the overlap of myTile are not cube-face Edges
341 IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
342
343 C- Advective flux in X
344 DO j=1-OLy,sNy+OLy
345 DO i=1-OLx,sNx+OLx
346 af(i,j) = 0.
347 ENDDO
348 ENDDO
349
350 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
351 C- Internal exchange for calculations in X
352 IF ( useCubedSphereExchange .AND.
353 & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
354 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
355 & localTij, bi,bj, myThid )
356 ENDIF
357 cph-exch2#endif
358
359 #ifdef ALLOW_AUTODIFF_TAMC
360 # ifndef DISABLE_MULTIDIM_ADVECTION
361 CADJ STORE localTij(:,:) =
362 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
363 # endif
364 #endif /* ALLOW_AUTODIFF_TAMC */
365
366 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
367 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
368 CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
369 I SEAICE_deltaTtherm, uTrans, uFld, localTij,
370 O af, myThid )
371 IF ( dBug .AND. bi.EQ.3 ) THEN
372 i=MIN(12,sNx)
373 j=MIN(11,sNy)
374 WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv: xFx=', af(i+1,j),
375 & localTij(i,j), uTrans(i+1,j), af(i+1,j)/uTrans(i+1,j)
376 ENDIF
377 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
378 CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE.,
379 I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
380 O af, myThid )
381 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
382 CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE.,
383 I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
384 O af, myThid )
385 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
386 CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE.,
387 I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
388 O af, myThid )
389 #ifndef ALLOW_AUTODIFF_TAMC
390 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
391 CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE.,
392 I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
393 O af, myThid )
394 #endif
395 ELSE
396 WRITE(msgBuf,'(A,I3,A)')
397 & 'SEAICE_ADVECTION: adv. scheme ', advectionScheme,
398 & ' incompatibale with multi-dim. adv.'
399 CALL PRINT_ERROR( msgBuf, myThid )
400 STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
401 ENDIF
402
403 C-- Advective flux in X : done
404 ENDIF
405
406 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
407 C-- Internal exchange for next calculations in Y
408 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
409 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
410 & localTij, bi,bj, myThid )
411 ENDIF
412 cph-exch2#endif
413
414 C- Update the local seaice field where needed:
415
416 C update in overlap-Only
417 IF ( overlapOnly ) THEN
418 iMinUpd = 1-OLx+1
419 iMaxUpd = sNx+OLx-1
420 C-- notes: these 2 lines below have no real effect (because recip_hFac=0
421 C in corner region) but safer to keep them.
422 IF ( W_edge ) iMinUpd = 1
423 IF ( E_edge ) iMaxUpd = sNx
424
425 IF ( S_edge .AND. extensiveFld ) THEN
426 DO j=1-OLy,0
427 DO i=iMinUpd,iMaxUpd
428 localTij(i,j)=localTij(i,j)
429 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
430 & *recip_rA(i,j,bi,bj)
431 & *( af(i+1,j)-af(i,j)
432 & )
433 ENDDO
434 ENDDO
435 ELSEIF ( S_edge ) THEN
436 DO j=1-OLy,0
437 DO i=iMinUpd,iMaxUpd
438 localTij(i,j)=localTij(i,j)
439 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
440 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
441 & *( (af(i+1,j)-af(i,j))
442 & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
443 & )
444 ENDDO
445 ENDDO
446 ENDIF
447 IF ( N_edge .AND. extensiveFld ) THEN
448 DO j=sNy+1,sNy+OLy
449 DO i=iMinUpd,iMaxUpd
450 localTij(i,j)=localTij(i,j)
451 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
452 & *recip_rA(i,j,bi,bj)
453 & *( af(i+1,j)-af(i,j)
454 & )
455 ENDDO
456 ENDDO
457 ELSEIF ( N_edge ) THEN
458 DO j=sNy+1,sNy+OLy
459 DO i=iMinUpd,iMaxUpd
460 localTij(i,j)=localTij(i,j)
461 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
462 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
463 & *( (af(i+1,j)-af(i,j))
464 & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
465 & )
466 ENDDO
467 ENDDO
468 ENDIF
469 C-- keep advective flux (for diagnostics)
470 IF ( S_edge ) THEN
471 DO j=1-OLy,0
472 DO i=1-OLx+1,sNx+OLx
473 afx(i,j) = af(i,j)
474 ENDDO
475 ENDDO
476 ENDIF
477 IF ( N_edge ) THEN
478 DO j=sNy+1,sNy+OLy
479 DO i=1-OLx+1,sNx+OLx
480 afx(i,j) = af(i,j)
481 ENDDO
482 ENDDO
483 ENDIF
484
485 ELSE
486 C do not only update the overlap
487 jMinUpd = 1-OLy
488 jMaxUpd = sNy+OLy
489 IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
490 IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
491 IF ( extensiveFld ) THEN
492 DO j=jMinUpd,jMaxUpd
493 DO i=1-OLx+1,sNx+OLx-1
494 localTij(i,j)=localTij(i,j)
495 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
496 & *recip_rA(i,j,bi,bj)
497 & *( af(i+1,j)-af(i,j)
498 & )
499 ENDDO
500 ENDDO
501 ELSE
502 DO j=jMinUpd,jMaxUpd
503 DO i=1-OLx+1,sNx+OLx-1
504 localTij(i,j)=localTij(i,j)
505 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
506 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
507 & *( (af(i+1,j)-af(i,j))
508 & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
509 & )
510 ENDDO
511 ENDDO
512 ENDIF
513 C-- keep advective flux (for diagnostics)
514 DO j=jMinUpd,jMaxUpd
515 DO i=1-OLx+1,sNx+OLx
516 afx(i,j) = af(i,j)
517 ENDDO
518 ENDDO
519
520 C- end if/else update overlap-Only
521 ENDIF
522
523 C-- End of X direction
524 ENDIF
525
526 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
527 C-- Y direction
528
529 #ifdef ALLOW_AUTODIFF_TAMC
530 # ifndef DISABLE_MULTIDIM_ADVECTION
531 CADJ STORE localTij(:,:) =
532 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
533 CADJ STORE af(:,:) =
534 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
535 # endif
536 #endif /* ALLOW_AUTODIFF_TAMC */
537
538 IF (calc_fluxes_Y) THEN
539
540 C- Do not compute fluxes if
541 C a) needed in overlap only
542 C and b) the overlap of myTile are not cube-face edges
543 IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
544
545 C- Advective flux in Y
546 DO j=1-OLy,sNy+OLy
547 DO i=1-OLx,sNx+OLx
548 af(i,j) = 0.
549 ENDDO
550 ENDDO
551
552 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
553 C- Internal exchange for calculations in Y
554 IF ( useCubedSphereExchange .AND.
555 & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
556 CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
557 & localTij, bi,bj, myThid )
558 ENDIF
559 cph-exch2#endif
560
561 #ifdef ALLOW_AUTODIFF_TAMC
562 #ifndef DISABLE_MULTIDIM_ADVECTION
563 CADJ STORE localTij(:,:) =
564 CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
565 #endif
566 #endif /* ALLOW_AUTODIFF_TAMC */
567
568 IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
569 & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
570 CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
571 I SEAICE_deltaTtherm, vTrans, vFld, localTij,
572 O af, myThid )
573 IF ( dBug .AND. bi.EQ.3 ) THEN
574 i=MIN(12,sNx)
575 j=MIN(11,sNy)
576 WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv: yFx=', af(i,j+1),
577 & localTij(i,j), vTrans(i,j+1), af(i,j+1)/vTrans(i,j+1)
578 ENDIF
579 ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
580 CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE.,
581 I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
582 O af, myThid )
583 ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
584 CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE.,
585 I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
586 O af, myThid )
587 ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
588 CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE.,
589 I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
590 O af, myThid )
591 #ifndef ALLOW_AUTODIFF_TAMC
592 ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
593 CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE.,
594 I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
595 O af, myThid )
596 #endif
597 ELSE
598 WRITE(msgBuf,'(A,I3,A)')
599 & 'SEAICE_ADVECTION: adv. scheme ', advectionScheme,
600 & ' incompatibale with multi-dim. adv.'
601 CALL PRINT_ERROR( msgBuf, myThid )
602 STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
603 ENDIF
604
605 C- Advective flux in Y : done
606 ENDIF
607
608 cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
609 C- Internal exchange for next calculations in X
610 IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
611 CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
612 & localTij, bi,bj, myThid )
613 ENDIF
614 cph-exch2#endif
615
616 C- Update the local seaice field where needed:
617
618 C update in overlap-Only
619 IF ( overlapOnly ) THEN
620 jMinUpd = 1-OLy+1
621 jMaxUpd = sNy+OLy-1
622 C- notes: these 2 lines below have no real effect (because recip_hFac=0
623 C in corner region) but safer to keep them.
624 IF ( S_edge ) jMinUpd = 1
625 IF ( N_edge ) jMaxUpd = sNy
626
627 IF ( W_edge .AND. extensiveFld ) THEN
628 DO j=jMinUpd,jMaxUpd
629 DO i=1-OLx,0
630 localTij(i,j)=localTij(i,j)
631 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
632 & *recip_rA(i,j,bi,bj)
633 & *( af(i,j+1)-af(i,j)
634 & )
635 ENDDO
636 ENDDO
637 ELSEIF ( W_edge ) THEN
638 DO j=jMinUpd,jMaxUpd
639 DO i=1-OLx,0
640 localTij(i,j)=localTij(i,j)
641 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
642 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
643 & *( (af(i,j+1)-af(i,j))
644 & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
645 & )
646 ENDDO
647 ENDDO
648 ENDIF
649 IF ( E_edge .AND. extensiveFld ) THEN
650 DO j=jMinUpd,jMaxUpd
651 DO i=sNx+1,sNx+OLx
652 localTij(i,j)=localTij(i,j)
653 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
654 & *recip_rA(i,j,bi,bj)
655 & *( af(i,j+1)-af(i,j)
656 & )
657 ENDDO
658 ENDDO
659 ELSEIF ( E_edge ) THEN
660 DO j=jMinUpd,jMaxUpd
661 DO i=sNx+1,sNx+OLx
662 localTij(i,j)=localTij(i,j)
663 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
664 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
665 & *( (af(i,j+1)-af(i,j))
666 & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
667 & )
668 ENDDO
669 ENDDO
670 ENDIF
671 C-- keep advective flux (for diagnostics)
672 IF ( W_edge ) THEN
673 DO j=1-OLy+1,sNy+OLy
674 DO i=1-OLx,0
675 afy(i,j) = af(i,j)
676 ENDDO
677 ENDDO
678 ENDIF
679 IF ( E_edge ) THEN
680 DO j=1-OLy+1,sNy+OLy
681 DO i=sNx+1,sNx+OLx
682 afy(i,j) = af(i,j)
683 ENDDO
684 ENDDO
685 ENDIF
686
687 ELSE
688 C do not only update the overlap
689 iMinUpd = 1-OLx
690 iMaxUpd = sNx+OLx
691 IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
692 IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
693 IF ( extensiveFld ) THEN
694 DO j=1-OLy+1,sNy+OLy-1
695 DO i=iMinUpd,iMaxUpd
696 localTij(i,j)=localTij(i,j)
697 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
698 & *recip_rA(i,j,bi,bj)
699 & *( af(i,j+1)-af(i,j)
700 & )
701 ENDDO
702 ENDDO
703 ELSE
704 DO j=1-OLy+1,sNy+OLy-1
705 DO i=iMinUpd,iMaxUpd
706 localTij(i,j)=localTij(i,j)
707 & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
708 & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
709 & *( (af(i,j+1)-af(i,j))
710 & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
711 & )
712 ENDDO
713 ENDDO
714 ENDIF
715 C-- keep advective flux (for diagnostics)
716 DO j=1-OLy+1,sNy+OLy
717 DO i=iMinUpd,iMaxUpd
718 afy(i,j) = af(i,j)
719 ENDDO
720 ENDDO
721
722 C end if/else update overlap-Only
723 ENDIF
724
725 C-- End of Y direction
726 ENDIF
727
728 C-- End of ipass loop
729 ENDDO
730
731 C- explicit advection is done ; store tendency in gFld:
732 DO j=1-OLy,sNy+OLy
733 DO i=1-OLx,sNx+OLx
734 gFld(i,j)=(localTij(i,j)-iceFld(i,j))/SEAICE_deltaTtherm
735 ENDDO
736 ENDDO
737 IF ( dBug .AND. bi.EQ.3 ) THEN
738 i=MIN(12,sNx)
739 j=MIN(11,sNy)
740 tmpFac= SEAICE_deltaTtherm*recip_rA(i,j,bi,bj)
741 WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv:',
742 & afx(i,j)*tmpFac,afx(i+1,j)*tmpFac,
743 & afy(i,j)*tmpFac,afy(i,j+1)*tmpFac
744 ENDIF
745
746 #ifdef ALLOW_DIAGNOSTICS
747 IF ( useDiagnostics ) THEN
748 diagName = 'ADVx'//diagSufx
749 CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
750 diagName = 'ADVy'//diagSufx
751 CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
752 ENDIF
753 #endif
754
755 #ifdef ALLOW_DEBUG
756 IF ( debugLevel .GE. debLevC
757 & .AND. tracerIdentity.EQ.GAD_HEFF
758 & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
759 & .AND. nPx.EQ.1 .AND. nPy.EQ.1
760 & .AND. useCubedSphereExchange ) THEN
761 CALL DEBUG_CS_CORNER_UV( ' afx,afy from SEAICE_ADVECTION',
762 & afx,afy, k, standardMessageUnit,bi,bj,myThid )
763 ENDIF
764 #endif /* ALLOW_DEBUG */
765
766 RETURN
767 END
  ViewVC Help
Powered by ViewVC 1.1.22