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

Contents of /MITgcm/pkg/generic_advdiff/gad_ppm_hat_r.F

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


Revision 1.1 - (show 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 C $Header: $
2 C $Name: $
3
4 # include "GAD_OPTIONS.h"
5
6 SUBROUTINE GAD_PPM_HAT_R(bi,bj,ix,iy,
7 & method,mask,fbar,edge,
8 & ohat,fhat,myThid)
9 C |================================================================|
10 C | PPM_HAT_R: reconstruct grid-cell PPM 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 indexing.
22 C ix,iy :: x-,y-indexing.
23 C method :: advection scheme.
24 C mask :: row of cell-wise masking values.
25 C fbar :: row of cell-wise values.
26 C edge :: row of edge-wise values.
27 C ohat :: row of oscl. coeff.
28 C - OHAT(1,:) = D^1F/DS^1.
29 C - OHAT(2,:) = D^2F/DS^2.
30 C fhat :: row of poly. coeff.
31 C - FHAT(:,I) = PPM coeff.
32 C myThid :: thread number.
33 C ================================================================
34 integer bi,bj,ix,iy
35 integer method
36 _RL mask(1-3:Nr+3)
37 _RL fbar(1-3:Nr+3)
38 _RL edge(1-0:Nr+1)
39 _RL ohat(1:2,
40 & 1-3:Nr+3)
41 _RL fhat(1:3,
42 & 1-0:Nr+0)
43 integer myThid
44
45 C ====================================================== variables
46 C ii,ir :: local r-indexing.
47 C ff00 :: centre-biased cell mean.
48 C ffll,ffrr :: left-, and right-biased cell values.
49 C fell,ferr :: left-, and right-biased edge values.
50 C dfds :: linear slope estimates.
51 C uhat :: "NULL" limited poly. coeff.
52 C lhat :: "MONO" limited poly. coeff.
53 C scal :: "WENO" weights.
54 C fmag,fdel :: local perturbation indicators.
55 C ================================================================
56 integer ii,ir
57 _RL ff00,ffll,ffrr
58 _RL fell,ferr
59 _RL dfds(-1:+1)
60 _RL uhat(+1:+5)
61 _RL lhat(+1:+5)
62 _RL scal(+1:+2)
63 _RL fmag,fdel
64 integer mono
65
66 do ir = +1, Nr
67
68 C =============================== assemble cell mean + edge values
69 ff00 = fbar(ir+0)
70 ffll = ff00
71 & + mask(ir-1)*(fbar(ir-1)-ff00)
72 ffrr = ff00
73 & + mask(ir+1)*(fbar(ir+1)-ff00)
74
75 fell = edge(ir-0)
76 ferr = edge(ir+1)
77
78 c select case(method)
79 c case(ENUM_PPM_NULL_LIMIT)
80 if ( method.eq.ENUM_PPM_NULL_LIMIT ) then
81 C =============================== "NULL" limited grid-cell profile
82 CALL GAD_PPM_FUN_NULL ( ff00,fell,ferr,
83 & lhat,mono)
84
85 c case(ENUM_PPM_MONO_LIMIT)
86 elseif ( method.eq.ENUM_PPM_MONO_LIMIT ) then
87 C =============================== "MONO" limited grid-cell profile
88 CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
89
90 CALL GAD_PPM_FUN_MONO ( ff00,ffll,ffrr,
91 & fell,ferr,dfds,lhat,mono)
92
93 c case(ENUM_PPM_WENO_LIMIT)
94 elseif ( method.eq.ENUM_PPM_WENO_LIMIT ) then
95 C =============================== "WENO" limited grid-cell profile
96 CALL GAD_PLM_FUN_U(ffll,ff00,ffrr,dfds)
97
98 CALL GAD_PPM_FUN_NULL ( ff00,fell,ferr,
99 & uhat,mono)
100
101 CALL GAD_PPM_FUN_MONO ( ff00,ffll,ffrr,
102 & fell,ferr,dfds,lhat,mono)
103
104 if ( mono .gt. 0) then
105
106 C =============================== only apply WENO if it is worth it
107 fdel = abs(ffrr-ff00)+abs(ff00-ffll)
108 fmag = abs(ffll)+abs(ff00)+abs(ffrr)
109
110 if (fdel .gt. 1. _d -6 * fmag) then
111
112 C =============================== calc. WENO oscillation weighting
113 CALL GAD_OSC_MUL_R(ir,+2,mask,
114 & ohat,scal)
115
116 do ii = +1, +3
117 C =============================== blend limited/un-limited profile
118 lhat(ii) = scal(1) * uhat(ii)
119 & + scal(2) * lhat(ii)
120 end do
121
122 end if
123
124 end if
125
126 c end select
127 endif
128
129 do ii = +1, +3
130 C =============================== copy polynomial onto output data
131 fhat(ii,ir) = lhat(ii)
132 end do
133
134 end do
135
136 return
137
138 c end subroutine GAD_PPM_HAT_R
139 end

  ViewVC Help
Powered by ViewVC 1.1.22