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

Annotation of /MITgcm/pkg/generic_advdiff/gad_pqm_flx_r.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_FLX_R(bi,bj,ix,iy,
7     & delT,wvel,wfac,
8     & fhat,flux,myThid)
9     C |================================================================|
10     C | PQM_FLX_R: evaluate PQM flux on grid-cell edges. |
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 ================================================================
21     C bi,bj :: tile indexing.
22     C ix :: x-index.
23     C iy :: y-index.
24     C delT :: time-step.
25     C wvel :: vel.-comp in r-direction.
26     C wfac :: vel.-flux in r-direction.
27     C fhat :: row of poly. coeff.
28     C flux :: adv.-flux in r-direction.
29     C myThid :: thread number.
30     C ================================================================
31     integer bi,bj,ix,iy
32     _RL delT(1:Nr)
33     _RL wvel(1-OLx:sNx+OLx,
34     & 1-OLy:sNy+OLy,
35     & 1:Nr)
36     _RL wfac(1-OLx:sNx+OLx,
37     & 1-OLy:sNy+OLy,
38     & 1:Nr)
39     _RL fhat(1:5 ,
40     & 1:Nr)
41     _RL flux(1-OLx:sNx+OLx,
42     & 1-OLy:sNy+OLy,
43     & 1:Nr)
44     integer myThid
45    
46     C ================================================================
47     C ir :: r-indexing.
48     C wCFL :: CFL number.
49     C intF :: upwind tracer edge-value.
50     C ss11,ss22 :: int. endpoints.
51     C ivec :: int. basis vec.
52     C ================================================================
53     integer ir
54     _RL wCFL,intF
55     _RL ss11,ss22
56     _RL ivec(1:5)
57    
58     C ================================================================
59     C (1): calc. "departure-points" for each grid-cell edge by int-
60     C egrating edge position backward in time over one single
61     C time-step. This is a "single-cell" implementation: requ-
62     C ires CFL <= 1.0.
63     C (2): calc. flux as the integral of the upwind grid-cell poly-
64     C nomial over the deformation interval found in (1).
65     C ================================================================
66    
67     do ir = +2, Nr
68    
69     if (wfac(ix,iy,ir) .eq. 0. _d 0) then
70    
71     flux(ix,iy,ir) = 0. _d 0
72    
73     else
74    
75     if (wfac(ix,iy,ir) .lt. 0. _d 0) then
76    
77     C ==================== integrate PQM profile over upwind cell IR-1
78     wCFL = wvel(ix,iy,ir)
79     & * delT(ir-1)*recip_drF(ir-1)
80    
81     ss11 = +1. _d 0 + 2. _d 0 * wCFL
82     ss22 = +1. _d 0
83    
84     C ==================== integrate profile over region swept by face
85     ivec(1) = ss22 - ss11
86     ivec(2) =(ss22 ** 2
87     & - ss11 ** 2)*(1. _d 0 / 2. _d 0)
88     ivec(3) =(ss22 ** 3
89     & - ss11 ** 3)*(1. _d 0 / 3. _d 0)
90     ivec(4) =(ss22 ** 4
91     & - ss11 ** 4)*(1. _d 0 / 4. _d 0)
92     ivec(5) =(ss22 ** 5
93     & - ss11 ** 5)*(1. _d 0 / 5. _d 0)
94    
95     intF = ivec(1) * fhat(1,ir-1)
96     & + ivec(2) * fhat(2,ir-1)
97     & + ivec(3) * fhat(3,ir-1)
98     & + ivec(4) * fhat(4,ir-1)
99     & + ivec(5) * fhat(5,ir-1)
100    
101     intF = intF / (ss22 - ss11)
102    
103     else
104    
105     C ==================== integrate PQM profile over upwind cell IR+0
106     wCFL = wvel(ix,iy,ir)
107     & * delT(ir-0)*recip_drF(ir-0)
108    
109     ss11 = -1. _d 0 + 2. _d 0 * wCFL
110     ss22 = -1. _d 0
111    
112     C ==================== integrate profile over region swept by face
113     ivec(1) = ss22 - ss11
114     ivec(2) =(ss22 ** 2
115     & - ss11 ** 2)*(1. _d 0 / 2. _d 0)
116     ivec(3) =(ss22 ** 3
117     & - ss11 ** 3)*(1. _d 0 / 3. _d 0)
118     ivec(4) =(ss22 ** 4
119     & - ss11 ** 4)*(1. _d 0 / 4. _d 0)
120     ivec(5) =(ss22 ** 5
121     & - ss11 ** 5)*(1. _d 0 / 5. _d 0)
122    
123     intF = ivec(1) * fhat(1,ir-0)
124     & + ivec(2) * fhat(2,ir-0)
125     & + ivec(3) * fhat(3,ir-0)
126     & + ivec(4) * fhat(4,ir-0)
127     & + ivec(5) * fhat(5,ir-0)
128    
129     intF = intF / (ss22 - ss11)
130    
131     end if
132    
133     C ==================== calc. flux = upwind tracer * face-transport
134     flux(ix,iy,ir) = + wfac(ix,iy,ir) * intF
135    
136     end if
137    
138     end do
139    
140     return
141    
142     c end subroutine GAD_PQM_FLX_R
143     end
  ViewVC Help
Powered by ViewVC 1.1.22