/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_som_adv_y.F
ViewVC logotype

Diff of /MITgcm/pkg/generic_advdiff/gad_som_adv_y.F

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

revision 1.5 by jmc, Fri May 9 21:43:16 2008 UTC revision 1.7 by jmc, Sat Mar 2 00:29:20 2013 UTC
# Line 11  C !INTERFACE: ========================== Line 11  C !INTERFACE: ==========================
11       I           bi,bj,k, limiter,       I           bi,bj,k, limiter,
12       I           overlapOnly, interiorOnly,       I           overlapOnly, interiorOnly,
13       I           N_edge, S_edge, E_edge, W_edge,       I           N_edge, S_edge, E_edge, W_edge,
14       I           deltaTloc, vTrans,       I           deltaTloc, vTrans, maskIn,
15       U           sm_v, sm_o, sm_x, sm_y, sm_z,       U           sm_v, sm_o, sm_x, sm_y, sm_z,
16       U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,       U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,
17       O           vT,       O           vT,
# Line 36  C---+----1----+----2----+----3----+----4 Line 36  C---+----1----+----2----+----3----+----4
36  C !USES: ===============================================================  C !USES: ===============================================================
37        IMPLICIT NONE        IMPLICIT NONE
38  #include "SIZE.h"  #include "SIZE.h"
39    #include "EEPARAMS.h"
40  #include "GAD.h"  #include "GAD.h"
41    
42  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
# Line 46  C  overlapOnly   :: only update the edge Line 47  C  overlapOnly   :: only update the edge
47  C  interiorOnly  :: only update the interior of myTile, but not the edges  C  interiorOnly  :: only update the interior of myTile, but not the edges
48  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube  C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
49  C  vTrans        :: zonal volume transport  C  vTrans        :: zonal volume transport
50    C  maskIn        :: 2-D array Interior mask
51  C  myThid        :: my Thread Id. number  C  myThid        :: my Thread Id. number
52        INTEGER bi,bj,k        INTEGER bi,bj,k
53        INTEGER limiter        INTEGER limiter
# Line 53  C  myThid        :: my Thread Id. number Line 55  C  myThid        :: my Thread Id. number
55        LOGICAL N_edge, S_edge, E_edge, W_edge        LOGICAL N_edge, S_edge, E_edge, W_edge
56        _RL deltaTloc        _RL deltaTloc
57        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58          _RS maskIn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
59        INTEGER myThid        INTEGER myThid
60    
61  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
# Line 81  C  vLoc          :: volume transported ( Line 84  C  vLoc          :: volume transported (
84  C [iMin,iMax]Upd :: loop range to update tracer field  C [iMin,iMax]Upd :: loop range to update tracer field
85  C [jMin,jMax]Upd :: loop range to update tracer field  C [jMin,jMax]Upd :: loop range to update tracer field
86  C  nbStrips      :: number of strips (if region to update is splitted)  C  nbStrips      :: number of strips (if region to update is splitted)
87        _RL two, three        _RL three
       PARAMETER( two   = 2. _d 0 )  
88        PARAMETER( three = 3. _d 0 )        PARAMETER( three = 3. _d 0 )
89        INTEGER i,j        INTEGER i,j
90        INTEGER ns, nbStrips        INTEGER ns, nbStrips
# Line 119  C  nbStrips      :: number of strips (if Line 121  C  nbStrips      :: number of strips (if
121  CEOP  CEOP
122    
123        recip_dT = 0.        recip_dT = 0.
124        IF ( deltaTloc.GT.0. _d 0 ) recip_dT = 1.0 _d 0 / deltaTloc        IF ( deltaTloc.GT.zeroRL ) recip_dT = 1.0 _d 0 / deltaTloc
125    
126  C-    Set loop ranges for updating tracer field (splitted in 2 strips)  C-    Set loop ranges for updating tracer field (splitted in 2 strips)
127        nbStrips   = 1        nbStrips   = 1
128        iMinUpd(1) = 1-Olx        iMinUpd(1) = 1-OLx
129        iMaxUpd(1) = sNx+Olx        iMaxUpd(1) = sNx+OLx
130        jMinUpd(1) = 1-Oly+1        jMinUpd(1) = 1-OLy+1
131        jMaxUpd(1) = sNy+Oly-1        jMaxUpd(1) = sNy+OLy-1
132        IF ( overlapOnly ) THEN        IF ( overlapOnly ) THEN
133  C     update in overlap-Only  C     update in overlap-Only
134          IF ( S_edge ) jMinUpd(1) = 1          IF ( S_edge ) jMinUpd(1) = 1
135          IF ( N_edge ) jMaxUpd(1) = sNy          IF ( N_edge ) jMaxUpd(1) = sNy
136          IF ( W_edge ) THEN          IF ( W_edge ) THEN
137            iMinUpd(1) = 1-Olx            iMinUpd(1) = 1-OLx
138            iMaxUpd(1) = 0            iMaxUpd(1) = 0
139          ENDIF          ENDIF
140          IF ( E_edge ) THEN          IF ( E_edge ) THEN
141            IF ( W_edge ) nbStrips = 2            IF ( W_edge ) nbStrips = 2
142            iMinUpd(nbStrips) = sNx+1            iMinUpd(nbStrips) = sNx+1
143            iMaxUpd(nbStrips) = sNx+Olx            iMaxUpd(nbStrips) = sNx+OLx
144          ENDIF          ENDIF
145        ELSE        ELSE
146  C     do not only update the overlap  C     do not only update the overlap
# Line 146  C     do not only update the overlap Line 148  C     do not only update the overlap
148          IF ( interiorOnly .AND. E_edge ) iMaxUpd(1) = sNx          IF ( interiorOnly .AND. E_edge ) iMaxUpd(1) = sNx
149        ENDIF        ENDIF
150    
 C-    Internal exchange for calculations in Y  
 c     IF ( overlapOnly ) THEN  
 c         CALL GAD_SOM_FILL_CS_CORNER( .FALSE.,  
 c    U           sm_v,  sm_o,  sm_x,  sm_y,  sm_z,  
 c    U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,  
 c    I           bi, bj, myThid )  
 c     ENDIF  
   
