pkg/generic_advdiff/gad_som_adv_r.F

d7ce0d34f8 Jean*0001 #include "GAD_OPTIONS.h"
7448700841 Mart*0002 #ifdef ALLOW_PTRACERS
                0003 # include "PTRACERS_OPTIONS.h"
                0004 #endif
1574069d50 Mart*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
d7ce0d34f8 Jean*0008
                0009 CBOP
                0010 C !ROUTINE: GAD_SOM_ADV_R
                0011
                0012 C !INTERFACE: ==========================================================
                0013       SUBROUTINE GAD_SOM_ADV_R(
                0014      I           bi,bj,k, kUp, kDw,
72de869c1b Jean*0015      I           deltaTloc, rTrans, maskUp, maskIn,
d7ce0d34f8 Jean*0016      U           sm_v,  sm_o,  sm_x,  sm_y,  sm_z,
                0017      U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,
                0018      U           alp,   aln,   fp_v,  fn_v,  fp_o,  fn_o,
                0019      U           fp_x,  fn_x,  fp_y,  fn_y,  fp_z,  fn_z,
                0020      U           fp_xx, fn_xx, fp_yy, fn_yy, fp_zz, fn_zz,
                0021      U           fp_xy, fn_xy, fp_xz, fn_xz, fp_yz, fn_yz,
                0022      O           wT,
                0023      I           myThid )
                0024
                0025 C !DESCRIPTION:
                0026 C  Calculates the area integrated vertical flux due to advection
                0027 C  of a tracer using
                0028 C--
                0029 C        Second-Order Moments Advection of tracer in Z-direction
                0030 C        ref: M.J.Prather, 1986, JGR, 91, D6, pp 6671-6681.
                0031 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0032 C      The 3-D grid has dimension  (Nx,Ny,Nz) with corresponding
                0033 C      velocity field (U,V,W).  Parallel subroutine calculate
                0034 C      advection in the X- and Y- directions.
                0035 C      The moment [Si] are as defined in the text, Sm refers to
                0036 C      the total mass in each grid box
                0037 C      the moments [Fi] are similarly defined and used as temporary
                0038 C      storage for portions of the grid boxes in transit.
                0039 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0040
                0041 C !USES: ===============================================================
                0042       IMPLICIT NONE
                0043 #include "SIZE.h"
                0044 #include "EEPARAMS.h"
                0045 #include "PARAMS.h"
f0a22e969e Jean*0046 #include "GRID.h"
d7ce0d34f8 Jean*0047 #include "GAD.h"
1574069d50 Mart*0048 #ifdef ALLOW_AUTODIFF_TAMC
                0049 # include "tamc.h"
                0050 #endif
d7ce0d34f8 Jean*0051
                0052 C !INPUT PARAMETERS: ===================================================
                0053 C  bi,bj        :: tile indices
                0054 C  k            :: vertical level
                0055 C  kUp          :: index into 2 1/2D array, toggles between 1 and 2
                0056 C  kDw          :: index into 2 1/2D array, toggles between 2 and 1
                0057 C  rTrans       :: vertical volume transport
                0058 C  maskUp       :: 2-D array mask for W points
