-
Notifications
You must be signed in to change notification settings - Fork 237
/
gad_dst2u1_impl_r.F
109 lines (97 loc) · 3.74 KB
/
gad_dst2u1_impl_r.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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
#include "GAD_OPTIONS.h"
CBOP
C !ROUTINE: GAD_DST2U1_IMPL_R
C !INTERFACE:
SUBROUTINE GAD_DST2U1_IMPL_R(
I bi,bj,k, iMin,iMax,jMin,jMax,
I advectionScheme, deltaTarg, rTrans, recip_hFac,
O a3d, b3d, c3d,
I myThid )
C !DESCRIPTION:
C Compute matrix element to solve vertical advection implicitly
C using DST 2nd.Order (=Lax-Wendroff) or 1rst Order Upwind scheme.
C Method:
C contribution of vertical transport at interface k is added
C to matrix lines k and k-1.
C !USES:
IMPLICIT NONE
C == Global variables ===
#include "SIZE.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GAD.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
C bi,bj :: tile indices
C k :: vertical level
C advectionScheme :: advection scheme to use: either 2nd Order DST
C or 1rst Order Upwind
C iMin,iMax :: computation domain
C jMin,jMax :: computation domain
C deltaTarg :: time step
C rTrans :: vertical volume transport
C recip_hFac :: inverse of cell open-depth factor
C a3d :: lower diagonal of the tridiagonal matrix
C b3d :: main diagonal of the tridiagonal matrix
C c3d :: upper diagonal of the tridiagonal matrix
C myThid :: thread number
INTEGER bi,bj,k
INTEGER iMin,iMax,jMin,jMax
INTEGER advectionScheme
_RL deltaTarg(Nr)
_RL rTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
_RS recip_hFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL a3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL b3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
_RL c3d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
INTEGER myThid
C == Local Variables ==
C i,j :: loop indices
C w_CFL :: Courant-Friedrich-Levy number
C rLimit :: centered (vs upwind) fraction
C rCenter :: centered contribution
C rUpwind :: upwind contribution
INTEGER i,j
c _RL w_CFL
_RL rLimit
_RL rCenter, rUpwind
_RL deltaTcfl
CEOP
rLimit = 0. _d 0
IF ( advectionScheme.EQ.ENUM_DST2 ) rLimit = 1. _d 0
C-- process interior interface only:
IF ( k.GT.1 .AND. k.LE.Nr ) THEN
C-- Add centered & upwind contributions
deltaTcfl = deltaTarg(k)
DO j=jMin,jMax
DO i=iMin,iMax
c w_CFL = deltaTcfl*ABS(rTrans(i,j))
c & *recip_rA(i,j,bi,bj)*recip_drC(k)
c & *recip_deepFac2F(k)*recip_rhoFacF(k)
rCenter = 0.5 _d 0 *rTrans(i,j)*recip_rA(i,j,bi,bj)*rkSign
rUpwind = ABS(rCenter)
& * ( 1. _d 0 - rLimit )
c & * ( 1. _d 0 - rLimit*( 1. _d 0 + w_CFL ) )
a3d(i,j,k) = a3d(i,j,k)
& - (rCenter+rUpwind)*deltaTarg(k)
& *recip_hFac(i,j,k)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
b3d(i,j,k) = b3d(i,j,k)
& - (rCenter-rUpwind)*deltaTarg(k)
& *recip_hFac(i,j,k)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
b3d(i,j,k-1) = b3d(i,j,k-1)
& + (rCenter+rUpwind)*deltaTarg(k-1)
& *recip_hFac(i,j,k-1)*recip_drF(k-1)
& *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
c3d(i,j,k-1) = c3d(i,j,k-1)
& + (rCenter-rUpwind)*deltaTarg(k-1)
& *recip_hFac(i,j,k-1)*recip_drF(k-1)
& *recip_deepFac2C(k-1)*recip_rhoFacC(k-1)
ENDDO
ENDDO
C-- process interior interface only: end
ENDIF
RETURN
END