151  C--   start 1rst loop on strip number "ns"  C--   start 1rst loop on strip number "ns"
152        DO ns=1,nbStrips        DO ns=1,nbStrips
153    
# Line 182  C---  part.1 : calculate flux for all mo Line 176  C---  part.1 : calculate flux for all mo
176         DO i=iMinUpd(ns),iMaxUpd(ns)         DO i=iMinUpd(ns),iMaxUpd(ns)
177          vLoc = vTrans(i,j)*deltaTloc          vLoc = vTrans(i,j)*deltaTloc
178  C--    Flux from (j-1) to (j) when V>0 (i.e., take right side of box j-1)  C--    Flux from (j-1) to (j) when V>0 (i.e., take right side of box j-1)
179          fp_v (i,j) = MAX( 0. _d 0,  vLoc )          fp_v (i,j) = MAX( zeroRL,  vLoc )
180          alp  (i,j) = fp_v(i,j)/sm_v(i,j-1)          alp  (i,j) = fp_v(i,j)/sm_v(i,j-1)
181          alpq       = alp(i,j)*alp(i,j)          alpq       = alp(i,j)*alp(i,j)
182          alp1       = 1. _d 0 - alp(i,j)          alp1       = 1. _d 0 - alp(i,j)
# Line 202  C       use same indexing as velocity, " Line 196  C       use same indexing as velocity, "
196          fp_zz(i,j) = alp(i,j)*sm_zz(i,j-1)          fp_zz(i,j) = alp(i,j)*sm_zz(i,j-1)
197          fp_xz(i,j) = alp(i,j)*sm_xz(i,j-1)          fp_xz(i,j) = alp(i,j)*sm_xz(i,j-1)
198  C--    Flux from (j) to (j-1) when V<0 (i.e., take left side of box j)  C--    Flux from (j) to (j-1) when V<0 (i.e., take left side of box j)
199          fn_v (i,j) = MAX( 0. _d 0, -vLoc )          fn_v (i,j) = MAX( zeroRL, -vLoc )
200          aln  (i,j) = fn_v(i,j)/sm_v(i, j )          aln  (i,j) = fn_v(i,j)/sm_v(i, j )
201          alnq       = aln(i,j)*aln(i,j)          alnq       = aln(i,j)*aln(i,j)
202          aln1       = 1. _d 0 - aln(i,j)          aln1       = 1. _d 0 - aln(i,j)
# Line 228  C--    Save zero-order flux: Line 222  C--    Save zero-order flux:
222  C--   end 1rst loop on strip number "ns"  C--   end 1rst loop on strip number "ns"
223  c     ENDDO  c     ENDDO
224    
 C-    Internal exchange for next calculations in X  
 c     IF ( overlapOnly ) THEN  
 c         CALL GAD_SOM_FILL_CS_CORNER( .TRUE.,  
 c    U           sm_v,  sm_o,  sm_x,  sm_y,  sm_z,  
 c    U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,  
 c    I           bi, bj, myThid )  
 c     ENDIF  
   
