pkg/generic_advdiff/gad_ppm_flx_x.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PPM_FLX_X(bi,bj,kk,iy,
                0004      &           calc_CFL,delT,uvel,
                0005      &           ufac,fhat,flux,myThid)
                0006 C     |================================================================|
                0007 C     | PPM_FLX_X: evaluate PPM flux on grid-cell edges.               |
                0008 C     |================================================================|
                0009
                0010           implicit none
                0011
                0012 C     =============================================== global variables
                0013 #         include "SIZE.h"
                0014 #         include "GRID.h"
                0015 #         include "GAD.h"
                0016
                0017 C     ================================================================
                0018 C       bi,bj     :: tile indexing.
                0019 C       kk        :: r-index.
                0020 C       iy        :: y-index.
                0021 C       calc_CFL  :: TRUE to calc. CFL from vel.
                0022 C       delT      :: time-step.
                0023 C       uvel      :: vel.-comp in x-direction.
                0024 C       ufac      :: vel.-flux in x-direction.
                0025 C       fhat      :: row of poly. coeff.
                0026 C       flux      :: adv.-flux in x-direction.
                0027 C       myThid    :: thread number.
                0028 C     ================================================================
                0029           integer bi,bj,kk,iy
                0030           logical calc_CFL
                0031           _RL delT
                0032           _RL uvel(1-OLx:sNx+OLx,
                0033      &             1-OLy:sNy+OLy)
                0034           _RL ufac(1-OLx:sNx+OLx,
                0035      &             1-OLy:sNy+OLy)
                0036           _RL fhat(1:3,
                0037      &             1-OLx:sNx+OLx)
                0038           _RL flux(1-OLx:sNx+OLx,
                0039      &             1-OLy:sNy+OLy)
                0040           integer myThid
                0041
                0042 C     ================================================================
                0043 C       ix        :: x-indexing.
                0044 C       uCFL      :: CFL number.
                0045 C       intF      :: upwind tracer edge-value.
                0046 C       ss11,ss22 :: int. endpoints.
                0047 C       ivec      :: int. basis vec.
                0048 C     ================================================================
                0049           integer ix
                0050           _RL uCFL,intF
                0051           _RL ss11,ss22
                0052           _RL ivec(1:3)
                0053
                0054 C     ================================================================
                0055 C       (1): calc. "departure-points" for each grid-cell edge by int-
                0056 C            egrating edge position backward in time over one single
                0057 C            time-step. This is a "single-cell" implementation: requ-
                0058 C            ires CFL <= 1.0.
                0059 C       (2): calc. flux as the integral of the upwind grid-cell poly-
                0060 C            nomial over the deformation interval found in (1).
                0061 C     ================================================================
                0062
                0063           do  ix = 1-OLx+3, sNx+OLx-2
                0064
                0065               if (uvel(ix,iy) .eq. 0. _d 0) then
                0066
                0067                   flux(ix,iy)    = 0. _d 0
                0068
                0069               else
                0070
                0071               if (uvel(ix,iy) .gt. 0. _d 0) then
                0072
                0073 C     ==================== integrate PPM profile over upwind cell IX-1
                0074               if ( calc_CFL ) then
                0075               uCFL = uvel(ix,iy) * delT
                0076      &             * recip_dxF(ix-1,iy,bi,bj)
                0077      &             * recip_deepFacC(kk)
                0078               else
                0079               uCFL = uvel(ix,iy)
                0080               end if
                0081
                0082               ss11 = +1. _d 0 - 2. _d 0 * uCFL
                0083               ss22 = +1. _d 0
                0084
                0085 C     ==================== integrate profile over region swept by face
                0086               ivec(1) = ss22 - ss11
                0087               ivec(2) =(ss22 ** 2
                0088      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0089               ivec(3) =(ss22 ** 3
                0090      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0091
                0092               intF = ivec(1) * fhat(1,ix-1)
                0093      &             + ivec(2) * fhat(2,ix-1)
                0094      &             + ivec(3) * fhat(3,ix-1)
                0095
83ddf5a6c6 Mart*0096 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0097
                0098               else
                0099
                0100 C     ==================== integrate PPM profile over upwind cell IX+0
                0101               if ( calc_CFL ) then
                0102               uCFL = uvel(ix,iy) * delT
                0103      &             * recip_dxF(ix-0,iy,bi,bj)
                0104      &             * recip_deepFacC(kk)
                0105               else
                0106               uCFL = uvel(ix,iy)
                0107               end if
                0108
                0109               ss11 = -1. _d 0 - 2. _d 0 * uCFL
                0110               ss22 = -1. _d 0
                0111
                0112 C     ==================== integrate profile over region swept by face
                0113               ivec(1) = ss22 - ss11
                0114               ivec(2) =(ss22 ** 2
                0115      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0116               ivec(3) =(ss22 ** 3
                0117      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0118
                0119               intF = ivec(1) * fhat(1,ix-0)
                0120      &             + ivec(2) * fhat(2,ix-0)
                0121      &             + ivec(3) * fhat(3,ix-0)
                0122
83ddf5a6c6 Mart*0123 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0124
                0125               end if
                0126
83ddf5a6c6 Mart*0127 c             intF = 0.5 _d 0 * intF / uCFL
                0128 C-    to avoid potential underflow:
                0129               intF = 0.5 _d 0 * intF / sign(max(abs(uCFL),1.d-20),uCFL)
                0130
8e4c181d69 Jean*0131 C     ==================== calc. flux = upwind tracer * face-transport
                0132               flux(ix,iy) = + ufac(ix,iy) * intF
                0133
                0134               end if
                0135
                0136           end do
                0137
                0138           return
                0139
                0140 c     end subroutine GAD_PPM_FLX_X
                0141       end