/[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.2 by jmc, Wed Aug 22 00:26:01 2007 UTC revision 1.6 by jmc, Mon Mar 5 17:59:15 2012 UTC
# Line 9  C !ROUTINE: GAD_SOM_ADV_Y Line 9  C !ROUTINE: GAD_SOM_ADV_Y
9  C !INTERFACE: ==========================================================  C !INTERFACE: ==========================================================
10        SUBROUTINE GAD_SOM_ADV_Y(        SUBROUTINE GAD_SOM_ADV_Y(
11       I           bi,bj,k, limiter,       I           bi,bj,k, limiter,
12       I           deltaTloc, vTrans,       I           overlapOnly, interiorOnly,
13         I           N_edge, S_edge, E_edge, W_edge,
14         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 34  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"
 c #include "GRID.h"  
39  #include "GAD.h"  #include "GAD.h"
40    
41  C !INPUT PARAMETERS: ===================================================  C !INPUT PARAMETERS: ===================================================
42  C  bi,bj        :: tile indices  C  bi,bj         :: tile indices
43  C  k            :: vertical level  C  k             :: vertical level
44  C  limiter      :: 0: no limiter ; 1: Prather, 1986 limiter  C  limiter       :: 0: no limiter ; 1: Prather, 1986 limiter
45  C  vTrans       :: zonal volume transport  C  overlapOnly   :: only update the edges of myTile, but not the interior
46  C  myThid       :: my Thread Id. number  C  interiorOnly  :: only update the interior of myTile, but not the edges
47    C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
48    C  vTrans        :: zonal volume transport
49    C  maskIn        :: 2-D array Interior mask
50    C  myThid        :: my Thread Id. number
51        INTEGER bi,bj,k        INTEGER bi,bj,k
52        INTEGER limiter        INTEGER limiter
53          LOGICAL overlapOnly, interiorOnly
54          LOGICAL N_edge, S_edge, E_edge, W_edge
55        _RL deltaTloc        _RL deltaTloc
56        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
57  c     _RL tracer(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RS maskIn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
58        INTEGER myThid        INTEGER myThid
59    
60  C !OUTPUT PARAMETERS: ==================================================  C !OUTPUT PARAMETERS: ==================================================
# Line 70  C  vT           :: meridional advective Line 77  C  vT           :: meridional advective
77        _RL sm_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL sm_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
78        _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
79    
 #ifdef GAD_ALLOW_SOM_ADVECT  
80  C !LOCAL VARIABLES: ====================================================  C !LOCAL VARIABLES: ====================================================
81  C  i,j          :: loop indices  C  i,j           :: loop indices
82  C  vLoc         :: volume transported (per time step)  C  vLoc          :: volume transported (per time step)
83    C [iMin,iMax]Upd :: loop range to update tracer field
84    C [jMin,jMax]Upd :: loop range to update tracer field
85    C  nbStrips      :: number of strips (if region to update is splitted)
86        _RL two, three        _RL two, three
87        PARAMETER( two   = 2. _d 0 )        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
91          INTEGER iMinUpd(2), iMaxUpd(2), jMinUpd(2), jMaxUpd(2)
92          _RL  recip_dT
93        _RL  slpmax, s1max, s1new, s2new        _RL  slpmax, s1max, s1new, s2new
94        _RL  vLoc, alf1, alf1q, alpmn        _RL  vLoc, alf1, alf1q, alpmn
95        _RL  alfp, alpq, alp1, locTp        _RL  alfp, alpq, alp1, locTp
# Line 108  C  vLoc         :: volume transported (p Line 120  C  vLoc         :: volume transported (p
120        _RL  fn_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL  fn_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
121  CEOP  CEOP
122    
123          recip_dT = 0.
124          IF ( deltaTloc.GT.0. _d 0 ) recip_dT = 1.0 _d 0 / deltaTloc
125    
126    C-    Set loop ranges for updating tracer field (splitted in 2 strips)
127          nbStrips   = 1
128          iMinUpd(1) = 1-OLx
129          iMaxUpd(1) = sNx+OLx
130          jMinUpd(1) = 1-OLy+1
131          jMaxUpd(1) = sNy+OLy-1
132          IF ( overlapOnly ) THEN
133    C     update in overlap-Only
134            IF ( S_edge ) jMinUpd(1) = 1
135            IF ( N_edge ) jMaxUpd(1) = sNy
136            IF ( W_edge ) THEN
137              iMinUpd(1) = 1-OLx
138              iMaxUpd(1) = 0
139            ENDIF
140            IF ( E_edge ) THEN
141              IF ( W_edge ) nbStrips = 2
142              iMinUpd(nbStrips) = sNx+1
143              iMaxUpd(nbStrips) = sNx+OLx
144            ENDIF
145          ELSE
146    C     do not only update the overlap
147            IF ( interiorOnly .AND. W_edge ) iMinUpd(1) = 1
148            IF ( interiorOnly .AND. E_edge ) iMaxUpd(1) = sNx
149          ENDIF
150    
151    C--   start 1rst loop on strip number "ns"
152          DO ns=1,nbStrips
153    
154        IF ( limiter.EQ.1 ) THEN        IF ( limiter.EQ.1 ) THEN
155         DO j=1-OLy,sNy+OLy         DO j=jMinUpd(1)-1,jMaxUpd(1)+1
156          DO i=1-OLx,sNx+OLx          DO i=iMinUpd(ns),iMaxUpd(ns)
157  C     If flux-limiting transport is to be applied, place limits on  C     If flux-limiting transport is to be applied, place limits on
158  C     appropriate moments before transport.  C     appropriate moments before transport.
159           slpmax = 0.           slpmax = 0.
# Line 129  C     appropriate moments before transpo Line 172  C     appropriate moments before transpo
172    
173  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
174  C---  part.1 : calculate flux for all moments  C---  part.1 : calculate flux for all moments
175        DO i=1-OLx,sNx+OLx        DO j=jMinUpd(1),jMaxUpd(1)+1
176         vT(i,1-OLy)=0.         DO i=iMinUpd(ns),iMaxUpd(ns)
       ENDDO  
       DO j=1-OLy+1,sNy+OLy  
        DO i=1-OLx,sNx+OLx  
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( 0. _d 0,  vLoc )
# Line 175  C       use same indexing as velocity, " Line 215  C       use same indexing as velocity, "
215          fn_zz(i,j) = aln(i,j)*sm_zz(i, j )          fn_zz(i,j) = aln(i,j)*sm_zz(i, j )
216          fn_xz(i,j) = aln(i,j)*sm_xz(i, j )          fn_xz(i,j) = aln(i,j)*sm_xz(i, j )
217  C--    Save zero-order flux:  C--    Save zero-order flux:
218          vT(i,j)    = fp_o(i,j) - fn_o(i,j)          vT(i,j) = ( fp_o(i,j) - fn_o(i,j) )*recip_dT
219         ENDDO         ENDDO
220        ENDDO        ENDDO
221    
222    C--   end 1rst loop on strip number "ns"
223    c     ENDDO
224    
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"
227    c     DO ns=1,nbStrips
228    
229  C---  part.2 : re-adjust moments remaining in the box  C---  part.2 : re-adjust moments remaining in the box
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=1-OLy+1,sNy+OLy-1        DO j=jMinUpd(1),jMaxUpd(1)
232         DO i=1-OLx,sNx+OLx         DO i=iMinUpd(ns),iMaxUpd(ns)
233    #ifdef ALLOW_OBCS
234            IF ( maskIn(i,j).NE.0. ) 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 198  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    
256  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
257  C---  part.3 : Put the temporary moments into appropriate neighboring boxes  C---  part.3 : Put the temporary moments into appropriate neighboring boxes
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=1-OLy+1,sNy+OLy-1        DO j=jMinUpd(1),jMaxUpd(1)
260         DO i=1-OLx,sNx+OLx         DO i=iMinUpd(ns),iMaxUpd(ns)
261    #ifdef ALLOW_OBCS
262            IF ( maskIn(i,j).NE.0. ) 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 240  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    
304  #endif /* GAD_ALLOW_SOM_ADVECT */  C--   end 2nd loop on strip number "ns"
305          ENDDO
306    
307        RETURN        RETURN
308        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.6
  ViewVC Help
Powered by ViewVC 1.1.22