225  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
226  C--   start 2nd loop on strip number "ns"  C--   start 2nd loop on strip number "ns"
227  c     DO ns=1,nbStrips  c     DO ns=1,nbStrips
# Line 244  C---  part.2 : re-adjust moments remaini Line 230  C---  part.2 : re-adjust moments remaini
230  C      take off from grid box (j): negative V(j) and positive V(j+1)  C      take off from grid box (j): negative V(j) and positive V(j+1)
231        DO j=jMinUpd(1),jMaxUpd(1)        DO j=jMinUpd(1),jMaxUpd(1)
232         DO i=iMinUpd(ns),iMaxUpd(ns)         DO i=iMinUpd(ns),iMaxUpd(ns)
233    #ifdef ALLOW_OBCS
234            IF ( maskIn(i,j).NE.zeroRS ) THEN
235    #endif /* ALLOW_OBCS */
236          alf1  = 1. _d 0 - aln(i,j) - alp(i,j+1)          alf1  = 1. _d 0 - aln(i,j) - alp(i,j+1)
237          alf1q = alf1*alf1          alf1q = alf1*alf1
238          alpmn = alp(i,j+1) - aln(i,j)          alpmn = alp(i,j+1) - aln(i,j)
# Line 258  C      take off from grid box (j): negat Line 247  C      take off from grid box (j): negat
247          sm_z (i,j) = sm_z (i,j) - fn_z (i,j) - fp_z (i,j+1)          sm_z (i,j) = sm_z (i,j) - fn_z (i,j) - fp_z (i,j+1)
248          sm_zz(i,j) = sm_zz(i,j) - fn_zz(i,j) - fp_zz(i,j+1)          sm_zz(i,j) = sm_zz(i,j) - fn_zz(i,j) - fp_zz(i,j+1)
249          sm_xz(i,j) = sm_xz(i,j) - fn_xz(i,j) - fp_xz(i,j+1)          sm_xz(i,j) = sm_xz(i,j) - fn_xz(i,j) - fp_xz(i,j+1)
250    #ifdef ALLOW_OBCS
251            ENDIF
252    #endif /* ALLOW_OBCS */
253         ENDDO         ENDDO
254        ENDDO        ENDDO
255    
# Line 266  C---  part.3 : Put the temporary moments Line 258  C---  part.3 : Put the temporary moments
258  C      add into grid box (j): positive V(j) and negative V(j+1)  C      add into grid box (j): positive V(j) and negative V(j+1)
259        DO j=jMinUpd(1),jMaxUpd(1)        DO j=jMinUpd(1),jMaxUpd(1)
260         DO i=iMinUpd(ns),iMaxUpd(ns)         DO i=iMinUpd(ns),iMaxUpd(ns)
261    #ifdef ALLOW_OBCS
262            IF ( maskIn(i,j).NE.zeroRS ) THEN
263    #endif /* ALLOW_OBCS */
264          sm_v (i,j) = sm_v (i,j) + fp_v (i,j) + fn_v (i,j+1)          sm_v (i,j) = sm_v (i,j) + fp_v (i,j) + fn_v (i,j+1)
265          alfp = fp_v(i, j )/sm_v(i,j)          alfp = fp_v(i, j )/sm_v(i,j)
266          alfn = fn_v(i,j+1)/sm_v(i,j)          alfn = fn_v(i,j+1)/sm_v(i,j)
# Line 279  C      add into grid box (j): positive V Line 274  C      add into grid box (j): positive V
274       &                                    + alfn*alfn*fn_yy(i,j+1)       &                                    + alfn*alfn*fn_yy(i,j+1)
275       &   - 5. _d 0*(-alpmn*alf1*sm_y(i,j) + alfp*alp1*fp_y(i,j)       &   - 5. _d 0*(-alpmn*alf1*sm_y(i,j) + alfp*alp1*fp_y(i,j)
276       &                                    - alfn*aln1*fn_y(i,j+1)       &                                    - alfn*aln1*fn_y(i,j+1)
277       &             + two*alfp*alfn*sm_o(i,j) + (alp1-alfp)*locTp       &             + twoRL*alfp*alfn*sm_o(i,j) + (alp1-alfp)*locTp
278       &                                       + (aln1-alfn)*locTn       &                                         + (aln1-alfn)*locTn
279       &             )       &             )
280          sm_xy(i,j) = alf1*sm_xy(i,j) + alfp*fp_xy(i,j)          sm_xy(i,j) = alf1*sm_xy(i,j) + alfp*fp_xy(i,j)
281       &                               + alfn*fn_xy(i,j+1)       &                               + alfn*fn_xy(i,j+1)
# Line 300  C      add into grid box (j): positive V Line 295  C      add into grid box (j): positive V
295          sm_z (i,j) = sm_z (i,j) + fp_z (i,j) + fn_z (i,j+1)          sm_z (i,j) = sm_z (i,j) + fp_z (i,j) + fn_z (i,j+1)
296          sm_zz(i,j) = sm_zz(i,j) + fp_zz(i,j) + fn_zz(i,j+1)          sm_zz(i,j) = sm_zz(i,j) + fp_zz(i,j) + fn_zz(i,j+1)
297          sm_xz(i,j) = sm_xz(i,j) + fp_xz(i,j) + fn_xz(i,j+1)          sm_xz(i,j) = sm_xz(i,j) + fp_xz(i,j) + fn_xz(i,j+1)
298    #ifdef ALLOW_OBCS
299            ENDIF
300    #endif /* ALLOW_OBCS */
301         ENDDO         ENDDO
302        ENDDO        ENDDO
303    
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22