Back to home page

MITgcm

 
 

    


File indexing completed on 2023-09-03 05:10:24 UTC

view on githubraw file Latest commit 74487008 on 2023-09-03 01:50:18 UTC
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_Y
                0011 
                0012 C !INTERFACE: ==========================================================
                0013       SUBROUTINE GAD_SOM_ADV_Y(
                0014      I           bi,bj,k, limiter,
b79a2b44f2 Jean*0015      I           overlapOnly, interiorOnly,
                0016      I           N_edge, S_edge, E_edge, W_edge,
72de869c1b Jean*0017      I           deltaTloc, vTrans, maskIn,
d7ce0d34f8 Jean*0018      U           sm_v, sm_o, sm_x, sm_y, sm_z,
                0019      U           sm_xx, sm_yy, sm_zz, sm_xy, sm_xz, sm_yz,
                0020      O           vT,
                0021      I           myThid )
                0022 
                0023 C !DESCRIPTION:
                0024 C  Calculates the area integrated meridional flux due to advection
                0025 C  of a tracer using
                0026 C--
                0027 C        Second-Order Moments Advection of tracer in Y-direction
                0028 C        ref: M.J.Prather, 1986, JGR, 91, D6, pp 6671-6681.
                0029 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0030 C      The 3-D grid has dimension  (Nx,Ny,Nz) with corresponding
                0031 C      velocity field (U,V,W).  Parallel subroutine calculate
                0032 C      advection in the X- and Z- directions.
                0033 C      The moment [Si] are as defined in the text, Sm refers to
                0034 C      the total mass in each grid box
                0035 C      the moments [Fi] are similarly defined and used as temporary
                0036 C      storage for portions of the grid boxes in transit.
                0037 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0038 
                0039 C !USES: ===============================================================
                0040       IMPLICIT NONE
                0041 #include "SIZE.h"
00fdbdcbd5 Jean*0042 #include "EEPARAMS.h"
d7ce0d34f8 Jean*0043 #include "GAD.h"
1574069d50 Mart*0044 #ifdef ALLOW_AUTODIFF_TAMC
                0045 # include "tamc.h"
                0046 #endif
d7ce0d34f8 Jean*0047 
                0048 C !INPUT PARAMETERS: ===================================================
b79a2b44f2 Jean*0049 C  bi,bj         :: tile indices
                0050 C  k             :: vertical level
                0051 C  limiter       :: 0: no limiter ; 1: Prather, 1986 limiter
                0052 C  overlapOnly   :: only update the edges of myTile, but not the interior
                0053 C  interiorOnly  :: only update the interior of myTile, but not the edges
                0054 C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
                0055 C  vTrans        :: zonal volume transport
72de869c1b Jean*0056 C  maskIn        :: 2-D array Interior mask
b79a2b44f2 Jean*0057 C  myThid        :: my Thread Id. number
d7ce0d34f8 Jean*0058       INTEGER bi,bj,k
                0059       INTEGER limiter
b79a2b44f2 Jean*0060       LOGICAL overlapOnly, interiorOnly
                0061       LOGICAL N_edge, S_edge, E_edge, W_edge
d7ce0d34f8 Jean*0062       _RL deltaTloc
                0063       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
