pkg/generic_advdiff/gad_pqm_flx_x.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PQM_FLX_X(bi,bj,kk,iy,
                0004      &           calc_CFL,delT,uvel,
                0005      &           ufac,fhat,flux,myThid)
                0006 C     |================================================================|
                0007 C     | PQM_FLX_X: 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       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:5,
                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:5)
                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+4, sNx+OLx-3
                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 PQM 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               ivec(4) =(ss22 ** 4
                0092      &                - ss11 ** 4)*(1. _d 0 / 4. _d 0)
                0093               ivec(5) =(ss22 ** 5
                0094      &                - ss11 ** 5)*(1. _d 0 / 5. _d 0)
                0095
                0096               intF = ivec(1) * fhat(1,ix-1)
                0097      &             + ivec(2) * fhat(2,ix-1)
                0098      &             + ivec(3) * fhat(3,ix-1)
                0099      &             + ivec(4) * fhat(4,ix-1)
                0100      &             + ivec(5) * fhat(5,ix-1)
                0101
83ddf5a6c6 Mart*0102 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0103
                0104               else
                0105
                0106 C     ==================== integrate PQM profile over upwind cell IX+0
                0107               if ( calc_CFL ) then
                0108               uCFL = uvel(ix,iy) * delT
                0109      &             * recip_dxF(ix-0,iy,bi,bj)
                0110      &             * recip_deepFacC(kk)
                0111               else
                0112               uCFL = uvel(ix,iy)
                0113               end if
                0114
                0115               ss11 = -1. _d 0 - 2. _d 0 * uCFL
                0116               ss22 = -1. _d 0
                0117
                0118 C     ==================== integrate profile over region swept by face
                0119               ivec(1) = ss22 - ss11
                0120               ivec(2) =(ss22 ** 2
                0121      &                - ss11 ** 2)*(1. _d 0 / 2. _d 0)
                0122               ivec(3) =(ss22 ** 3
                0123      &                - ss11 ** 3)*(1. _d 0 / 3. _d 0)
                0124               ivec(4) =(ss22 ** 4
                0125      &                - ss11 ** 4)*(1. _d 0 / 4. _d 0)
                0126               ivec(5) =(ss22 ** 5
                0127      &                - ss11 ** 5)*(1. _d 0 / 5. _d 0)
                0128
                0129               intF = ivec(1) * fhat(1,ix-0)
                0130      &             + ivec(2) * fhat(2,ix-0)
                0131      &             + ivec(3) * fhat(3,ix-0)
                0132      &             + ivec(4) * fhat(4,ix-0)
                0133      &             + ivec(5) * fhat(5,ix-0)
                0134
83ddf5a6c6 Mart*0135 CML              intF = intF / (ss22 - ss11)
8e4c181d69 Jean*0136
                0137               end if
                0138
83ddf5a6c6 Mart*0139 c             intF = 0.5 _d 0 * intF / uCFL
                0140 C-    to avoid potential underflow:
                0141               intF = 0.5 _d 0 * intF / sign(max(abs(uCFL),1.d-20),uCFL)
                0142
8e4c181d69 Jean*0143 C     ==================== calc. flux = upwind tracer * face-transport
                0144               flux(ix,iy) = + ufac(ix,iy) * intF
                0145
                0146               end if
                0147
                0148           end do
                0149
                0150           return
                0151
                0152 c     end subroutine GAD_PQM_FLX_X
                0153       end