/[MITgcm]/MITgcm/pkg/generic_advdiff/gad_pqm_hat_y.F
ViewVC logotype

Annotation of /MITgcm/pkg/generic_advdiff/gad_pqm_hat_y.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Sun Mar 13 01:44:03 2016 UTC (8 years, 1 month ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65v, checkpoint65w, checkpoint65u, HEAD
- from Darren: add PPM and PQM advection schemes (number 40-42 and 50-52)
  with 2 types of limiter (see: Engwirda & Kelley, submit. to JCP);
  Note (from Darren): unlimited PPM/PQM scheme (40 & 50) are just for
  testing and not for actual use.

1 jmc 1.1 C $Header: $
2     C $Name: $
3    
4     # include "GAD_OPTIONS.h"
5    
6     SUBROUTINE GAD_PQM_HAT_Y(bi,bj,kk,ix,
7     & method,mask,fbar,edge,
8     & ohat,fhat,myThid)
9     C |================================================================|
10     C | PQM_HAT_Y: reconstruct grid-cell PQM polynomials. |
11     C |================================================================|
12    
13     implicit none
14    
15     C =============================================== global variables
16     # include "SIZE.h"
17     # include "GRID.h"
18     # include "GAD.h"
19    
20     C ====================================================== arguments
21     C bi,bj :: tile indices.
22     C kk :: r-index.
23     C ix :: x-index.
24     C method :: advection scheme.
25     C mask :: row of cell-wise mask values.
26     C fbar :: row of cell-wise values.
27     C edge :: row of edge-wise values/slopes.
28     C - EDGE(1,:) = VALUE.
29     C - EDGE(2,:) = DF/DY.
30     C ohat :: row of oscl. coeff.
31     C - OHAT(1,:) = D^1F/DS^1.
32     C - OHAT(2,:) = D^2F/DS^2.
33     C fhat :: row of poly. coeff.
34     C - FHAT(:,I) = PQM coeff.
35     C myThid :: thread number.
36     C ================================================================
37     integer bi,bj,kk,ix
38     integer method
39     _RL mask(1-OLy:sNy+OLy)
40     _RL fbar(1-OLy:sNy+OLy)
41     _RL edge(1:2,
42     & 1-OLy:sNy+OLy)
43     _RL ohat(1:2,
44     & 1-OLy:sNy+OLy)
45     _RL fhat(1:5,
46     & 1-OLy:sNy+OLy)
47     integer myThid
48    
49     C ====================================================== variables
50     C ii,iy :: local x-indexing.
51     C ff00 :: centre-biased cell value.
52     C ffll,ffrr :: left-, and right-biased cell values.
53     C xhat :: local coord. scaling.
54     C fell,ferr :: left-, and right-biased edge values.
55     C dell,derr :: left-, and right-biased edge slopes.
56     C dfds :: linear slope estimates.
57     C uhat :: "NULL" limited poly. coeff.
58     C lhat :: "MONO" limited poly. coeff.
59     C scal :: "WENO" weights.
60     C fmag,fdel :: local perturbation indicators.
61     C ================================================================
62     integer ii,iy
63     _RL ff00
64     _RL ffll,ffrr
65     _RL yhat
66     _RL fell,ferr
67     _RL dell,derr
68     _RL dfds(-1:+1)
69     _RL uhat(+1:+5)
70     _RL lhat(+1:+5)
71     _RL scal(+1:+2)
72     _RL fmag,fdel
73     integer mono
74    
75     C ==================== reconstruct coeff. for grid-cell poynomials
76     do iy = 1-OLy+3, sNy+OLy-3
77    
78     if (mask(iy) .gt. 0. _d 0) then
79    
80     C =============================== scale to local grid-cell co-ords
81     yhat = dyF(ix,iy,bi,bj) * 0.5 _d 0
82    
83     C =============================== assemble cell mean + edge values
84     ff00 = fbar(iy+0)
85     ffll = ff00
86     & + mask(iy-1)*(fbar(iy-1)-ff00)
87     ffrr = ff00
88     & + mask(iy+1)*(fbar(iy+1)-ff00)
89    
90     fell = edge(+1,iy-0)
91     ferr = edge(+1,iy+1)
92    
93     dell = edge(+2,iy-0)
94     derr = edge(+2,iy+1)
95    
96     dell = dell * yhat
97     derr = derr * yhat
98    
99     c select case(method)
100     c case(ENUM_PQM_NULL_LIMIT)
101     if ( method.eq.ENUM_PQM_NULL_LIMIT ) then
102     C =============================== "NULL" limited grid-cell profile
103     CALL GAD_PQM_FUN_NULL ( ff00,
104     & fell,ferr,dell,derr,lhat,mono)
105    
106     c case(ENUM_PQM_MONO_LIMIT)
107     elseif ( method.eq.ENUM_PQM_MONO_LIMIT ) then
108     C =============================== "MONO" limited grid-cell profile
109     CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
110    
111     CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
112     & fell,ferr,dell,derr,dfds,lhat,
113     & mono)
114    
115     c case(ENUM_PQM_WENO_LIMIT)
116     elseif ( method.eq.ENUM_PQM_WENO_LIMIT ) then
117     C =============================== "WENO" limited grid-cell profile
118     CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
119    
120     CALL GAD_PQM_FUN_NULL ( ff00,
121     & fell,ferr,dell,derr,uhat,mono)
122    
123     CALL GAD_PQM_FUN_MONO ( ff00,ffll,ffrr,
124     & fell,ferr,dell,derr,dfds,lhat,
125     & mono)
126    
127     if ( mono .gt. 0) then
128    
129     C =============================== only apply WENO if it is worth it
130     fdel = abs(ffrr-ff00)+abs(ff00-ffll)
131     fmag = abs(ffll)+abs(ff00)+abs(ffrr)
132    
133     if (fdel .gt. 1. _d -6 * fmag) then
134    
135     C =============================== calc. WENO oscillation weighting
136     CALL GAD_OSC_MUL_Y(iy,+2,mask,
137     & ohat,scal)
138    
139     do ii = +1, +5
140     C =============================== blend limited/un-limited profile
141     lhat(ii) = scal(1) * uhat(ii)
142     & + scal(2) * lhat(ii)
143     end do
144    
145     end if
146    
147     end if
148    
149     c end select
150     endif
151    
152     do ii = +1, +5
153     C =============================== copy polynomial onto output data
154     fhat(ii,iy) = lhat(ii)
155     end do
156    
157     else
158    
159     do ii = +1, +5
160     fhat(ii,iy) = 0.0 _d 0
161     end do
162    
163     end if
164    
165     end do
166    
167     return
168    
169     c end subroutine GAD_PQM_HAT_Y
170     end

  ViewVC Help
Powered by ViewVC 1.1.22