72de869c1b Jean*0059 C  maskIn       :: 2-D array Interior mask
d7ce0d34f8 Jean*0060 C  myThid       :: my Thread Id. number
                0061       INTEGER bi,bj,k, kUp, kDw
                0062       _RL deltaTloc
                0063       _RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0064       _RS maskUp(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72de869c1b Jean*0065       _RS maskIn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d7ce0d34f8 Jean*0066       INTEGER myThid
                0067
                0068 C !OUTPUT PARAMETERS: ==================================================
                0069 C  sm_v         :: volume of grid cell
                0070 C  sm_o         :: tracer content of grid cell (zero order moment)
                0071 C  sm_x,y,z     :: 1rst order moment of tracer distribution, in x,y,z direction
                0072 C  sm_xx,yy,zz  ::  2nd order moment of tracer distribution, in x,y,z direction
                0073 C  sm_xy,xz,yz  ::  2nd order moment of tracer distr., in cross direction xy,xz,yz
                0074 C  wT           :: vertical advective flux
                0075       _RL sm_v  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0076       _RL sm_o  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0077       _RL sm_x  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0078       _RL sm_y  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0079       _RL sm_z  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0080       _RL sm_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0081       _RL sm_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0082       _RL sm_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0083       _RL sm_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0084       _RL sm_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0085       _RL sm_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
                0086       _RL  alp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0087       _RL  aln  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0088       _RL  fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0089       _RL  fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0090       _RL  fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0091       _RL  fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0092       _RL  fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0093       _RL  fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0094       _RL  fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0095       _RL  fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0096       _RL  fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0097       _RL  fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0098       _RL  fp_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0099       _RL  fn_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0100       _RL  fp_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0101       _RL  fn_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0102       _RL  fp_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0103       _RL  fn_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0104       _RL  fp_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0105       _RL  fn_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0106       _RL  fp_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0107       _RL  fn_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0108       _RL  fp_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0109       _RL  fn_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
                0110       _RL wT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111
7448700841 Mart*0112 #if ( defined GAD_ALLOW_TS_SOM_ADV || defined PTRACERS_ALLOW_DYN_STATE )
d7ce0d34f8 Jean*0113 C !LOCAL VARIABLES: ====================================================
                0114 C  i,j          :: loop indices
                0115 C  wLoc         :: volume transported (per time step)
0733e7e26d Jean*0116       _RL three
d7ce0d34f8 Jean*0117       PARAMETER( three = 3. _d 0 )
                0118       INTEGER i,j
                0119       INTEGER km1
0733e7e26d Jean*0120       LOGICAL noFlowAcrossSurf
9822905e7f Jean*0121       _RL  recip_dT
d7ce0d34f8 Jean*0122       _RL  wLoc, alf1, alf1q, alpmn
                0123       _RL  alfp, alpq, alp1, locTp
                0124       _RL  alfn, alnq, aln1, locTn
                0125 CEOP
                0126
64442af1fe Jean*0127 #ifdef ALLOW_AUTODIFF
1574069d50 Mart*0128 # ifdef ALLOW_AUTODIFF_TAMC
                0129 CADJ INIT somtape_r = COMMON, 1
                0130 # endif
7a28262a98 Patr*0131       alp(1,1,kDw) = alp(1,1,kDw)
                0132       aln(1,1,kDw) = aln(1,1,kDw)
                0133       fp_v(1,1,kDw) = fp_v(1,1,kDw)
                0134       fn_v(1,1,kDw) = fn_v(1,1,kDw)
                0135       fp_o(1,1,kDw) = fp_o(1,1,kDw)
                0136       fn_o(1,1,kDw) = fn_o(1,1,kDw)
                0137       fp_x(1,1,kDw) = fp_x(1,1,kDw)
                0138       fn_x(1,1,kDw) = fn_x(1,1,kDw)
                0139       fp_y(1,1,kDw) = fp_y(1,1,kDw)
                0140       fn_y(1,1,kDw) = fn_y(1,1,kDw)
                0141       fp_z(1,1,kDw) = fp_z(1,1,kDw)
                0142       fn_z(1,1,kDw) = fn_z(1,1,kDw)
                0143       fp_xx(1,1,kDw) = fp_xx(1,1,kDw)
                0144       fn_xx(1,1,kDw) = fn_xx(1,1,kDw)
                0145       fp_yy(1,1,kDw) = fp_yy(1,1,kDw)
                0146       fn_yy(1,1,kDw) = fn_yy(1,1,kDw)
                0147       fp_zz(1,1,kDw) = fp_zz(1,1,kDw)
                0148       fn_zz(1,1,kDw) = fn_zz(1,1,kDw)
                0149       fp_xy(1,1,kDw) = fp_xy(1,1,kDw)
                0150       fn_xy(1,1,kDw) = fn_xy(1,1,kDw)
                0151       fp_xz(1,1,kDw) = fp_xz(1,1,kDw)
                0152       fn_xz(1,1,kDw) = fn_xz(1,1,kDw)
                0153       fp_yz(1,1,kDw) = fp_yz(1,1,kDw)
                0154       fn_yz(1,1,kDw) = fn_yz(1,1,kDw)
                0155 #endif
                0156
0733e7e26d Jean*0157       recip_dT = zeroRL
                0158       IF ( deltaTloc.GT.zeroRL ) recip_dT = 1.0 _d 0 / deltaTloc
                0159       noFlowAcrossSurf = rigidLid .OR. nonlinFreeSurf.GE.1
                0160      &                            .OR. select_rStar.NE.0
9822905e7f Jean*0161
1574069d50 Mart*0162 #ifdef ALLOW_AUTODIFF_TAMC
                0163 CADJ STORE sm_o,sm_v,sm_x,sm_xz,sm_y,sm_yz,sm_z,sm_zz
                0164 CADJ &     = somtape_r, key = 1, kind = isbyte
                0165 #endif
d7ce0d34f8 Jean*0166 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0167 C---  part.1 : calculate flux for all moments
8aa8d99107 Jean*0168       DO j=jMinAdvR,jMaxAdvR
                0169         DO i=iMinAdvR,iMaxAdvR
d7ce0d34f8 Jean*0170           wLoc = rTrans(i,j)*deltaTloc
                0171 C--    Flux from (k) to (k-1) when W>0 (i.e., take upper side of box k)
                0172 C- note: Linear free surface case: this takes care of w_surf advection out
                0173 C       of the domain since for this particular case, rTrans is not masked
0733e7e26d Jean*0174           fp_v (i,j,kUp) = MAX( zeroRL,  wLoc )
d7ce0d34f8 Jean*0175           alp  (i,j,kUp) = fp_v(i,j,kUp)/sm_v(i,j,k)
                0176           alpq           = alp(i,j,kUp)*alp(i,j,kUp)
                0177           alp1           = 1. _d 0 - alp(i,j,kUp)
                0178 C-     Create temporary moments/masses for partial boxes in transit
                0179 C       use same indexing as velocity, "p" for positive W
                0180           fp_o (i,j,kUp) = alp(i,j,kUp)*
                0181      &                   ( sm_o(i,j, k ) + alp1*sm_z(i,j, k )
                0182      &                   + alp1*(alp1-alp(i,j,kUp))*sm_zz(i,j, k )
                0183      &                   )
                0184           fp_z (i,j,kUp) = alpq*
                0185      &                   ( sm_z(i,j, k ) + three*alp1*sm_zz(i,j, k ) )
                0186           fp_zz(i,j,kUp) = alp(i,j,kUp)*alpq*sm_zz(i,j, k )
                0187           fp_x (i,j,kUp) = alp(i,j,kUp)*
                0188      &                   ( sm_x(i,j, k ) + alp1*sm_xz(i,j, k ) )
                0189           fp_y (i,j,kUp) = alp(i,j,kUp)*
                0190      &                   ( sm_y(i,j, k ) + alp1*sm_yz(i,j, k ) )
                0191           fp_xz(i,j,kUp) = alpq        *sm_xz(i,j, k )
                0192           fp_yz(i,j,kUp) = alpq        *sm_yz(i,j, k )
                0193           fp_xx(i,j,kUp) = alp(i,j,kUp)*sm_xx(i,j, k )
                0194           fp_yy(i,j,kUp) = alp(i,j,kUp)*sm_yy(i,j, k )
                0195           fp_xy(i,j,kUp) = alp(i,j,kUp)*sm_xy(i,j, k )
                0196         ENDDO
                0197       ENDDO
                0198       IF ( k.EQ.1 ) THEN
                0199 C--   Linear free surface, calculate w_surf (<0) advection term
                0200        km1 = 1
8aa8d99107 Jean*0201        DO j=jMinAdvR,jMaxAdvR
                0202         DO i=iMinAdvR,iMaxAdvR
d7ce0d34f8 Jean*0203           wLoc = rTrans(i,j)*deltaTloc
                0204 C-     Flux from above to (k) when W<0 , surface case:
                0205 C      take box k=1, assuming zero 1rst & 2nd moment in Z dir.
0733e7e26d Jean*0206           fn_v (i,j,kUp) = MAX( zeroRL, -wLoc )
d7ce0d34f8 Jean*0207           aln  (i,j,kUp) = fn_v(i,j,kUp)/sm_v(i,j,km1)
                0208           alnq           = aln(i,j,kUp)*aln(i,j,kUp)
                0209           aln1           = 1. _d 0 - aln(i,j,kUp)
                0210 C-     Create temporary moments/masses for partial boxes in transit
                0211 C       use same indexing as velocity, "n" for negative W
                0212           fn_o (i,j,kUp) = aln(i,j,kUp)*sm_o(i,j,km1)
0733e7e26d Jean*0213           fn_z (i,j,kUp) = zeroRL
                0214           fn_zz(i,j,kUp) = zeroRL
d7ce0d34f8 Jean*0215           fn_x (i,j,kUp) = aln(i,j,kUp)*sm_x(i,j,km1)
                0216           fn_y (i,j,kUp) = aln(i,j,kUp)*sm_y(i,j,km1)
0733e7e26d Jean*0217           fn_xz(i,j,kUp) = zeroRL
                0218           fn_yz(i,j,kUp) = zeroRL
d7ce0d34f8 Jean*0219           fn_xx(i,j,kUp) = aln(i,j,kUp)*sm_xx(i,j,km1)
                0220           fn_yy(i,j,kUp) = aln(i,j,kUp)*sm_yy(i,j,km1)
                0221           fn_xy(i,j,kUp) = aln(i,j,kUp)*sm_xy(i,j,km1)
                0222 C--    Save zero-order flux:
9822905e7f Jean*0223           wT(i,j) = ( fp_o(i,j,kUp) - fn_o(i,j,kUp) )*recip_dT
d7ce0d34f8 Jean*0224         ENDDO
                0225        ENDDO
                0226       ELSE
                0227 C--   Interior only: mask rTrans (if not already done)
                0228        km1 = k-1
8aa8d99107 Jean*0229        DO j=jMinAdvR,jMaxAdvR
                0230         DO i=iMinAdvR,iMaxAdvR
d7ce0d34f8 Jean*0231           wLoc = maskUp(i,j)*rTrans(i,j)*deltaTloc
                0232 C-     Flux from (k-1) to (k) when W<0 (i.e., take lower side of box k-1)
0733e7e26d Jean*0233           fn_v (i,j,kUp) = MAX( zeroRL, -wLoc )
d7ce0d34f8 Jean*0234           aln  (i,j,kUp) = fn_v(i,j,kUp)/sm_v(i,j,km1)
                0235           alnq           = aln(i,j,kUp)*aln(i,j,kUp)
                0236           aln1           = 1. _d 0 - aln(i,j,kUp)
                0237 C-     Create temporary moments/masses for partial boxes in transit
                0238 C       use same indexing as velocity, "n" for negative W
                0239           fn_o (i,j,kUp) = aln(i,j,kUp)*
                0240      &                   ( sm_o(i,j,km1) - aln1*sm_z(i,j,km1)
                0241      &                   + aln1*(aln1-aln(i,j,kUp))*sm_zz(i,j,km1)
                0242      &                   )
                0243           fn_z (i,j,kUp) = alnq*
                0244      &                   ( sm_z(i,j,km1) - three*aln1*sm_zz(i,j,km1) )
                0245           fn_zz(i,j,kUp) = aln(i,j,kUp)*alnq*sm_zz(i,j,km1)
                0246           fn_x (i,j,kUp) = aln(i,j,kUp)*
                0247      &                   ( sm_x(i,j,km1) - aln1*sm_xz(i,j,km1) )
                0248           fn_y (i,j,kUp) = aln(i,j,kUp)*
                0249      &                   ( sm_y(i,j,km1) - aln1*sm_yz(i,j,km1) )
                0250           fn_xz(i,j,kUp) = alnq        *sm_xz(i,j,km1)
                0251           fn_yz(i,j,kUp) = alnq        *sm_yz(i,j,km1)
                0252           fn_xx(i,j,kUp) = aln(i,j,kUp)*sm_xx(i,j,km1)
                0253           fn_yy(i,j,kUp) = aln(i,j,kUp)*sm_yy(i,j,km1)
                0254           fn_xy(i,j,kUp) = aln(i,j,kUp)*sm_xy(i,j,km1)
                0255 C--    Save zero-order flux:
9822905e7f Jean*0256           wT(i,j) = ( fp_o(i,j,kUp) - fn_o(i,j,kUp) )*recip_dT
d7ce0d34f8 Jean*0257         ENDDO
                0258        ENDDO
                0259 C--   end surface/interior cases for W<0 advective fluxes
                0260       ENDIF
0733e7e26d Jean*0261       IF ( .NOT.uniformFreeSurfLev .AND. k.NE.1 .AND.
                0262      &     .NOT.noFlowAcrossSurf ) THEN
d7ce0d34f8 Jean*0263 C--   Linear free surface, but surface not @ k=1 :
                0264 C     calculate w_surf (<0) advection term from current level
                0265 C     moments assuming zero 1rst & 2nd moment in Z dir. ;
                0266 C     and add to previous fluxes; note: identical to resetting fluxes
                0267 C     since previous fluxes are zero in this case (=> let TAF decide)
                0268        km1 = k
8aa8d99107 Jean*0269        DO j=jMinAdvR,jMaxAdvR
                0270         DO i=iMinAdvR,iMaxAdvR
d7ce0d34f8 Jean*0271          wLoc = rTrans(i,j)*deltaTloc
0733e7e26d Jean*0272          IF ( k.EQ.kSurfC(i,j,bi,bj) ) THEN
d7ce0d34f8 Jean*0273 C-     Flux from (k-1) to (k) when W<0 (special surface case, take box k)
0733e7e26d Jean*0274           fn_v (i,j,kUp) = MAX( zeroRL, -wLoc )
d7ce0d34f8 Jean*0275           aln  (i,j,kUp) = fn_v(i,j,kUp)/sm_v(i,j,km1)
                0276 C-     Create temporary moments/masses for partial boxes in transit
                0277 C       use same indexing as velocity, "n" for negative W
                0278           fn_o (i,j,kUp) = aln(i,j,kUp)*sm_o(i,j,km1)
                0279           fn_x (i,j,kUp) = aln(i,j,kUp)*sm_x(i,j,km1)
                0280           fn_y (i,j,kUp) = aln(i,j,kUp)*sm_y(i,j,km1)
                0281           fn_xx(i,j,kUp) = aln(i,j,kUp)*sm_xx(i,j,km1)
                0282           fn_yy(i,j,kUp) = aln(i,j,kUp)*sm_yy(i,j,km1)
                0283           fn_xy(i,j,kUp) = aln(i,j,kUp)*sm_xy(i,j,km1)
                0284 C--    Save zero-order flux:
9822905e7f Jean*0285           wT(i,j) = ( fp_o(i,j,kUp) - fn_o(i,j,kUp) )*recip_dT
d7ce0d34f8 Jean*0286          ENDIF
                0287         ENDDO
                0288        ENDDO
                0289       ENDIF
                0290
                0291 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0292 C---  part.2 : re-adjust moments remaining in the box
                0293 C      take off from grid box (k): negative W(kDw) and positive W(kUp)
8aa8d99107 Jean*0294       DO j=jMinAdvR,jMaxAdvR
                0295        DO i=iMinAdvR,iMaxAdvR
72de869c1b Jean*0296 #ifdef ALLOW_OBCS
0733e7e26d Jean*0297         IF ( maskIn(i,j).NE.zeroRS ) THEN
72de869c1b Jean*0298 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0299         alf1  = 1. _d 0 - aln(i,j,kDw) - alp(i,j,kUp)
                0300         alf1q = alf1*alf1
                0301         alpmn = alp(i,j,kUp) - aln(i,j,kDw)
                0302         sm_v (i,j,k) = sm_v (i,j,k) - fn_v (i,j,kDw) - fp_v (i,j,kUp)
                0303         sm_o (i,j,k) = sm_o (i,j,k) - fn_o (i,j,kDw) - fp_o (i,j,kUp)
                0304         sm_z (i,j,k) = alf1q*( sm_z(i,j,k) - three*alpmn*sm_zz(i,j,k) )
                0305         sm_zz(i,j,k) = alf1*alf1q*sm_zz(i,j,k)
                0306         sm_xz(i,j,k) = alf1q*sm_xz(i,j,k)
                0307         sm_yz(i,j,k) = alf1q*sm_yz(i,j,k)
                0308         sm_x (i,j,k) = sm_x (i,j,k) - fn_x (i,j,kDw) - fp_x (i,j,kUp)
                0309         sm_xx(i,j,k) = sm_xx(i,j,k) - fn_xx(i,j,kDw) - fp_xx(i,j,kUp)
                0310         sm_y (i,j,k) = sm_y (i,j,k) - fn_y (i,j,kDw) - fp_y (i,j,kUp)
                0311         sm_yy(i,j,k) = sm_yy(i,j,k) - fn_yy(i,j,kDw) - fp_yy(i,j,kUp)
                0312         sm_xy(i,j,k) = sm_xy(i,j,k) - fn_xy(i,j,kDw) - fp_xy(i,j,kUp)
72de869c1b Jean*0313 #ifdef ALLOW_OBCS
                0314         ENDIF
                0315 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0316        ENDDO
                0317       ENDDO
                0318
                0319 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0320 C---  part.3 : Put the temporary moments into appropriate neighboring boxes
                0321 C      add into grid box (k): positive W(kDw) and negative W(kUp)
8aa8d99107 Jean*0322       DO j=jMinAdvR,jMaxAdvR
                0323        DO i=iMinAdvR,iMaxAdvR
72de869c1b Jean*0324 #ifdef ALLOW_OBCS
0733e7e26d Jean*0325         IF ( maskIn(i,j).NE.zeroRS ) THEN
72de869c1b Jean*0326 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0327         sm_v (i,j,k) = sm_v (i,j,k) + fp_v (i,j,kDw) + fn_v (i,j,kUp)
                0328         alfp = fp_v(i,j,kDw)/sm_v(i,j,k)
                0329         alfn = fn_v(i,j,kUp)/sm_v(i,j,k)
                0330         alf1 = 1. _d 0 - alfp - alfn
                0331         alp1 = 1. _d 0 - alfp
                0332         aln1 = 1. _d 0 - alfn
                0333         alpmn = alfp - alfn
                0334         locTp = alfp*sm_o(i,j,k) - alp1*fp_o(i,j,kDw)
                0335         locTn = alfn*sm_o(i,j,k) - aln1*fn_o(i,j,kUp)
                0336         sm_zz(i,j,k) = alf1*alf1*sm_zz(i,j,k) + alfp*alfp*fp_zz(i,j,kDw)
                0337      &                                        + alfn*alfn*fn_zz(i,j,kUp)
                0338      &     - 5. _d 0*(-alpmn*alf1*sm_z(i,j,k) + alfp*alp1*fp_z(i,j,kDw)
                0339      &                                        - alfn*aln1*fn_z(i,j,kUp)
0733e7e26d Jean*0340      &               + twoRL*alfp*alfn*sm_o(i,j,k) + (alp1-alfp)*locTp
                0341      &                                             + (aln1-alfn)*locTn
d7ce0d34f8 Jean*0342      &               )
                0343         sm_xz(i,j,k) = alf1*sm_xz(i,j,k) + alfp*fp_xz(i,j,kDw)
                0344      &                                   + alfn*fn_xz(i,j,kUp)
                0345      &       + three*( alpmn*sm_x(i,j,k) - alp1*fp_x(i,j,kDw)
                0346      &                                   + aln1*fn_x(i,j,kUp)
                0347      &               )
                0348         sm_yz(i,j,k) = alf1*sm_yz(i,j,k) + alfp*fp_yz(i,j,kDw)
                0349      &                                   + alfn*fn_yz(i,j,kUp)
                0350      &       + three*( alpmn*sm_y(i,j,k) - alp1*fp_y(i,j,kDw)
                0351      &                                   + aln1*fn_y(i,j,kUp)
                0352      &               )
                0353         sm_z (i,j,k) = alf1*sm_z(i,j,k) + alfp*fp_z(i,j,kDw)
                0354      &                                  + alfn*fn_z(i,j,kUp)
                0355      &               + three*( locTp - locTn )
                0356         sm_o (i,j,k) = sm_o (i,j,k) + fp_o (i,j,kDw) + fn_o (i,j,kUp)
                0357         sm_x (i,j,k) = sm_x (i,j,k) + fp_x (i,j,kDw) + fn_x (i,j,kUp)
                0358         sm_xx(i,j,k) = sm_xx(i,j,k) + fp_xx(i,j,kDw) + fn_xx(i,j,kUp)
                0359         sm_y (i,j,k) = sm_y (i,j,k) + fp_y (i,j,kDw) + fn_y (i,j,kUp)
                0360         sm_yy(i,j,k) = sm_yy(i,j,k) + fp_yy(i,j,kDw) + fn_yy(i,j,kUp)
                0361         sm_xy(i,j,k) = sm_xy(i,j,k) + fp_xy(i,j,kDw) + fn_xy(i,j,kUp)
72de869c1b Jean*0362 #ifdef ALLOW_OBCS
                0363         ENDIF
                0364 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0365        ENDDO
                0366       ENDDO
7448700841 Mart*0367 #endif /* GAD_ALLOW_TS_SOM_ADV or PTRACERS_ALLOW_DYN_STATE */
d7ce0d34f8 Jean*0368
                0369       RETURN
                0370       END