pkg/generic_advdiff/gad_pqm_flx_r.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PQM_FLX_R(bi,bj,ix,iy,
                0004      &           delT,wvel,wfac,
                0005      &           fhat,flux,myThid)
                0006 C     |================================================================|
                0007 C     | PQM_FLX_R: evaluate PQM 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:5 ,
                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:5)
                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 (wfac(ix,iy,ir) .eq. 0. _d 0) then
                0067
                0068                   flux(ix,iy,ir)    = 0. _d 0
                0069
                0070               else
                0071
                0072               if (wfac(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               ivec(4) =(ss22 ** 4
                0088      &                - ss11 ** 4)*(1. _d 0 / 4. _d 0)
                0089               ivec(5) =(ss22 ** 5
                0090      &                - ss11 ** 5)*(1. _d 0 / 5. _d 0)
                0091
                0092               intF = ivec(1) * fhat(1,ir-1)
                0093      &             + ivec(2) * fhat(2,ir-1)
                0094      &             + ivec(3) * fhat(3,ir-1)
                0095      &             + ivec(4) * fhat(4,ir-1)
                0096      &             + ivec(5) * fhat(5,ir-1)
                0097
83ddf5a6c6 Mart*0098 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0099
                0100               else
                0101
                0102 C     ==================== integrate PQM profile over upwind cell IR+0
                0103               wCFL = wvel(ix,iy,ir)
                0104      &             * delT(ir-0)*recip_drF(ir-0)
                0105
                0106               ss11 = -1. _d 0 + 2. _d 0 * wCFL
                0107               ss22 = -1. _d 0
                0108
                0109 C     ==================== integrate profile over region swept by face
                0110               ivec(1) = ss22 - ss11
                0111               ivec(2) =(ss22 ** 2
                0112      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0113               ivec(3) =(ss22 ** 3
                0114      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0115               ivec(4) =(ss22 ** 4
                0116      &                - ss11 ** 4)*(1. _d 0 / 4. _d 0)
                0117               ivec(5) =(ss22 ** 5
                0118      &                - ss11 ** 5)*(1. _d 0 / 5. _d 0)
                0119
                0120               intF = ivec(1) * fhat(1,ir-0)
                0121      &             + ivec(2) * fhat(2,ir-0)
                0122      &             + ivec(3) * fhat(3,ir-0)
                0123      &             + ivec(4) * fhat(4,ir-0)
                0124      &             + ivec(5) * fhat(5,ir-0)
                0125
83ddf5a6c6 Mart*0126 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0127
                0128               end if
                0129
83ddf5a6c6 Mart*0130 c             intF = -0.5 _d 0 * intF / wCFL
                0131 C-    to avoid potential underflow:
                0132               intF = -0.5 _d 0 * intF/sign(max(abs(wCFL),1.d-20),wCFL)
                0133
8e4c181d69 Jean*0134 C     ==================== calc. flux = upwind tracer * face-transport
                0135               flux(ix,iy,ir) = + wfac(ix,iy,ir) * intF
                0136
                0137               end if
                0138
                0139           end do
                0140
                0141           return
                0142
                0143 c     end subroutine GAD_PQM_FLX_R
                0144       end