pkg/generic_advdiff/gad_ppm_flx_r.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PPM_FLX_R(bi,bj,ix,iy,
                0004      &           delT,wvel,wfac,fhat,
                0005      &           flux,myThid )
                0006 C     |================================================================|
                0007 C     | PPM_FLX_R: 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       ix        :: x-index.
                0020 C       iy        :: y-index.
                0021 C       delT      :: time-step.
                0022 C       wvel      :: vel.-comp in r-direction.
                0023 C       wfac      :: vel.-flux in r-direction.
                0024 C       fhat      :: row of poly. coeff.
                0025 C       flux      :: adv.-flux in r-direction.
                0026 C       myThid    :: thread number.
                0027 C     ================================================================
                0028           integer bi,bj,ix,iy
                0029           _RL delT(1:Nr)
                0030           _RL wvel(1-OLx:sNx+OLx,
                0031      &             1-OLy:sNy+OLy,
                0032      &             1:Nr)
                0033           _RL wfac(1-OLx:sNx+OLx,
                0034      &             1-OLy:sNy+OLy,
                0035      &             1:Nr)
                0036           _RL fhat(1:3 ,
                0037      &             1:Nr)
                0038           _RL flux(1-OLx:sNx+OLx,
                0039      &             1-OLy:sNy+OLy,
                0040      &             1:Nr)
                0041           integer myThid
                0042
                0043 C     ================================================================
                0044 C       ir        :: r-indexing.
                0045 C       wCFL      :: CFL number.
                0046 C       intF      :: upwind tracer edge-value.
                0047 C       ss11,ss22 :: int. endpoints.
                0048 C       ivec      :: int. basis vec.
                0049 C     ================================================================
                0050           integer ir
                0051           _RL wCFL,intF
                0052           _RL ss11,ss22
                0053           _RL ivec(1:3)
                0054
                0055 C     ================================================================
                0056 C       (1): calc. "departure-points" for each grid-cell edge by int-
                0057 C            egrating edge position backward in time over one single
                0058 C            time-step. This is a "single-cell" implementation: requ-
                0059 C            ires CFL <= 1.0.
                0060 C       (2): calc. flux as the integral of the upwind grid-cell poly-
                0061 C            nomial over the deformation interval found in (1).
                0062 C     ================================================================
                0063
                0064           do  ir = +2, Nr
                0065
                0066               if (wvel(ix,iy,ir) .eq. 0. _d 0) then
                0067
                0068                   flux(ix,iy,ir)    = 0. _d 0
                0069
                0070               else
                0071
                0072               if (wvel(ix,iy,ir) .lt. 0. _d 0) then
                0073
                0074 C     ==================== integrate PQM profile over upwind cell IR-1
                0075               wCFL = wvel(ix,iy,ir)
                0076      &             * delT(ir-1)*recip_drF(ir-1)
                0077
                0078               ss11 = +1. _d 0 + 2. _d 0 * wCFL
                0079               ss22 = +1. _d 0
                0080
                0081 C     ==================== integrate profile over region swept by face
                0082               ivec(1) = ss22 - ss11
                0083               ivec(2) =(ss22 ** 2
                0084      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0085               ivec(3) =(ss22 ** 3
                0086      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0087
                0088               intF = ivec(1) * fhat(1,ir-1)
                0089      &             + ivec(2) * fhat(2,ir-1)
                0090      &             + ivec(3) * fhat(3,ir-1)
                0091
83ddf5a6c6 Mart*0092 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0093
                0094               else
                0095
                0096 C     ==================== integrate PQM profile over upwind cell IR+0
                0097               wCFL = wvel(ix,iy,ir)
                0098      &             * delT(ir-0)*recip_drF(ir-0)
                0099
                0100               ss11 = -1. _d 0 + 2. _d 0 * wCFL
                0101               ss22 = -1. _d 0
                0102
                0103 C     ==================== integrate profile over region swept by face
                0104               ivec(1) = ss22 - ss11
                0105               ivec(2) =(ss22 ** 2
                0106      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0107               ivec(3) =(ss22 ** 3
                0108      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0109
                0110               intF = ivec(1) * fhat(1,ir-0)
                0111      &             + ivec(2) * fhat(2,ir-0)
                0112      &             + ivec(3) * fhat(3,ir-0)
                0113
83ddf5a6c6 Mart*0114 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0115
                0116               end if
                0117
83ddf5a6c6 Mart*0118 c             intF = -0.5 _d 0 * intF / wCFL
                0119 C-    to avoid potential underflow:
                0120               intF = -0.5 _d 0 * intF/sign(max(abs(wCFL),1.d-20),wCFL)
                0121
8e4c181d69 Jean*0122 C     ==================== calc. flux = upwind tracer * face-transport
                0123               flux(ix,iy,ir) = + wfac(ix,iy,ir) * intF
                0124
                0125               end if
                0126
                0127           end do
                0128
                0129           return
                0130
                0131 c     end subroutine GAD_PPM_FLX_R
                0132       end