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

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

  ViewVC Help
Powered by ViewVC 1.1.22