/[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.6 by jmc, Mon Mar 5 17:59:15 2012 UTC revision 1.7 by jmc, Sat Mar 2 00:29:20 2013 UTC
# 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 83  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 121  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
# Line 176  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 196  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 231  C      take off from grid box (j): negat Line 231  C      take off from grid box (j): negat
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  #ifdef ALLOW_OBCS
234          IF ( maskIn(i,j).NE.0. ) THEN          IF ( maskIn(i,j).NE.zeroRS ) THEN
235  #endif /* ALLOW_OBCS */  #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
# Line 259  C      add into grid box (j): positive V Line 259  C      add into grid box (j): positive V
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  #ifdef ALLOW_OBCS
262          IF ( maskIn(i,j).NE.0. ) THEN          IF ( maskIn(i,j).NE.zeroRS ) THEN
263  #endif /* ALLOW_OBCS */  #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)
# Line 274  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)
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22