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

Annotation of /MITgcm_contrib/torge/itd/code/seaice_advection.F

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


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

  ViewVC Help
Powered by ViewVC 1.1.22