pkg/generic_advdiff/gad_pqm_hat_y.F

8e4c181d69 Jean*0001 #     include "GAD_OPTIONS.h"
                0002
                0003       SUBROUTINE GAD_PQM_HAT_Y(bi,bj,kk,ix,
                0004      &           method,mask,fbar,edge,
                0005      &           ohat,fhat,myThid)
                0006 C     |================================================================|
                0007 C     | PQM_HAT_Y: reconstruct grid-cell PQM 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 indices.
                0019 C       kk        :: r-index.
                0020 C       ix        :: x-index.
                0021 C       method    :: advection scheme.
                0022 C       mask      :: row of cell-wise mask values.
                0023 C       fbar      :: row of cell-wise values.
                0024 C       edge      :: row of edge-wise values/slopes.
                0025 C                    - EDGE(1,:) = VALUE.
                0026 C                    - EDGE(2,:) = DF/DY.
                0027 C       ohat      :: row of oscl. coeff.
                0028 C                    - OHAT(1,:) = D^1F/DS^1.
                0029 C                    - OHAT(2,:) = D^2F/DS^2.
                0030 C       fhat      :: row of poly. coeff.
                0031 C                    - FHAT(:,I) = PQM coeff.
                0032 C       myThid    :: thread number.
                0033 C     ================================================================
                0034           integer bi,bj,kk,ix
                0035           integer method
                0036           _RL mask(1-OLy:sNy+OLy)
                0037           _RL fbar(1-OLy:sNy+OLy)
                0038           _RL edge(1:2,
                0039      &             1-OLy:sNy+OLy)
                0040           _RL ohat(1:2,
                0041      &             1-OLy:sNy+OLy)
                0042           _RL fhat(1:5,
                0043      &             1-OLy:sNy+OLy)
                0044           integer myThid
                0045
                0046 C     ====================================================== variables
                0047 C       ii,iy     :: local x-indexing.
                0048 C       ff00      :: centre-biased cell value.
                0049 C       ffll,ffrr :: left-, and right-biased cell values.
                0050 C       xhat      :: local coord. scaling.
                0051 C       fell,ferr :: left-, and right-biased edge values.
                0052 C       dell,derr :: left-, and right-biased edge slopes.
                0053 C       dfds      :: linear slope estimates.
                0054 C       uhat      :: "NULL" limited poly. coeff.
                0055 C       lhat      :: "MONO" limited poly. coeff.
                0056 C       scal      :: "WENO" weights.
                0057 C       fmag,fdel :: local perturbation indicators.
                0058 C     ================================================================
                0059           integer ii,iy
                0060           _RL ff00
                0061           _RL ffll,ffrr
                0062           _RL yhat
                0063           _RL fell,ferr
                0064           _RL dell,derr
                0065           _RL dfds(-1:+1)
                0066           _RL uhat(+1:+5)
                0067           _RL lhat(+1:+5)
                0068           _RL scal(+1:+2)
                0069           _RL fmag,fdel
                0070           integer  mono
                0071
                0072 C     ==================== reconstruct coeff. for grid-cell poynomials
                0073           do  iy = 1-OLy+3, sNy+OLy-3
                0074
                0075              if (mask(iy) .gt. 0. _d 0) then
                0076
                0077 C     =============================== scale to local grid-cell co-ords
                0078               yhat = dyF(ix,iy,bi,bj) * 0.5 _d 0
                0079
                0080 C     =============================== assemble cell mean + edge values
                0081               ff00 = fbar(iy+0)
                0082               ffll = ff00
                0083      &         + mask(iy-1)*(fbar(iy-1)-ff00)
                0084               ffrr = ff00
                0085      &         + mask(iy+1)*(fbar(iy+1)-ff00)
                0086
                0087               fell = edge(+1,iy-0)
                0088               ferr = edge(+1,iy+1)
                0089
                0090               dell = edge(+2,iy-0)
                0091               derr = edge(+2,iy+1)
                0092
                0093               dell = dell * yhat
                0094               derr = derr * yhat
                0095
                0096 c             select case(method)
                0097 c             case(ENUM_PQM_NULL_LIMIT)
                0098               if ( method.eq.ENUM_PQM_NULL_LIMIT ) then
                0099 C     =============================== "NULL" limited grid-cell profile
                0100               CALL GAD_PQM_FUN_NULL ( ff00,
                0101      &             fell,ferr,dell,derr,lhat,mono)
                0102
                0103 c             case(ENUM_PQM_MONO_LIMIT)
                0104               elseif ( method.eq.ENUM_PQM_MONO_LIMIT ) then
                0105 C     =============================== "MONO" limited grid-cell profile
                0106               CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
                0107
                0108               CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
                0109      &             fell,ferr,dell,derr,dfds,lhat,
                0110      &             mono)
                0111
                0112 c             case(ENUM_PQM_WENO_LIMIT)
                0113               elseif ( method.eq.ENUM_PQM_WENO_LIMIT ) then
                0114 C     =============================== "WENO" limited grid-cell profile
                0115               CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
                0116
                0117               CALL GAD_PQM_FUN_NULL ( ff00,
                0118      &             fell,ferr,dell,derr,uhat,mono)
                0119
                0120               CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
                0121      &             fell,ferr,dell,derr,dfds,lhat,
                0122      &             mono)
                0123
                0124               if ( mono .gt. 0) then
                0125
                0126 C     =============================== only apply WENO if it is worth it
                0127               fdel = abs(ffrr-ff00)+abs(ff00-ffll)
                0128               fmag = abs(ffll)+abs(ff00)+abs(ffrr)
                0129
                0130               if (fdel .gt. 1. _d -6 * fmag) then
                0131
                0132 C     =============================== calc. WENO oscillation weighting
                0133                   CALL GAD_OSC_MUL_Y(iy,+2,mask,
                0134      &                 ohat,scal)
                0135
                0136                   do  ii = +1, +5
                0137 C     =============================== blend limited/un-limited profile
                0138                   lhat(ii) = scal(1) * uhat(ii)
                0139      &                     + scal(2) * lhat(ii)
                0140                   end do
                0141
                0142               end if
                0143
                0144               end if
                0145
                0146 c             end select
                0147               endif
                0148
                0149               do  ii = +1, +5
                0150 C     =============================== copy polynomial onto output data
                0151                   fhat(ii,iy) = lhat(ii)
                0152               end do
                0153
                0154               else
                0155
                0156               do  ii = +1, +5
                0157                   fhat(ii,iy) = 0.0 _d 0
                0158               end do
                0159
                0160               end if
                0161
                0162           end do
                0163
                0164           return
                0165
                0166 c     end subroutine GAD_PQM_HAT_Y
                0167       end