forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
smooth_impldiff.F
179 lines (152 loc) · 4.28 KB
/
smooth_impldiff.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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
C $Header: /u/gcmpack/MITgcm/pkg/smooth/smooth_impldiff.F,v 1.3 2012/09/04 14:37:18 gforget Exp $
C $Name: $
#include "SMOOTH_OPTIONS.h"
SUBROUTINE SMOOTH_IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
I deltaTX, KappaRX, recip_hFac,
U gXNm1,
I myThid )
C *==========================================================*
C | SUBROUTINE smooth_impldiff
C | o Copy of model/src/impldiff
C | (simplified and with specified time step)
C *==========================================================*
C !USES:
IMPLICIT NONE
C == Global data ==
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine Arguments ==
INTEGER bi,bj,iMin,iMax,jMin,jMax
_RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
INTEGER i,j,k
_RL deltaTX
_RL gYnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
_RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
_RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
CEOP
c deltaTX = 1. _d 0
C-- Initialise
DO k=1,Nr
DO j=jMin,jMax
DO i=iMin,iMax
gYNm1(i,j,k,bi,bj) = 0. _d 0
ENDDO
ENDDO
ENDDO
C-- Old aLower
DO j=jMin,jMax
DO i=iMin,iMax
a(i,j,1) = 0. _d 0
ENDDO
ENDDO
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
a(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
& *KappaRX(i,j, k )*recip_drC( k )
& *deepFac2F(k)*rhoFacF(k)
IF (recip_hFac(i,j,k-1,bi,bj).EQ.0.) a(i,j,k)=0.
ENDDO
ENDDO
ENDDO
C-- Old aUpper
DO k=1,Nr-1
DO j=jMin,jMax
DO i=iMin,iMax
c(i,j,k) = -deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
& *recip_deepFac2C(k)*recip_rhoFacC(k)
& *KappaRX(i,j,k+1)*recip_drC(k+1)
& *deepFac2F(k+1)*rhoFacF(k+1)
IF (recip_hFac(i,j,k+1,bi,bj).EQ.0.) c(i,j,k)=0.
ENDDO
ENDDO
ENDDO
DO j=jMin,jMax
DO i=iMin,iMax
c(i,j,Nr) = 0. _d 0
ENDDO
ENDDO
C-- Old aCenter
DO k=1,Nr
DO j=jMin,jMax
DO i=iMin,iMax
b(i,j,k) = 1. _d 0 - c(i,j,k) - a(i,j,k)
ENDDO
ENDDO
ENDDO
C-- Old and new gam, bet are the same
DO k=1,Nr
DO j=jMin,jMax
DO i=iMin,iMax
bet(i,j,k) = 1. _d 0
gam(i,j,k) = 0. _d 0
ENDDO
ENDDO
ENDDO
C-- Only need do anything if Nr>1
IF (Nr.GT.1) THEN
k = 1
C-- Beginning of forward sweep (top level)
DO j=jMin,jMax
DO i=iMin,iMax
IF (b(i,j,1).NE.0.) bet(i,j,1) = 1. _d 0 / b(i,j,1)
ENDDO
ENDDO
ENDIF
C-- Middle of forward sweep
IF (Nr.GE.2) THEN
CADJ loop = sequential
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
gam(i,j,k) = c(i,j,k-1)*bet(i,j,k-1)
IF ( ( b(i,j,k) - a(i,j,k)*gam(i,j,k) ) .NE. 0.)
& bet(i,j,k) = 1. _d 0 / ( b(i,j,k) - a(i,j,k)*gam(i,j,k) )
ENDDO
ENDDO
ENDDO
ENDIF
DO j=jMin,jMax
DO i=iMin,iMax
gYNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j,1)
ENDDO
ENDDO
DO k=2,Nr
DO j=jMin,jMax
DO i=iMin,iMax
gYnm1(i,j,k,bi,bj) = bet(i,j,k)*
& (gXnm1(i,j,k,bi,bj) - a(i,j,k)*gYnm1(i,j,k-1,bi,bj))
ENDDO
ENDDO
ENDDO
C-- Backward sweep
CADJ loop = sequential
DO k=Nr-1,1,-1
DO j=jMin,jMax
DO i=iMin,iMax
gYnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
& -gam(i,j,k+1)*gYnm1(i,j,k+1,bi,bj)
ENDDO
ENDDO
ENDDO
DO k=1,Nr
DO j=jMin,jMax
DO i=iMin,iMax
gXnm1(i,j,k,bi,bj)=gYnm1(i,j,k,bi,bj)
ENDDO
ENDDO
ENDDO
RETURN
END