| 1 | jmc | 1.7 | C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_som_adv_y.F,v 1.6 2012/03/05 17:59:15 jmc Exp $ | 
| 2 | jmc | 1.1 | C $Name:  $ | 
| 3 |  |  |  | 
| 4 |  |  | #include "GAD_OPTIONS.h" | 
| 5 |  |  |  | 
| 6 |  |  | CBOP | 
| 7 |  |  | C !ROUTINE: GAD_SOM_ADV_Y | 
| 8 |  |  |  | 
| 9 |  |  | C !INTERFACE: ========================================================== | 
| 10 |  |  | SUBROUTINE GAD_SOM_ADV_Y( | 
| 11 |  |  | I           bi,bj,k, limiter, | 
| 12 | jmc | 1.4 | I           overlapOnly, interiorOnly, | 
| 13 |  |  | I           N_edge, S_edge, E_edge, W_edge, | 
| 14 | jmc | 1.6 | I           deltaTloc, vTrans, maskIn, | 
| 15 | jmc | 1.1 | 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, | 
| 17 |  |  | O           vT, | 
| 18 |  |  | I           myThid ) | 
| 19 |  |  |  | 
| 20 |  |  | C !DESCRIPTION: | 
| 21 |  |  | C  Calculates the area integrated meridional flux due to advection | 
| 22 |  |  | C  of a tracer using | 
| 23 |  |  | C-- | 
| 24 |  |  | C        Second-Order Moments Advection of tracer in Y-direction | 
| 25 |  |  | C        ref: M.J.Prather, 1986, JGR, 91, D6, pp 6671-6681. | 
| 26 |  |  | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 27 |  |  | C      The 3-D grid has dimension  (Nx,Ny,Nz) with corresponding | 
| 28 |  |  | C      velocity field (U,V,W).  Parallel subroutine calculate | 
| 29 |  |  | C      advection in the X- and Z- directions. | 
| 30 |  |  | C      The moment [Si] are as defined in the text, Sm refers to | 
| 31 |  |  | C      the total mass in each grid box | 
| 32 |  |  | C      the moments [Fi] are similarly defined and used as temporary | 
| 33 |  |  | C      storage for portions of the grid boxes in transit. | 
| 34 |  |  | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 35 |  |  |  | 
| 36 |  |  | C !USES: =============================================================== | 
| 37 |  |  | IMPLICIT NONE | 
| 38 |  |  | #include "SIZE.h" | 
| 39 | jmc | 1.7 | #include "EEPARAMS.h" | 
| 40 | jmc | 1.1 | #include "GAD.h" | 
| 41 |  |  |  | 
| 42 |  |  | C !INPUT PARAMETERS: =================================================== | 
| 43 | jmc | 1.4 | C  bi,bj         :: tile indices | 
| 44 |  |  | C  k             :: vertical level | 
| 45 |  |  | C  limiter       :: 0: no limiter ; 1: Prather, 1986 limiter | 
| 46 |  |  | C  overlapOnly   :: only update the edges of myTile, but not the interior | 
| 47 |  |  | 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 | 
| 49 |  |  | C  vTrans        :: zonal volume transport | 
| 50 | jmc | 1.6 | C  maskIn        :: 2-D array Interior mask | 
| 51 | jmc | 1.4 | C  myThid        :: my Thread Id. number | 
| 52 | jmc | 1.1 | INTEGER bi,bj,k | 
| 53 |  |  | INTEGER limiter | 
| 54 | jmc | 1.4 | LOGICAL overlapOnly, interiorOnly | 
| 55 |  |  | LOGICAL N_edge, S_edge, E_edge, W_edge | 
| 56 | jmc | 1.1 | _RL deltaTloc | 
| 57 |  |  | _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 58 | jmc | 1.6 | _RS maskIn(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 59 | jmc | 1.1 | INTEGER myThid | 
| 60 |  |  |  | 
| 61 |  |  | C !OUTPUT PARAMETERS: ================================================== | 
| 62 |  |  | C  sm_v         :: volume of grid cell | 
| 63 |  |  | C  sm_o         :: tracer content of grid cell (zero order moment) | 
| 64 |  |  | C  sm_x,y,z     :: 1rst order moment of tracer distribution, in x,y,z direction | 
| 65 |  |  | C  sm_xx,yy,zz  ::  2nd order moment of tracer distribution, in x,y,z direction | 
| 66 |  |  | C  sm_xy,xz,yz  ::  2nd order moment of tracer distr., in cross direction xy,xz,yz | 
| 67 |  |  | C  vT           :: meridional advective flux | 
| 68 |  |  | _RL sm_v  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 69 |  |  | _RL sm_o  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 70 |  |  | _RL sm_x  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 71 |  |  | _RL sm_y  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 72 |  |  | _RL sm_z  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 73 |  |  | _RL sm_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 74 |  |  | _RL sm_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 75 |  |  | _RL sm_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 76 |  |  | _RL sm_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 77 |  |  | _RL sm_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 78 |  |  | _RL sm_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 79 |  |  | _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 80 |  |  |  | 
| 81 |  |  | C !LOCAL VARIABLES: ==================================================== | 
| 82 | jmc | 1.4 | C  i,j           :: loop indices | 
| 83 |  |  | C  vLoc          :: volume transported (per time step) | 
| 84 |  |  | C [iMin,iMax]Upd :: loop range to update tracer field | 
| 85 |  |  | C [jMin,jMax]Upd :: loop range to update tracer field | 
| 86 |  |  | C  nbStrips      :: number of strips (if region to update is splitted) | 
| 87 | jmc | 1.7 | _RL three | 
| 88 | jmc | 1.1 | PARAMETER( three = 3. _d 0 ) | 
| 89 |  |  | INTEGER i,j | 
| 90 | jmc | 1.4 | INTEGER ns, nbStrips | 
| 91 |  |  | INTEGER iMinUpd(2), iMaxUpd(2), jMinUpd(2), jMaxUpd(2) | 
| 92 | jmc | 1.3 | _RL  recip_dT | 
| 93 | jmc | 1.1 | _RL  slpmax, s1max, s1new, s2new | 
| 94 |  |  | _RL  vLoc, alf1, alf1q, alpmn | 
| 95 |  |  | _RL  alfp, alpq, alp1, locTp | 
| 96 |  |  | _RL  alfn, alnq, aln1, locTn | 
| 97 |  |  | _RL  alp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 98 |  |  | _RL  aln  (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 99 |  |  | _RL  fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 100 |  |  | _RL  fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 101 |  |  | _RL  fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 102 |  |  | _RL  fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 103 |  |  | _RL  fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 104 |  |  | _RL  fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 105 |  |  | _RL  fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 106 |  |  | _RL  fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 107 |  |  | _RL  fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 108 |  |  | _RL  fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 109 |  |  | _RL  fp_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 110 |  |  | _RL  fn_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 111 |  |  | _RL  fp_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 112 |  |  | _RL  fn_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 113 |  |  | _RL  fp_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 114 |  |  | _RL  fn_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 115 |  |  | _RL  fp_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 116 |  |  | _RL  fn_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 117 |  |  | _RL  fp_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 118 |  |  | _RL  fn_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 119 |  |  | _RL  fp_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 120 |  |  | _RL  fn_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy) | 
| 121 |  |  | CEOP | 
| 122 |  |  |  | 
| 123 | jmc | 1.3 | recip_dT = 0. | 
| 124 | jmc | 1.7 | IF ( deltaTloc.GT.zeroRL ) recip_dT = 1.0 _d 0 / deltaTloc | 
| 125 | jmc | 1.3 |  | 
| 126 | jmc | 1.4 | C-    Set loop ranges for updating tracer field (splitted in 2 strips) | 
| 127 |  |  | nbStrips   = 1 | 
| 128 | jmc | 1.6 | iMinUpd(1) = 1-OLx | 
| 129 |  |  | iMaxUpd(1) = sNx+OLx | 
| 130 |  |  | jMinUpd(1) = 1-OLy+1 | 
| 131 |  |  | jMaxUpd(1) = sNy+OLy-1 | 
| 132 | jmc | 1.4 | 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 | jmc | 1.6 | iMinUpd(1) = 1-OLx | 
| 138 | jmc | 1.4 | iMaxUpd(1) = 0 | 
| 139 |  |  | ENDIF | 
| 140 |  |  | IF ( E_edge ) THEN | 
| 141 |  |  | IF ( W_edge ) nbStrips = 2 | 
| 142 |  |  | iMinUpd(nbStrips) = sNx+1 | 
| 143 | jmc | 1.6 | iMaxUpd(nbStrips) = sNx+OLx | 
| 144 | jmc | 1.4 | 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 | jmc | 1.1 | IF ( limiter.EQ.1 ) THEN | 
| 155 | jmc | 1.4 | DO j=jMinUpd(1)-1,jMaxUpd(1)+1 | 
| 156 |  |  | DO i=iMinUpd(ns),iMaxUpd(ns) | 
| 157 | jmc | 1.1 | C     If flux-limiting transport is to be applied, place limits on | 
| 158 |  |  | C     appropriate moments before transport. | 
| 159 |  |  | slpmax = 0. | 
| 160 |  |  | IF ( sm_o(i,j).GT.0. ) slpmax = sm_o(i,j) | 
| 161 |  |  | s1max = slpmax*1.5 _d 0 | 
| 162 |  |  | s1new = MIN(  s1max, MAX(-s1max,sm_y(i,j)) ) | 
| 163 |  |  | s2new = MIN( (slpmax+slpmax-ABS(s1new)/three), | 
| 164 |  |  | &                MAX(ABS(s1new)-slpmax,sm_yy(i,j))  ) | 
| 165 |  |  | sm_xy(i,j) = MIN( slpmax, MAX(-slpmax,sm_xy(i,j)) ) | 
| 166 |  |  | sm_yz(i,j) = MIN( slpmax, MAX(-slpmax,sm_yz(i,j)) ) | 
| 167 | jmc | 1.2 | sm_y (i,j) = s1new | 
| 168 |  |  | sm_yy(i,j) = s2new | 
| 169 | jmc | 1.1 | ENDDO | 
| 170 |  |  | ENDDO | 
| 171 |  |  | ENDIF | 
| 172 |  |  |  | 
| 173 |  |  | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 174 |  |  | C---  part.1 : calculate flux for all moments | 
| 175 | jmc | 1.4 | DO j=jMinUpd(1),jMaxUpd(1)+1 | 
| 176 |  |  | DO i=iMinUpd(ns),iMaxUpd(ns) | 
| 177 | jmc | 1.1 | 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) | 
| 179 | jmc | 1.7 | fp_v (i,j) = MAX( zeroRL,  vLoc ) | 
| 180 | jmc | 1.1 | alp  (i,j) = fp_v(i,j)/sm_v(i,j-1) | 
| 181 |  |  | alpq       = alp(i,j)*alp(i,j) | 
| 182 |  |  | alp1       = 1. _d 0 - alp(i,j) | 
| 183 |  |  | C-     Create temporary moments/masses for partial boxes in transit | 
| 184 |  |  | C       use same indexing as velocity, "p" for positive V | 
| 185 |  |  | fp_o (i,j) = alp(i,j)*( sm_o(i,j-1) + alp1*sm_y(i,j-1) | 
| 186 |  |  | &                        + alp1*(alp1-alp(i,j))*sm_yy(i,j-1) | 
| 187 |  |  | &                        ) | 
| 188 |  |  | fp_y (i,j) = alpq    *( sm_y(i,j-1) + three*alp1*sm_yy(i,j-1) ) | 
| 189 |  |  | fp_yy(i,j) = alp(i,j)*alpq*sm_yy(i,j-1) | 
| 190 |  |  | fp_x (i,j) = alp(i,j)*( sm_x(i,j-1) + alp1*sm_xy(i,j-1) ) | 
| 191 |  |  | fp_z (i,j) = alp(i,j)*( sm_z(i,j-1) + alp1*sm_yz(i,j-1) ) | 
| 192 |  |  |  | 
| 193 |  |  | fp_xy(i,j) = alpq    *sm_xy(i,j-1) | 
| 194 |  |  | fp_yz(i,j) = alpq    *sm_yz(i,j-1) | 
| 195 |  |  | fp_xx(i,j) = alp(i,j)*sm_xx(i,j-1) | 
| 196 |  |  | 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) | 
| 198 |  |  | C--    Flux from (j) to (j-1) when V<0 (i.e., take left side of box j) | 
| 199 | jmc | 1.7 | fn_v (i,j) = MAX( zeroRL, -vLoc ) | 
| 200 | jmc | 1.1 | aln  (i,j) = fn_v(i,j)/sm_v(i, j ) | 
| 201 |  |  | alnq       = aln(i,j)*aln(i,j) | 
| 202 |  |  | aln1       = 1. _d 0 - aln(i,j) | 
| 203 |  |  | C-     Create temporary moments/masses for partial boxes in transit | 
| 204 |  |  | C       use same indexing as velocity, "n" for negative V | 
| 205 |  |  | fn_o (i,j) = aln(i,j)*( sm_o(i, j ) - aln1*sm_y(i, j ) | 
| 206 |  |  | &                        + aln1*(aln1-aln(i,j))*sm_yy(i, j ) | 
| 207 |  |  | &                        ) | 
| 208 |  |  | fn_y (i,j) = alnq    *( sm_y(i, j ) - three*aln1*sm_yy(i, j ) ) | 
| 209 |  |  | fn_yy(i,j) = aln(i,j)*alnq*sm_yy(i, j ) | 
| 210 |  |  | fn_x (i,j) = aln(i,j)*( sm_x(i, j ) - aln1*sm_xy(i, j ) ) | 
| 211 |  |  | fn_z (i,j) = aln(i,j)*( sm_z(i, j ) - aln1*sm_yz(i, j ) ) | 
| 212 |  |  | fn_xy(i,j) = alnq    *sm_xy(i, j ) | 
| 213 |  |  | fn_yz(i,j) = alnq    *sm_yz(i, j ) | 
| 214 |  |  | fn_xx(i,j) = aln(i,j)*sm_xx(i, j ) | 
| 215 |  |  | fn_zz(i,j) = aln(i,j)*sm_zz(i, j ) | 
| 216 |  |  | fn_xz(i,j) = aln(i,j)*sm_xz(i, j ) | 
| 217 |  |  | C--    Save zero-order flux: | 
| 218 | jmc | 1.3 | vT(i,j) = ( fp_o(i,j) - fn_o(i,j) )*recip_dT | 
| 219 | jmc | 1.1 | ENDDO | 
| 220 |  |  | ENDDO | 
| 221 |  |  |  | 
| 222 | jmc | 1.4 | C--   end 1rst loop on strip number "ns" | 
| 223 |  |  | c     ENDDO | 
| 224 |  |  |  | 
| 225 | jmc | 1.1 | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 226 | jmc | 1.4 | C--   start 2nd loop on strip number "ns" | 
| 227 |  |  | c     DO ns=1,nbStrips | 
| 228 |  |  |  | 
| 229 | jmc | 1.1 | 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) | 
| 231 | jmc | 1.4 | DO j=jMinUpd(1),jMaxUpd(1) | 
| 232 |  |  | DO i=iMinUpd(ns),iMaxUpd(ns) | 
| 233 | jmc | 1.6 | #ifdef ALLOW_OBCS | 
| 234 | jmc | 1.7 | IF ( maskIn(i,j).NE.zeroRS ) THEN | 
| 235 | jmc | 1.6 | #endif /* ALLOW_OBCS */ | 
| 236 | jmc | 1.1 | alf1  = 1. _d 0 - aln(i,j) - alp(i,j+1) | 
| 237 |  |  | alf1q = alf1*alf1 | 
| 238 |  |  | alpmn = alp(i,j+1) - aln(i,j) | 
| 239 |  |  | sm_v (i,j) = sm_v (i,j) - fn_v (i,j) - fp_v (i,j+1) | 
| 240 |  |  | sm_o (i,j) = sm_o (i,j) - fn_o (i,j) - fp_o (i,j+1) | 
| 241 |  |  | sm_y (i,j) = alf1q*( sm_y(i,j) - three*alpmn*sm_yy(i,j) ) | 
| 242 |  |  | sm_yy(i,j) = alf1*alf1q*sm_yy(i,j) | 
| 243 |  |  | sm_xy(i,j) = alf1q*sm_xy(i,j) | 
| 244 |  |  | sm_yz(i,j) = alf1q*sm_yz(i,j) | 
| 245 |  |  | sm_x (i,j) = sm_x (i,j) - fn_x (i,j) - fp_x (i,j+1) | 
| 246 |  |  | sm_xx(i,j) = sm_xx(i,j) - fn_xx(i,j) - fp_xx(i,j+1) | 
| 247 |  |  | 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) | 
| 249 |  |  | sm_xz(i,j) = sm_xz(i,j) - fn_xz(i,j) - fp_xz(i,j+1) | 
| 250 | jmc | 1.6 | #ifdef ALLOW_OBCS | 
| 251 |  |  | ENDIF | 
| 252 |  |  | #endif /* ALLOW_OBCS */ | 
| 253 | jmc | 1.1 | ENDDO | 
| 254 |  |  | ENDDO | 
| 255 |  |  |  | 
| 256 |  |  | C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| | 
| 257 |  |  | 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) | 
| 259 | jmc | 1.4 | DO j=jMinUpd(1),jMaxUpd(1) | 
| 260 |  |  | DO i=iMinUpd(ns),iMaxUpd(ns) | 
| 261 | jmc | 1.6 | #ifdef ALLOW_OBCS | 
| 262 | jmc | 1.7 | IF ( maskIn(i,j).NE.zeroRS ) THEN | 
| 263 | jmc | 1.6 | #endif /* ALLOW_OBCS */ | 
| 264 | jmc | 1.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) | 
| 266 |  |  | alfn = fn_v(i,j+1)/sm_v(i,j) | 
| 267 |  |  | alf1 = 1. _d 0 - alfp - alfn | 
| 268 |  |  | alp1 = 1. _d 0 - alfp | 
| 269 |  |  | aln1 = 1. _d 0 - alfn | 
| 270 |  |  | alpmn = alfp - alfn | 
| 271 |  |  | locTp = alfp*sm_o(i,j) - alp1*fp_o(i,j) | 
| 272 |  |  | locTn = alfn*sm_o(i,j) - aln1*fn_o(i,j+1) | 
| 273 |  |  | sm_yy(i,j) = alf1*alf1*sm_yy(i,j) + alfp*alfp*fp_yy(i,j) | 
| 274 |  |  | &                                    + alfn*alfn*fn_yy(i,j+1) | 
| 275 |  |  | &   - 5. _d 0*(-alpmn*alf1*sm_y(i,j) + alfp*alp1*fp_y(i,j) | 
| 276 |  |  | &                                    - alfn*aln1*fn_y(i,j+1) | 
| 277 | jmc | 1.7 | &             + twoRL*alfp*alfn*sm_o(i,j) + (alp1-alfp)*locTp | 
| 278 |  |  | &                                         + (aln1-alfn)*locTn | 
| 279 | jmc | 1.1 | &             ) | 
| 280 |  |  | sm_xy(i,j) = alf1*sm_xy(i,j) + alfp*fp_xy(i,j) | 
| 281 |  |  | &                               + alfn*fn_xy(i,j+1) | 
| 282 |  |  | &     + three*( alpmn*sm_x(i,j) - alp1*fp_x(i,j) | 
| 283 |  |  | &                               + aln1*fn_x(i,j+1) | 
| 284 |  |  | &             ) | 
| 285 |  |  | sm_yz(i,j) = alf1*sm_yz(i,j) + alfp*fp_yz(i,j) | 
| 286 |  |  | &                               + alfn*fn_yz(i,j+1) | 
| 287 |  |  | &     + three*( alpmn*sm_z(i,j) - alp1*fp_z(i,j) | 
| 288 |  |  | &                               + aln1*fn_z(i,j+1) | 
| 289 |  |  | &             ) | 
| 290 |  |  | sm_y (i,j) = alf1*sm_y(i,j) + alfp*fp_y(i,j) + alfn*fn_y(i,j+1) | 
| 291 |  |  | &             + three*( locTp - locTn ) | 
| 292 |  |  | sm_o (i,j) = sm_o (i,j) + fp_o (i,j) + fn_o (i,j+1) | 
| 293 |  |  | sm_x (i,j) = sm_x (i,j) + fp_x (i,j) + fn_x (i,j+1) | 
| 294 |  |  | sm_xx(i,j) = sm_xx(i,j) + fp_xx(i,j) + fn_xx(i,j+1) | 
| 295 |  |  | 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) | 
| 297 |  |  | sm_xz(i,j) = sm_xz(i,j) + fp_xz(i,j) + fn_xz(i,j+1) | 
| 298 | jmc | 1.6 | #ifdef ALLOW_OBCS | 
| 299 |  |  | ENDIF | 
| 300 |  |  | #endif /* ALLOW_OBCS */ | 
| 301 | jmc | 1.1 | ENDDO | 
| 302 |  |  | ENDDO | 
| 303 |  |  |  | 
| 304 | jmc | 1.4 | C--   end 2nd loop on strip number "ns" | 
| 305 |  |  | ENDDO | 
| 306 |  |  |  | 
| 307 | jmc | 1.1 | RETURN | 
| 308 |  |  | END |