pkg/generic_advdiff/gad_ppm_hat_x.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PPM_HAT_X(bi,bj,kk,iy,
                0004      &           method,mask,fbar,edge,
                0005      &           ohat,fhat,myThid)
                0006 C     |================================================================|
                0007 C     | PPM_HAT_X: reconstruct grid-cell PPM polynomials.              |
                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     ====================================================== arguments
                0018 C       bi,bj     :: tile indexing.
                0019 C       kk        :: r-indexing.
                0020 C       iy        :: y-indexing.
                0021 C       method    :: advection scheme.
                0022 C       mask      :: row of cell-wise masking values.
                0023 C       fbar      :: row of cell-wise values.
                0024 C       edge      :: row of edge-wise values.
                0025 C       ohat      :: row of oscl. coeff.
                0026 C                    - OHAT(1,:) = D^1F/DS^1.
                0027 C                    - OHAT(2,:) = D^2F/DS^2.
                0028 C       fhat      :: row of poly. coeff.
                0029 C                    - FHAT(:,I) = PPM coeff.
                0030 C       myThid    :: thread number.
                0031 C     ================================================================
                0032           integer bi,bj,kk,iy
                0033           integer method
                0034           _RL mask(1-OLx:sNx+OLx)
                0035           _RL fbar(1-OLx:sNx+OLx)
                0036           _RL edge(1-OLx:sNx+OLx)
                0037           _RL ohat(1:2,
                0038      &             1-OLx:sNx+OLx)
                0039           _RL fhat(1:3,
                0040      &             1-OLx:sNx+OLx)
                0041           integer myThid
                0042
                0043 C     ====================================================== variables
                0044 C       ii,ix     :: local x-indexing.
                0045 C       ff00      :: centre-biased cell mean.
                0046 C       ffll,ffrr :: left-, and right-biased cell values.
                0047 C       fell,ferr :: left-, and right-biased edge values.
                0048 C       dfds      :: linear slope estimates.
                0049 C       uhat      :: "NULL" limited poly. coeff.
                0050 C       lhat      :: "MONO" limited poly. coeff.
                0051 C       scal      :: "WENO" weights.
                0052 C       fmag,fdel :: local perturbation indicators.
                0053 C     ================================================================
                0054           integer ii,ix
                0055           _RL ff00,ffll,ffrr
                0056           _RL fell,ferr
                0057           _RL dfds(-1:+1)
                0058           _RL uhat(+1:+3)
                0059           _RL lhat(+1:+3)
                0060           _RL scal(+1:+2)
                0061           _RL fmag,fdel
                0062           integer  mono
                0063
                0064           do  ix = 1-OLx+2, sNx+OLx-2
                0065
                0066 C     =============================== assemble cell mean + edge values
                0067               ff00 = fbar(ix+0)
                0068               ffll = ff00
                0069      &             + mask(ix-1)*(fbar(ix-1)-ff00)
                0070               ffrr = ff00
                0071      &             + mask(ix+1)*(fbar(ix+1)-ff00)
                0072
                0073               fell = edge(ix-0)
                0074               ferr = edge(ix+1)
                0075
                0076 c             select case(method)
                0077 c             case(ENUM_PPM_NULL_LIMIT)
                0078               if ( method.eq.ENUM_PPM_NULL_LIMIT ) then
                0079 C     =============================== "NULL" limited grid-cell profile
                0080               CALL GAD_PPM_FUN_NULL ( ff00,fell,ferr,
                0081      &             lhat,mono)
                0082
                0083 c             case(ENUM_PPM_MONO_LIMIT)
                0084               elseif ( method.eq.ENUM_PPM_MONO_LIMIT ) then
                0085 C     =============================== "MONO" limited grid-cell profile
                0086               CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
                0087
                0088               CALL GAD_PPM_FUN_MONO ( ff00,ffll,ffrr,
                0089      &             fell,ferr,dfds,lhat,mono)
                0090
                0091 c             case(ENUM_PPM_WENO_LIMIT)
                0092               elseif ( method.eq.ENUM_PPM_WENO_LIMIT ) then
                0093 C     =============================== "WENO" limited grid-cell profile
                0094               CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
                0095
                0096               CALL GAD_PPM_FUN_NULL ( ff00,fell,ferr,
                0097      &             uhat,mono)
                0098
                0099               CALL GAD_PPM_FUN_MONO ( ff00,ffll,ffrr,
                0100      &             fell,ferr,dfds,lhat,mono)
                0101
                0102               if ( mono .gt. 0) then
                0103
                0104 C     =============================== only apply WENO if it is worth it
                0105               fdel = abs(ffrr-ff00)+abs(ff00-ffll)
                0106               fmag = abs(ffll)+abs(ff00)+abs(ffrr)
                0107
                0108               if (fdel .gt. 1. _d -6 * fmag) then
                0109
                0110 C     =============================== calc. WENO oscillation weighting
                0111                   CALL GAD_OSC_MUL_X(ix,+2,mask,
                0112      &                 ohat,scal)
                0113
                0114                   do  ii = +1, +3
                0115 C     =============================== blend limited/un-limited profile
                0116                   lhat(ii) = scal(1) * uhat(ii)
                0117      &                     + scal(2) * lhat(ii)
                0118                   end do
                0119
                0120               end if
                0121
                0122               end if
                0123
                0124 c             end select
                0125               endif
                0126
                0127               do  ii = +1, +3
                0128 C     =============================== copy polynomial onto output data
                0129                   fhat(ii,ix) = lhat(ii)
                0130               end do
                0131
                0132           end do
                0133
                0134           return
                0135
                0136 c     end subroutine GAD_PPM_HAT_X
                0137       end