forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
gad_pqm_p5e_x.F
82 lines (64 loc) · 2.67 KB
/
gad_pqm_p5e_x.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
C $Header: /u/gcmpack/MITgcm/pkg/generic_advdiff/gad_pqm_p5e_x.F,v 1.1 2016/03/13 01:44:03 jmc Exp $
C $Name: $
# include "GAD_OPTIONS.h"
SUBROUTINE GAD_PQM_P5E_X(bi,bj,kk,iy,
& mask,fbar,edge,myThid)
C |================================================================|
C | PQM_P5E_X: approximate edge values with degree-5 polynomials. |
C | Fixed grid-spacing variant in X. |
C |================================================================|
implicit none
C =============================================== global variables
# include "SIZE.h"
# include "GRID.h"
# include "GAD.h"
C ====================================================== arguments
integer bi,bj,kk,iy
_RL mask(1-OLx:sNx+OLx)
_RL fbar(1-OLx:sNx+OLx)
_RL edge(1:2,
& 1-OLx:sNx+OLx)
integer myThid
C ====================================================== variables
integer ix
_RL mloc(-3:+2)
_RL floc(-3:+2)
_RL ftmp
do ix = 1-OLx+3, sNx+OLx-2
C ================ mask local stencil: expand from centre outwards
mloc(-1) = mask(ix-1)
mloc(+0) = mask(ix+0)
floc(-1) = fbar(ix+0)
& + mloc(-1)*(fbar(ix-1)-fbar(ix+0))
floc(+0) = fbar(ix-1)
& + mloc(+0)*(fbar(ix+0)-fbar(ix-1))
mloc(-2) = mask(ix-2) * mloc(-1)
mloc(-3) = mask(ix-3) * mloc(-2)
ftmp = 2. _d 0 * floc(-1) - floc(+0)
floc(-2) = ftmp
& + mloc(-2)*(fbar(ix-2)-ftmp)
ftmp = 2. _d 0 * floc(-2) - floc(-1)
floc(-3) = ftmp
& + mloc(-3)*(fbar(ix-3)-ftmp)
mloc(+1) = mask(ix+1) * mloc(+0)
mloc(+2) = mask(ix+2) * mloc(+1)
ftmp = 2. _d 0 * floc(+0) - floc(-1)
floc(+1) = ftmp
& + mloc(+1)*(fbar(ix+1)-ftmp)
ftmp = 2. _d 0 * floc(+1) - floc(+0)
floc(+2) = ftmp
& + mloc(+2)*(fbar(ix+2)-ftmp)
C ================ centred, 5th-order interpolation for edge value
edge(1,ix) =
& +( 1. _d 0/60. _d 0)*(floc(-3)+floc(+2))
& -( 8. _d 0/60. _d 0)*(floc(-2)+floc(+1))
& +(37. _d 0/60. _d 0)*(floc(-1)+floc(+0))
edge(2,ix) = (
& -( 1. _d 0/90. _d 0)*(floc(-3)-floc(+2))
& +( 5. _d 0/36. _d 0)*(floc(-2)-floc(+1))
& -(49. _d 0/36. _d 0)*(floc(-1)-floc(+0))
& ) * recip_dxC(ix,iy,bi,bj)
end do
return
c end subroutine GAD_PQM_P5E_X
end