72de869c1b Jean*0064       _RS maskIn(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
d7ce0d34f8 Jean*0065       INTEGER myThid
                0066 
                0067 C !OUTPUT PARAMETERS: ==================================================
                0068 C  sm_v         :: volume of grid cell
                0069 C  sm_o         :: tracer content of grid cell (zero order moment)
                0070 C  sm_x,y,z     :: 1rst order moment of tracer distribution, in x,y,z direction
                0071 C  sm_xx,yy,zz  ::  2nd order moment of tracer distribution, in x,y,z direction
                0072 C  sm_xy,xz,yz  ::  2nd order moment of tracer distr., in cross direction xy,xz,yz
                0073 C  vT           :: meridional advective flux
                0074       _RL sm_v  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0075       _RL sm_o  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0076       _RL sm_x  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0077       _RL sm_y  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0078       _RL sm_z  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL sm_xx (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL sm_yy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL sm_zz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL sm_xy (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL sm_xz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL sm_yz (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL vT    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086 
7448700841 Mart*0087 #if ( defined GAD_ALLOW_TS_SOM_ADV || defined PTRACERS_ALLOW_DYN_STATE )
d7ce0d34f8 Jean*0088 C !LOCAL VARIABLES: ====================================================
b79a2b44f2 Jean*0089 C  i,j           :: loop indices
                0090 C  vLoc          :: volume transported (per time step)
                0091 C [iMin,iMax]Upd :: loop range to update tracer field
                0092 C [jMin,jMax]Upd :: loop range to update tracer field
                0093 C  nbStrips      :: number of strips (if region to update is splitted)
00fdbdcbd5 Jean*0094       _RL three
d7ce0d34f8 Jean*0095       PARAMETER( three = 3. _d 0 )
                0096       INTEGER i,j
b79a2b44f2 Jean*0097       INTEGER ns, nbStrips
                0098       INTEGER iMinUpd(2), iMaxUpd(2), jMinUpd(2), jMaxUpd(2)
9822905e7f Jean*0099       _RL  recip_dT
d7ce0d34f8 Jean*0100       _RL  slpmax, s1max, s1new, s2new
                0101       _RL  vLoc, alf1, alf1q, alpmn
                0102       _RL  alfp, alpq, alp1, locTp
                0103       _RL  alfn, alnq, aln1, locTn
                0104       _RL  alp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0105       _RL  aln  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0106       _RL  fp_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0107       _RL  fn_v (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0108       _RL  fp_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0109       _RL  fn_o (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0110       _RL  fp_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0111       _RL  fn_x (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0112       _RL  fp_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0113       _RL  fn_y (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0114       _RL  fp_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0115       _RL  fn_z (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0116       _RL  fp_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0117       _RL  fn_xx(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0118       _RL  fp_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0119       _RL  fn_yy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0120       _RL  fp_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0121       _RL  fn_zz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0122       _RL  fp_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0123       _RL  fn_xy(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0124       _RL  fp_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0125       _RL  fn_xz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0126       _RL  fp_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0127       _RL  fn_yz(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0128 CEOP
                0129 
1574069d50 Mart*0130 #ifdef ALLOW_AUTODIFF_TAMC
                0131 CADJ INIT somtape_y = COMMON, 2
                0132 #endif
9822905e7f Jean*0133       recip_dT = 0.
00fdbdcbd5 Jean*0134       IF ( deltaTloc.GT.zeroRL ) recip_dT = 1.0 _d 0 / deltaTloc
9822905e7f Jean*0135 
b79a2b44f2 Jean*0136 C-    Set loop ranges for updating tracer field (splitted in 2 strips)
                0137       nbStrips   = 1
72de869c1b Jean*0138       iMinUpd(1) = 1-OLx
                0139       iMaxUpd(1) = sNx+OLx
                0140       jMinUpd(1) = 1-OLy+1
                0141       jMaxUpd(1) = sNy+OLy-1
b79a2b44f2 Jean*0142       IF ( overlapOnly ) THEN
                0143 C     update in overlap-Only
                0144         IF ( S_edge ) jMinUpd(1) = 1
                0145         IF ( N_edge ) jMaxUpd(1) = sNy
                0146         IF ( W_edge ) THEN
72de869c1b Jean*0147           iMinUpd(1) = 1-OLx
b79a2b44f2 Jean*0148           iMaxUpd(1) = 0
                0149         ENDIF
                0150         IF ( E_edge ) THEN
                0151           IF ( W_edge ) nbStrips = 2
                0152           iMinUpd(nbStrips) = sNx+1
72de869c1b Jean*0153           iMaxUpd(nbStrips) = sNx+OLx
b79a2b44f2 Jean*0154         ENDIF
                0155       ELSE
                0156 C     do not only update the overlap
                0157         IF ( interiorOnly .AND. W_edge ) iMinUpd(1) = 1
                0158         IF ( interiorOnly .AND. E_edge ) iMaxUpd(1) = sNx
                0159       ENDIF
                0160 
                0161 C--   start 1rst loop on strip number "ns"
                0162       DO ns=1,nbStrips
                0163 
d7ce0d34f8 Jean*0164       IF ( limiter.EQ.1 ) THEN
1574069d50 Mart*0165 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0166 CADJ STORE sm_o,sm_y,sm_yy,sm_xy,sm_yz
1574069d50 Mart*0167 CADJ &     = somtape_y, key = ns, kind = isbyte
                0168 #endif
b79a2b44f2 Jean*0169        DO j=jMinUpd(1)-1,jMaxUpd(1)+1
                0170         DO i=iMinUpd(ns),iMaxUpd(ns)
d7ce0d34f8 Jean*0171 C     If flux-limiting transport is to be applied, place limits on
                0172 C     appropriate moments before transport.
                0173          slpmax = 0.
                0174          IF ( sm_o(i,j).GT.0. ) slpmax = sm_o(i,j)
                0175          s1max = slpmax*1.5 _d 0
                0176          s1new = MIN(  s1max, MAX(-s1max,sm_y(i,j)) )
                0177          s2new = MIN( (slpmax+slpmax-ABS(s1new)/three),
                0178      &                MAX(ABS(s1new)-slpmax,sm_yy(i,j))  )
                0179          sm_xy(i,j) = MIN( slpmax, MAX(-slpmax,sm_xy(i,j)) )
                0180          sm_yz(i,j) = MIN( slpmax, MAX(-slpmax,sm_yz(i,j)) )
a5a86b736b Jean*0181          sm_y (i,j) = s1new
                0182          sm_yy(i,j) = s2new
d7ce0d34f8 Jean*0183         ENDDO
                0184        ENDDO
                0185       ENDIF
                0186 
1574069d50 Mart*0187 #ifdef ALLOW_AUTODIFF_TAMC
                0188 CADJ STORE sm_o,sm_v,sm_x,sm_xx,sm_xy,sm_xz,sm_y,sm_yy,sm_yz,sm_z,sm_zz
                0189 CADJ &     = somtape_y, key = ns, kind = isbyte
                0190 #endif
d7ce0d34f8 Jean*0191 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0192 C---  part.1 : calculate flux for all moments
b79a2b44f2 Jean*0193       DO j=jMinUpd(1),jMaxUpd(1)+1
                0194        DO i=iMinUpd(ns),iMaxUpd(ns)
d7ce0d34f8 Jean*0195         vLoc = vTrans(i,j)*deltaTloc
                0196 C--    Flux from (j-1) to (j) when V>0 (i.e., take right side of box j-1)
00fdbdcbd5 Jean*0197         fp_v (i,j) = MAX( zeroRL,  vLoc )
d7ce0d34f8 Jean*0198         alp  (i,j) = fp_v(i,j)/sm_v(i,j-1)
                0199         alpq       = alp(i,j)*alp(i,j)
                0200         alp1       = 1. _d 0 - alp(i,j)
                0201 C-     Create temporary moments/masses for partial boxes in transit
                0202 C       use same indexing as velocity, "p" for positive V
                0203         fp_o (i,j) = alp(i,j)*( sm_o(i,j-1) + alp1*sm_y(i,j-1)
                0204      &                        + alp1*(alp1-alp(i,j))*sm_yy(i,j-1)
                0205      &                        )
                0206         fp_y (i,j) = alpq    *( sm_y(i,j-1) + three*alp1*sm_yy(i,j-1) )
                0207         fp_yy(i,j) = alp(i,j)*alpq*sm_yy(i,j-1)
                0208         fp_x (i,j) = alp(i,j)*( sm_x(i,j-1) + alp1*sm_xy(i,j-1) )
                0209         fp_z (i,j) = alp(i,j)*( sm_z(i,j-1) + alp1*sm_yz(i,j-1) )
                0210 
                0211         fp_xy(i,j) = alpq    *sm_xy(i,j-1)
                0212         fp_yz(i,j) = alpq    *sm_yz(i,j-1)
                0213         fp_xx(i,j) = alp(i,j)*sm_xx(i,j-1)
                0214         fp_zz(i,j) = alp(i,j)*sm_zz(i,j-1)
                0215         fp_xz(i,j) = alp(i,j)*sm_xz(i,j-1)
                0216 C--    Flux from (j) to (j-1) when V<0 (i.e., take left side of box j)
00fdbdcbd5 Jean*0217         fn_v (i,j) = MAX( zeroRL, -vLoc )
d7ce0d34f8 Jean*0218         aln  (i,j) = fn_v(i,j)/sm_v(i, j )
                0219         alnq       = aln(i,j)*aln(i,j)
                0220         aln1       = 1. _d 0 - aln(i,j)
                0221 C-     Create temporary moments/masses for partial boxes in transit
                0222 C       use same indexing as velocity, "n" for negative V
                0223         fn_o (i,j) = aln(i,j)*( sm_o(i, j ) - aln1*sm_y(i, j )
                0224      &                        + aln1*(aln1-aln(i,j))*sm_yy(i, j )
                0225      &                        )
                0226         fn_y (i,j) = alnq    *( sm_y(i, j ) - three*aln1*sm_yy(i, j ) )
                0227         fn_yy(i,j) = aln(i,j)*alnq*sm_yy(i, j )
                0228         fn_x (i,j) = aln(i,j)*( sm_x(i, j ) - aln1*sm_xy(i, j ) )
                0229         fn_z (i,j) = aln(i,j)*( sm_z(i, j ) - aln1*sm_yz(i, j ) )
                0230         fn_xy(i,j) = alnq    *sm_xy(i, j )
                0231         fn_yz(i,j) = alnq    *sm_yz(i, j )
                0232         fn_xx(i,j) = aln(i,j)*sm_xx(i, j )
                0233         fn_zz(i,j) = aln(i,j)*sm_zz(i, j )
                0234         fn_xz(i,j) = aln(i,j)*sm_xz(i, j )
                0235 C--    Save zero-order flux:
9822905e7f Jean*0236         vT(i,j) = ( fp_o(i,j) - fn_o(i,j) )*recip_dT
d7ce0d34f8 Jean*0237        ENDDO
                0238       ENDDO
                0239 
b79a2b44f2 Jean*0240 C--   end 1rst loop on strip number "ns"
                0241 c     ENDDO
                0242 
1574069d50 Mart*0243 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0244 CADJ STORE sm_o,sm_v,sm_x,sm_y,sm_z,sm_xy,sm_yy,sm_yz
1574069d50 Mart*0245 CADJ &     = somtape_y, key = ns, kind = isbyte
                0246 #endif
d7ce0d34f8 Jean*0247 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
b79a2b44f2 Jean*0248 C--   start 2nd loop on strip number "ns"
                0249 c     DO ns=1,nbStrips
                0250 
d7ce0d34f8 Jean*0251 C---  part.2 : re-adjust moments remaining in the box
                0252 C      take off from grid box (j): negative V(j) and positive V(j+1)
b79a2b44f2 Jean*0253       DO j=jMinUpd(1),jMaxUpd(1)
                0254        DO i=iMinUpd(ns),iMaxUpd(ns)
72de869c1b Jean*0255 #ifdef ALLOW_OBCS
00fdbdcbd5 Jean*0256         IF ( maskIn(i,j).NE.zeroRS ) THEN
72de869c1b Jean*0257 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0258         alf1  = 1. _d 0 - aln(i,j) - alp(i,j+1)
                0259         alf1q = alf1*alf1
                0260         alpmn = alp(i,j+1) - aln(i,j)
                0261         sm_v (i,j) = sm_v (i,j) - fn_v (i,j) - fp_v (i,j+1)
                0262         sm_o (i,j) = sm_o (i,j) - fn_o (i,j) - fp_o (i,j+1)
                0263         sm_y (i,j) = alf1q*( sm_y(i,j) - three*alpmn*sm_yy(i,j) )
                0264         sm_yy(i,j) = alf1*alf1q*sm_yy(i,j)
                0265         sm_xy(i,j) = alf1q*sm_xy(i,j)
                0266         sm_yz(i,j) = alf1q*sm_yz(i,j)
                0267         sm_x (i,j) = sm_x (i,j) - fn_x (i,j) - fp_x (i,j+1)
                0268         sm_xx(i,j) = sm_xx(i,j) - fn_xx(i,j) - fp_xx(i,j+1)
                0269         sm_z (i,j) = sm_z (i,j) - fn_z (i,j) - fp_z (i,j+1)
                0270         sm_zz(i,j) = sm_zz(i,j) - fn_zz(i,j) - fp_zz(i,j+1)
                0271         sm_xz(i,j) = sm_xz(i,j) - fn_xz(i,j) - fp_xz(i,j+1)
72de869c1b Jean*0272 #ifdef ALLOW_OBCS
                0273         ENDIF
                0274 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0275        ENDDO
                0276       ENDDO
                0277 
                0278 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0279 C---  part.3 : Put the temporary moments into appropriate neighboring boxes
                0280 C      add into grid box (j): positive V(j) and negative V(j+1)
b79a2b44f2 Jean*0281       DO j=jMinUpd(1),jMaxUpd(1)
                0282        DO i=iMinUpd(ns),iMaxUpd(ns)
72de869c1b Jean*0283 #ifdef ALLOW_OBCS
00fdbdcbd5 Jean*0284         IF ( maskIn(i,j).NE.zeroRS ) THEN
72de869c1b Jean*0285 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0286         sm_v (i,j) = sm_v (i,j) + fp_v (i,j) + fn_v (i,j+1)
                0287         alfp = fp_v(i, j )/sm_v(i,j)
                0288         alfn = fn_v(i,j+1)/sm_v(i,j)
                0289         alf1 = 1. _d 0 - alfp - alfn
                0290         alp1 = 1. _d 0 - alfp
                0291         aln1 = 1. _d 0 - alfn
                0292         alpmn = alfp - alfn
                0293         locTp = alfp*sm_o(i,j) - alp1*fp_o(i,j)
                0294         locTn = alfn*sm_o(i,j) - aln1*fn_o(i,j+1)
                0295         sm_yy(i,j) = alf1*alf1*sm_yy(i,j) + alfp*alfp*fp_yy(i,j)
                0296      &                                    + alfn*alfn*fn_yy(i,j+1)
                0297      &   - 5. _d 0*(-alpmn*alf1*sm_y(i,j) + alfp*alp1*fp_y(i,j)
                0298      &                                    - alfn*aln1*fn_y(i,j+1)
00fdbdcbd5 Jean*0299      &             + twoRL*alfp*alfn*sm_o(i,j) + (alp1-alfp)*locTp
                0300      &                                         + (aln1-alfn)*locTn
d7ce0d34f8 Jean*0301      &             )
                0302         sm_xy(i,j) = alf1*sm_xy(i,j) + alfp*fp_xy(i,j)
                0303      &                               + alfn*fn_xy(i,j+1)
                0304      &     + three*( alpmn*sm_x(i,j) - alp1*fp_x(i,j)
                0305      &                               + aln1*fn_x(i,j+1)
                0306      &             )
                0307         sm_yz(i,j) = alf1*sm_yz(i,j) + alfp*fp_yz(i,j)
                0308      &                               + alfn*fn_yz(i,j+1)
                0309      &     + three*( alpmn*sm_z(i,j) - alp1*fp_z(i,j)
                0310      &                               + aln1*fn_z(i,j+1)
                0311      &             )
                0312         sm_y (i,j) = alf1*sm_y(i,j) + alfp*fp_y(i,j) + alfn*fn_y(i,j+1)
                0313      &             + three*( locTp - locTn )
                0314         sm_o (i,j) = sm_o (i,j) + fp_o (i,j) + fn_o (i,j+1)
                0315         sm_x (i,j) = sm_x (i,j) + fp_x (i,j) + fn_x (i,j+1)
                0316         sm_xx(i,j) = sm_xx(i,j) + fp_xx(i,j) + fn_xx(i,j+1)
                0317         sm_z (i,j) = sm_z (i,j) + fp_z (i,j) + fn_z (i,j+1)
                0318         sm_zz(i,j) = sm_zz(i,j) + fp_zz(i,j) + fn_zz(i,j+1)
                0319         sm_xz(i,j) = sm_xz(i,j) + fp_xz(i,j) + fn_xz(i,j+1)
72de869c1b Jean*0320 #ifdef ALLOW_OBCS
                0321         ENDIF
                0322 #endif /* ALLOW_OBCS */
d7ce0d34f8 Jean*0323        ENDDO
                0324       ENDDO
                0325 
b79a2b44f2 Jean*0326 C--   end 2nd loop on strip number "ns"
                0327       ENDDO
7448700841 Mart*0328 #endif /* GAD_ALLOW_TS_SOM_ADV or PTRACERS_ALLOW_DYN_STATE */
b79a2b44f2 Jean*0329 
d7ce0d34f8 Jean*0330       RETURN
                0331       END