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

Contents 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 - (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_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