/[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.1 - (hide annotations) (download)
Wed Oct 24 21:48:53 2012 UTC (12 years, 9 months ago) by torge
Branch: MAIN
adding "#include SEAICE_SIZE.h" where necessary,
i.e. where SEAICE_PARAMS.h is included but SEAICE_SIZE.h wasn't

1 torge 1.1 C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advection.F,v 1.29 2011/06/24 22:23:26 jmc Exp $
2     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     _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
124     _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
125     INTEGER iMin,iMax,jMin,jMax
126     INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
127     INTEGER i,j,k
128     _RL af (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
129     _RL localTij(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
130     LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
131     LOGICAL interiorOnly, overlapOnly
132     INTEGER nipass,ipass
133     INTEGER nCFace
134     LOGICAL N_edge, S_edge, E_edge, W_edge
135     #ifdef ALLOW_EXCH2
136     INTEGER myTile
137     #endif
138     #ifdef ALLOW_AUTODIFF_TAMC
139     C msgBuf :: Informational/error message buffer
140     CHARACTER*(MAX_LEN_MBUF) msgBuf
141     #endif
142     #ifdef ALLOW_DIAGNOSTICS
143     CHARACTER*8 diagName
144     CHARACTER*4 SEAICE_DIAG_SUFX, diagSufx
145     EXTERNAL SEAICE_DIAG_SUFX
146     #endif
147     LOGICAL dBug
148     INTEGER ioUnit
149     _RL tmpFac
150     CEOP
151    
152     #ifdef ALLOW_AUTODIFF_TAMC
153     act0 = tracerIdentity - 1
154     max0 = maxpass
155     act1 = bi - myBxLo(myThid)
156     max1 = myBxHi(myThid) - myBxLo(myThid) + 1
157     act2 = bj - myByLo(myThid)
158     max2 = myByHi(myThid) - myByLo(myThid) + 1
159     act3 = myThid - 1
160     max3 = nTx*nTy
161     act4 = ikey_dynamics - 1
162     igadkey = (act0 + 1)
163     & + act1*max0
164     & + act2*max0*max1
165     & + act3*max0*max1*max2
166     & + act4*max0*max1*max2*max3
167     if (tracerIdentity.GT.maxpass) then
168     WRITE(msgBuf,'(A,2I3)')
169     & 'SEAICE_ADVECTION: maxpass < tracerIdentity ',
170     & maxpass, tracerIdentity
171     CALL PRINT_ERROR( msgBuf, myThid )
172     STOP 'ABNORMAL END: S/R SEAICE_ADVECTION'
173     endif
174     #endif /* ALLOW_AUTODIFF_TAMC */
175    
176     #ifdef ALLOW_DIAGNOSTICS
177     C-- Set diagnostic suffix for the current tracer
178     IF ( useDiagnostics ) THEN
179     diagSufx = SEAICE_DIAG_SUFX( tracerIdentity, myThid )
180     ENDIF
181     #endif
182    
183     ioUnit = standardMessageUnit
184     dBug = debugLevel.GE.debLevC
185     & .AND. myIter.EQ.nIter0
186     & .AND. ( tracerIdentity.EQ.GAD_HEFF .OR.
187     & tracerIdentity.EQ.GAD_QICE2 )
188     c & .AND. tracerIdentity.EQ.GAD_HEFF
189    
190     C-- Set up work arrays with valid (i.e. not NaN) values
191     C These inital values do not alter the numerical results. They
192     C just ensure that all memory references are to valid floating
193     C point numbers. This prevents spurious hardware signals due to
194     C uninitialised but inert locations.
195     #ifdef ALLOW_AUTODIFF_TAMC
196     DO j=1-OLy,sNy+OLy
197     DO i=1-OLx,sNx+OLx
198     localTij(i,j) = 0. _d 0
199     ENDDO
200     ENDDO
201     #endif
202    
203     C-- Set tile-specific parameters for horizontal fluxes
204     IF (useCubedSphereExchange) THEN
205     nipass=3
206     #ifdef ALLOW_AUTODIFF_TAMC
207     IF ( nipass.GT.maxcube ) STOP 'maxcube needs to be = 3'
208     #endif
209     #ifdef ALLOW_EXCH2
210     myTile = W2_myTileList(bi,bj)
211     nCFace = exch2_myFace(myTile)
212     N_edge = exch2_isNedge(myTile).EQ.1
213     S_edge = exch2_isSedge(myTile).EQ.1
214     E_edge = exch2_isEedge(myTile).EQ.1
215     W_edge = exch2_isWedge(myTile).EQ.1
216     #else
217     nCFace = bi
218     N_edge = .TRUE.
219     S_edge = .TRUE.
220     E_edge = .TRUE.
221     W_edge = .TRUE.
222     #endif
223     ELSE
224     nipass=2
225     nCFace = bi
226     N_edge = .FALSE.
227     S_edge = .FALSE.
228     E_edge = .FALSE.
229     W_edge = .FALSE.
230     ENDIF
231    
232     iMin = 1-OLx
233     iMax = sNx+OLx
234     jMin = 1-OLy
235     jMax = sNy+OLy
236    
237     k = 1
238     C-- Start of k loop for horizontal fluxes
239     #ifdef ALLOW_AUTODIFF_TAMC
240     CADJ STORE iceFld =
241     CADJ & comlev1_bibj_k_gadice, key=igadkey, byte=isbyte
242     #endif /* ALLOW_AUTODIFF_TAMC */
243    
244     C Content of CALC_COMMON_FACTORS, adapted for 2D fields
245     C-- Get temporary terms used by tendency routines
246    
247     C-- Make local copy of sea-ice field and mask West & South
248     DO j=1-OLy,sNy+OLy
249     DO i=1-OLx,sNx+OLx
250     localTij(i,j)=iceFld(i,j)
251     #ifdef ALLOW_OBCS
252     maskLocW(i,j) = _maskW(i,j,k,bi,bj)*maskInW(i,j,bi,bj)
253     maskLocS(i,j) = _maskS(i,j,k,bi,bj)*maskInS(i,j,bi,bj)
254     #else /* ALLOW_OBCS */
255     maskLocW(i,j) = _maskW(i,j,k,bi,bj)
256     maskLocS(i,j) = _maskS(i,j,k,bi,bj)
257     #endif /* ALLOW_OBCS */
258     ENDDO
259     ENDDO
260    
261     #ifdef ALLOW_AUTODIFF_TAMC
262     C- Initialise Advective flux in X & Y
263     DO j=1-OLy,sNy+OLy
264     DO i=1-OLx,sNx+OLx
265     afx(i,j) = 0.
266     afy(i,j) = 0.
267     ENDDO
268     ENDDO
269     #endif
270    
271     cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
272     IF (useCubedSphereExchange) THEN
273     withSigns = .FALSE.
274     CALL FILL_CS_CORNER_UV_RS(
275     & withSigns, maskLocW,maskLocS, bi,bj, myThid )
276     ENDIF
277     cph-exch2#endif
278    
279     C-- Multiple passes for different directions on different tiles
280     C-- For cube need one pass for each of red, green and blue axes.
281     DO ipass=1,nipass
282     #ifdef ALLOW_AUTODIFF_TAMC
283     passkey = ipass + (igadkey-1)*maxpass
284     IF (nipass .GT. maxpass) THEN
285     STOP 'SEAICE_ADVECTION: nipass > maxcube. check tamc.h'
286     ENDIF
287     #endif /* ALLOW_AUTODIFF_TAMC */
288    
289     interiorOnly = .FALSE.
290     overlapOnly = .FALSE.
291     IF (useCubedSphereExchange) THEN
292     C-- CubedSphere : pass 3 times, with partial update of local seaice field
293     IF (ipass.EQ.1) THEN
294     overlapOnly = MOD(nCFace,3).EQ.0
295     interiorOnly = MOD(nCFace,3).NE.0
296     calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
297     calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
298     ELSEIF (ipass.EQ.2) THEN
299     overlapOnly = MOD(nCFace,3).EQ.2
300     calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
301     calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
302     ELSE
303     calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
304     calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
305     ENDIF
306     ELSE
307     C-- not CubedSphere
308     calc_fluxes_X = MOD(ipass,2).EQ.1
309     calc_fluxes_Y = .NOT.calc_fluxes_X
310     ENDIF
311     IF (dBug.AND.bi.EQ.3 ) WRITE(ioUnit,*)'ICE_adv:',tracerIdentity,
312     & ipass,calc_fluxes_X,calc_fluxes_Y,overlapOnly,interiorOnly
313    
314     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
315     C-- X direction
316    
317     #ifdef ALLOW_AUTODIFF_TAMC
318     CADJ STORE localTij(:,:) =
319     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
320     # ifndef DISABLE_MULTIDIM_ADVECTION
321     CADJ STORE af(:,:) =
322     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
323     # endif
324     #endif /* ALLOW_AUTODIFF_TAMC */
325     C
326     IF (calc_fluxes_X) THEN
327    
328     C- Do not compute fluxes if
329     C a) needed in overlap only
330     C and b) the overlap of myTile are not cube-face Edges
331     IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN
332    
333     C- Advective flux in X
334     DO j=1-Oly,sNy+Oly
335     DO i=1-Olx,sNx+Olx
336     af(i,j) = 0.
337     ENDDO
338     ENDDO
339    
340     cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
341     C- Internal exchange for calculations in X
342     IF ( useCubedSphereExchange .AND.
343     & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
344     CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
345     & localTij, bi,bj, myThid )
346     ENDIF
347     cph-exch2#endif
348    
349     #ifdef ALLOW_AUTODIFF_TAMC
350     # ifndef DISABLE_MULTIDIM_ADVECTION
351     CADJ STORE localTij(:,:) =
352     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
353     # endif
354     #endif /* ALLOW_AUTODIFF_TAMC */
355    
356     IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
357     & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
358     CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .TRUE.,
359     I SEAICE_deltaTtherm, uTrans, uFld, localTij,
360     O af, myThid )
361     IF ( dBug .AND. bi.EQ.3 ) THEN
362     i=MIN(12,sNx)
363     j=MIN(11,sNy)
364     WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv: xFx=', af(i+1,j),
365     & localTij(i,j), uTrans(i+1,j), af(i+1,j)/uTrans(i+1,j)
366     ENDIF
367     ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
368     CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .TRUE.,
369     I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
370     O af, myThid )
371     ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
372     CALL GAD_DST3_ADV_X( bi,bj,k, .TRUE.,
373     I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
374     O af, myThid )
375     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
376     CALL GAD_DST3FL_ADV_X( bi,bj,k, .TRUE.,
377     I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
378     O af, myThid )
379     #ifndef ALLOW_AUTODIFF_TAMC
380     ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
381     CALL GAD_OS7MP_ADV_X( bi,bj,k, .TRUE.,
382     I SEAICE_deltaTtherm, uTrans, uFld, maskLocW, localTij,
383     O af, myThid )
384     #endif
385     ELSE
386     STOP
387     & 'SEAICE_ADVECTION: adv. scheme incompatibale with multi-dim'
388     ENDIF
389    
390     C-- Advective flux in X : done
391     ENDIF
392    
393     cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
394     C-- Internal exchange for next calculations in Y
395     IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
396     CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
397     & localTij, bi,bj, myThid )
398     ENDIF
399     cph-exch2#endif
400    
401     C- Update the local seaice field where needed:
402    
403     C update in overlap-Only
404     IF ( overlapOnly ) THEN
405     iMinUpd = 1-OLx+1
406     iMaxUpd = sNx+OLx-1
407     C-- notes: these 2 lines below have no real effect (because recip_hFac=0
408     C in corner region) but safer to keep them.
409     IF ( W_edge ) iMinUpd = 1
410     IF ( E_edge ) iMaxUpd = sNx
411    
412     IF ( S_edge .AND. extensiveFld ) THEN
413     DO j=1-OLy,0
414     DO i=iMinUpd,iMaxUpd
415     localTij(i,j)=localTij(i,j)
416     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
417     & *recip_rA(i,j,bi,bj)
418     & *( af(i+1,j)-af(i,j)
419     & )
420     ENDDO
421     ENDDO
422     ELSEIF ( S_edge ) 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)*r_hFld(i,j)
428     & *( (af(i+1,j)-af(i,j))
429     & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
430     & )
431     ENDDO
432     ENDDO
433     ENDIF
434     IF ( N_edge .AND. extensiveFld ) THEN
435     DO j=sNy+1,sNy+OLy
436     DO i=iMinUpd,iMaxUpd
437     localTij(i,j)=localTij(i,j)
438     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
439     & *recip_rA(i,j,bi,bj)
440     & *( af(i+1,j)-af(i,j)
441     & )
442     ENDDO
443     ENDDO
444     ELSEIF ( N_edge ) 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)*r_hFld(i,j)
450     & *( (af(i+1,j)-af(i,j))
451     & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
452     & )
453     ENDDO
454     ENDDO
455     ENDIF
456     C-- keep advective flux (for diagnostics)
457     IF ( S_edge ) THEN
458     DO j=1-OLy,0
459     DO i=1-OLx+1,sNx+OLx
460     afx(i,j) = af(i,j)
461     ENDDO
462     ENDDO
463     ENDIF
464     IF ( N_edge ) THEN
465     DO j=sNy+1,sNy+OLy
466     DO i=1-OLx+1,sNx+OLx
467     afx(i,j) = af(i,j)
468     ENDDO
469     ENDDO
470     ENDIF
471    
472     ELSE
473     C do not only update the overlap
474     jMinUpd = 1-OLy
475     jMaxUpd = sNy+OLy
476     IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
477     IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
478     IF ( extensiveFld ) THEN
479     DO j=jMinUpd,jMaxUpd
480     DO i=1-OLx+1,sNx+OLx-1
481     localTij(i,j)=localTij(i,j)
482     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
483     & *recip_rA(i,j,bi,bj)
484     & *( af(i+1,j)-af(i,j)
485     & )
486     ENDDO
487     ENDDO
488     ELSE
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)*r_hFld(i,j)
494     & *( (af(i+1,j)-af(i,j))
495     & -(uTrans(i+1,j)-uTrans(i,j))*iceFld(i,j)
496     & )
497     ENDDO
498     ENDDO
499     ENDIF
500     C-- keep advective flux (for diagnostics)
501     DO j=jMinUpd,jMaxUpd
502     DO i=1-OLx+1,sNx+OLx
503     afx(i,j) = af(i,j)
504     ENDDO
505     ENDDO
506    
507     C- end if/else update overlap-Only
508     ENDIF
509    
510     C-- End of X direction
511     ENDIF
512    
513     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
514     C-- Y direction
515    
516     #ifdef ALLOW_AUTODIFF_TAMC
517     # ifndef DISABLE_MULTIDIM_ADVECTION
518     CADJ STORE localTij(:,:) =
519     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
520     CADJ STORE af(:,:) =
521     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
522     # endif
523     #endif /* ALLOW_AUTODIFF_TAMC */
524    
525     IF (calc_fluxes_Y) THEN
526    
527     C- Do not compute fluxes if
528     C a) needed in overlap only
529     C and b) the overlap of myTile are not cube-face edges
530     IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN
531    
532     C- Advective flux in Y
533     DO j=1-OLy,sNy+OLy
534     DO i=1-OLx,sNx+OLx
535     af(i,j) = 0.
536     ENDDO
537     ENDDO
538    
539     cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
540     C- Internal exchange for calculations in Y
541     IF ( useCubedSphereExchange .AND.
542     & ( overlapOnly .OR. ipass.EQ.1 ) ) THEN
543     CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
544     & localTij, bi,bj, myThid )
545     ENDIF
546     cph-exch2#endif
547    
548     #ifdef ALLOW_AUTODIFF_TAMC
549     #ifndef DISABLE_MULTIDIM_ADVECTION
550     CADJ STORE localTij(:,:) =
551     CADJ & comlev1_bibj_k_gadice_pass, key=passkey, byte=isbyte
552     #endif
553     #endif /* ALLOW_AUTODIFF_TAMC */
554    
555     IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
556     & .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
557     CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .TRUE.,
558     I SEAICE_deltaTtherm, vTrans, vFld, localTij,
559     O af, myThid )
560     IF ( dBug .AND. bi.EQ.3 ) THEN
561     i=MIN(12,sNx)
562     j=MIN(11,sNy)
563     WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv: yFx=', af(i,j+1),
564     & localTij(i,j), vTrans(i,j+1), af(i,j+1)/vTrans(i,j+1)
565     ENDIF
566     ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
567     CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .TRUE.,
568     I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
569     O af, myThid )
570     ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
571     CALL GAD_DST3_ADV_Y( bi,bj,k, .TRUE.,
572     I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
573     O af, myThid )
574     ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
575     CALL GAD_DST3FL_ADV_Y( bi,bj,k, .TRUE.,
576     I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
577     O af, myThid )
578     #ifndef ALLOW_AUTODIFF_TAMC
579     ELSEIF (advectionScheme.EQ.ENUM_OS7MP ) THEN
580     CALL GAD_OS7MP_ADV_Y( bi,bj,k, .TRUE.,
581     I SEAICE_deltaTtherm, vTrans, vFld, maskLocS, localTij,
582     O af, myThid )
583     #endif
584     ELSE
585     STOP
586     & 'SEAICE_ADVECTION: adv. scheme incompatibale with mutli-dim'
587     ENDIF
588    
589     C- Advective flux in Y : done
590     ENDIF
591    
592     cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
593     C- Internal exchange for next calculations in X
594     IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
595     CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
596     & localTij, bi,bj, myThid )
597     ENDIF
598     cph-exch2#endif
599    
600     C- Update the local seaice field where needed:
601    
602     C update in overlap-Only
603     IF ( overlapOnly ) THEN
604     jMinUpd = 1-OLy+1
605     jMaxUpd = sNy+OLy-1
606     C- notes: these 2 lines below have no real effect (because recip_hFac=0
607     C in corner region) but safer to keep them.
608     IF ( S_edge ) jMinUpd = 1
609     IF ( N_edge ) jMaxUpd = sNy
610    
611     IF ( W_edge .AND. extensiveFld ) THEN
612     DO j=jMinUpd,jMaxUpd
613     DO i=1-OLx,0
614     localTij(i,j)=localTij(i,j)
615     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
616     & *recip_rA(i,j,bi,bj)
617     & *( af(i,j+1)-af(i,j)
618     & )
619     ENDDO
620     ENDDO
621     ELSEIF ( W_edge ) THEN
622     DO j=jMinUpd,jMaxUpd
623     DO i=1-OLx,0
624     localTij(i,j)=localTij(i,j)
625     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
626     & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
627     & *( (af(i,j+1)-af(i,j))
628     & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
629     & )
630     ENDDO
631     ENDDO
632     ENDIF
633     IF ( E_edge .AND. extensiveFld ) THEN
634     DO j=jMinUpd,jMaxUpd
635     DO i=sNx+1,sNx+OLx
636     localTij(i,j)=localTij(i,j)
637     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
638     & *recip_rA(i,j,bi,bj)
639     & *( af(i,j+1)-af(i,j)
640     & )
641     ENDDO
642     ENDDO
643     ELSEIF ( E_edge ) THEN
644     DO j=jMinUpd,jMaxUpd
645     DO i=sNx+1,sNx+OLx
646     localTij(i,j)=localTij(i,j)
647     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
648     & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
649     & *( (af(i,j+1)-af(i,j))
650     & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
651     & )
652     ENDDO
653     ENDDO
654     ENDIF
655     C-- keep advective flux (for diagnostics)
656     IF ( W_edge ) THEN
657     DO j=1-OLy+1,sNy+OLy
658     DO i=1-OLx,0
659     afy(i,j) = af(i,j)
660     ENDDO
661     ENDDO
662     ENDIF
663     IF ( E_edge ) THEN
664     DO j=1-OLy+1,sNy+OLy
665     DO i=sNx+1,sNx+OLx
666     afy(i,j) = af(i,j)
667     ENDDO
668     ENDDO
669     ENDIF
670    
671     ELSE
672     C do not only update the overlap
673     iMinUpd = 1-OLx
674     iMaxUpd = sNx+OLx
675     IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
676     IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
677     IF ( extensiveFld ) THEN
678     DO j=1-OLy+1,sNy+OLy-1
679     DO i=iMinUpd,iMaxUpd
680     localTij(i,j)=localTij(i,j)
681     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
682     & *recip_rA(i,j,bi,bj)
683     & *( af(i,j+1)-af(i,j)
684     & )
685     ENDDO
686     ENDDO
687     ELSE
688     DO j=1-OLy+1,sNy+OLy-1
689     DO i=iMinUpd,iMaxUpd
690     localTij(i,j)=localTij(i,j)
691     & -SEAICE_deltaTtherm*maskInC(i,j,bi,bj)
692     & *recip_rA(i,j,bi,bj)*r_hFld(i,j)
693     & *( (af(i,j+1)-af(i,j))
694     & -(vTrans(i,j+1)-vTrans(i,j))*iceFld(i,j)
695     & )
696     ENDDO
697     ENDDO
698     ENDIF
699     C-- keep advective flux (for diagnostics)
700     DO j=1-OLy+1,sNy+OLy
701     DO i=iMinUpd,iMaxUpd
702     afy(i,j) = af(i,j)
703     ENDDO
704     ENDDO
705    
706     C end if/else update overlap-Only
707     ENDIF
708    
709     C-- End of Y direction
710     ENDIF
711    
712     C-- End of ipass loop
713     ENDDO
714    
715     C- explicit advection is done ; store tendency in gFld:
716     DO j=1-OLy,sNy+OLy
717     DO i=1-OLx,sNx+OLx
718     gFld(i,j)=(localTij(i,j)-iceFld(i,j))/SEAICE_deltaTtherm
719     ENDDO
720     ENDDO
721     IF ( dBug .AND. bi.EQ.3 ) THEN
722     i=MIN(12,sNx)
723     j=MIN(11,sNy)
724     tmpFac= SEAICE_deltaTtherm*recip_rA(i,j,bi,bj)
725     WRITE(ioUnit,'(A,1P4E14.6)') 'ICE_adv:',
726     & afx(i,j)*tmpFac,afx(i+1,j)*tmpFac,
727     & afy(i,j)*tmpFac,afy(i,j+1)*tmpFac
728     ENDIF
729    
730     #ifdef ALLOW_DIAGNOSTICS
731     IF ( useDiagnostics ) THEN
732     diagName = 'ADVx'//diagSufx
733     CALL DIAGNOSTICS_FILL(afx,diagName, k,1, 2,bi,bj, myThid)
734     diagName = 'ADVy'//diagSufx
735     CALL DIAGNOSTICS_FILL(afy,diagName, k,1, 2,bi,bj, myThid)
736     ENDIF
737     #endif
738    
739     #ifdef ALLOW_DEBUG
740     IF ( debugLevel .GE. debLevC
741     & .AND. tracerIdentity.EQ.GAD_HEFF
742     & .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
743     & .AND. nPx.EQ.1 .AND. nPy.EQ.1
744     & .AND. useCubedSphereExchange ) THEN
745     CALL DEBUG_CS_CORNER_UV( ' afx,afy from SEAICE_ADVECTION',
746     & afx,afy, k, standardMessageUnit,bi,bj,myThid )
747     ENDIF
748     #endif /* ALLOW_DEBUG */
749    
750     RETURN
751     END

  ViewVC Help
Powered by ViewVC 1.1.22