forked from altMITgcm/MITgcm
-
Notifications
You must be signed in to change notification settings - Fork 0
/
smooth_diff2d.F
202 lines (170 loc) · 5.69 KB
/
smooth_diff2d.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
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
C $Header: /u/gcmpack/MITgcm/pkg/smooth/smooth_diff2d.F,v 1.2 2010/04/13 04:43:45 gforget Exp $
C $Name: $
#include "SMOOTH_OPTIONS.h"
subroutine smooth_diff2D (
U fld_in,mask_in,nbt_in,mythid)
C *==========================================================*
C | SUBROUTINE smooth_diff2D
C | o Routine that smoothes a 2D field, using diffusion
C *==========================================================*
IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"
#include "PARAMS.h"
#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc.h"
#include "tamc_keys.h"
#endif /* ALLOW_AUTODIFF_TAMC */
#include "SMOOTH.h"
integer i,j,k, bi, bj
integer itlo,ithi
integer jtlo,jthi
integer myThid,myIter(nSx,nSy),key_in
_RL smooth2Dmask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL mask_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nR,nSx,nSy)
_RL fld_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL gt_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
_RL gtm1_in(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
integer nbt_in
integer iloop, ilev_1, ilev_2, ilev_3
integer max_lev2, max_lev3
_RL ab15,ab05
_RL gt_tmp
character*( 80) fnamegeneric
c for now: useless, because level 3 is recomputed anyway
c but : if level3 was computed during the fwd loop by callung
c mdsmooth_diff3D (assumes that it would be called
c directly by the_main_loop) then I would need to pass key_in
c as a parameter, with different values for T, S, ...
c in order not to overwrite the same tape
key_in=0
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
DO bj = jtlo,jthi
DO bi = itlo,ithi
DO j = 1,sNy
DO i = 1,sNx
gt_in(i,j,bi,bj) = 0. _d 0
gtm1_in(i,j,bi,bj) = 0. _d 0
smooth2Dmask(i,j,bi,bj) = mask_in(i,j,1,bi,bj)
ENDDO
ENDDO
ENDDO
ENDDO
_EXCH_XY_RL ( fld_in , myThid )
_EXCH_XY_RL ( gt_in , myThid )
_EXCH_XY_RL ( gtm1_in , myThid )
_EXCH_XY_RL ( smooth2Dmask , myThid )
#ifdef ALLOW_TAMC_CHECKPOINTING
c checkpointing:
max_lev3=nbt_in/(nchklev_1*nchklev_2)+1
max_lev2=nbt_in/nchklev_1+1
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_smooth2D_lev3 = USER
#endif /* ALLOW_AUTODIFF_TAMC */
do ilev_3 = 1,nchklev_3
if(ilev_3.le.max_lev3) then
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fld_in = tape_smooth2D_lev3 ,
CADJ & key = key_in*max_lev3 + ilev_3
CADJ STORE gTm1_in = tape_smooth2D_lev3 ,
CADJ & key = key_in*max_lev3 + ilev_3
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_smooth2D_lev2 = USER
#endif /* ALLOW_AUTODIFF_TAMC */
do ilev_2 = 1,nchklev_2
if(ilev_2.le.max_lev2) then
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE fld_in = tape_smooth2D_lev2 ,
CADJ & key = key_in*nchklev_2 + ilev_2
CADJ STORE gTm1_in = tape_smooth2D_lev2 ,
CADJ & key = key_in*nchklev_2 + ilev_2
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_AUTODIFF_TAMC
CADJ INIT tape_smooth2D_lev1 = COMMON,
CADJ & nchklev_1*nsx*nsy*nthreads_chkpt
#endif /* ALLOW_AUTODIFF_TAMC */
do ilev_1 = 1,nchklev_1
iloop = (ilev_2 - 1)*nchklev_1 + ilev_1
& + (ilev_3 - 1)*nchklev_2*nchklev_1
if ( iloop .le. nbt_in ) then
#ifdef ALLOW_AUTODIFF_TAMC
c needed?? CADJ STORE fld_in = tape_smooth2D_lev1 ,
c CADJ & key = key_in*nchklev_1 + ilev_1
CADJ STORE gtm1_in = tape_smooth2D_lev1 ,
CADJ & key = key_in*nchklev_1 + ilev_1
#endif /* ALLOW_AUTODIFF_TAMC */
#else /* ALLOW_TAMC_CHECKPOINTING */
do iloop=1,nbt_in
#endif
DO bj = jtlo,jthi
DO bi = itlo,ithi
DO j = 1,sNy
DO i = 1,sNx
gt_in(i,j,bi,bj)=0.
if (smooth2Dmask(i,j,bi,bj).NE.0.) then
gt_in(i,j,bi,bj)=gt_in(i,j,bi,bj)+
& smooth2D_Kux(i,j,bi,bj)*dyG(i,j,bi,bj)*
& smooth2Dmask(i,j,bi,bj)*smooth2Dmask(i-1,j,bi,bj)*
& (fld_in(i,j,bi,bj)-fld_in(i-1,j,bi,bj))*recip_dxC(i,j,bi,bj)
gt_in(i,j,bi,bj)=gt_in(i,j,bi,bj)+
& smooth2D_Kux(i+1,j,bi,bj)*dyG(i+1,j,bi,bj)*
& smooth2Dmask(i,j,bi,bj)*smooth2Dmask(i+1,j,bi,bj)*
& (fld_in(i,j,bi,bj)-fld_in(i+1,j,bi,bj))*recip_dxC(i+1,j,bi,bj)
gt_in(i,j,bi,bj)=gt_in(i,j,bi,bj)+
& smooth2D_Kvy(i,j,bi,bj)*dxG(i,j,bi,bj)*
& smooth2Dmask(i,j,bi,bj)*smooth2Dmask(i,j-1,bi,bj)*
& (fld_in(i,j,bi,bj)-fld_in(i,j-1,bi,bj))*recip_dyC(i,j,bi,bj)
gt_in(i,j,bi,bj)=gt_in(i,j,bi,bj)+
& smooth2D_Kvy(i,j+1,bi,bj)*dxG(i,j+1,bi,bj)*
& smooth2Dmask(i,j,bi,bj)*smooth2Dmask(i,j+1,bi,bj)*
& (fld_in(i,j,bi,bj)-fld_in(i,j+1,bi,bj))*recip_dyC(i,j+1,bi,bj)
endif
ENDDO
ENDDO
ENDDO
ENDDO
do bj = jtlo,jthi
do bi = itlo,ithi
c Adams-Bashforth timestepping
myIter(bi,bj)=iloop-1
IF ( myIter(bi,bj).EQ.0 ) THEN
ab15=1.0
ab05=0.0
ELSE
ab15=1.5+abEps
ab05=-(0.5+abEps)
ENDIF
DO j = 1,sNy
DO i = 1,sNx
c Compute effective G-term with Adams-Bashforth
gt_tmp = ab15*gt_in(i,j,bi,bj)
& + ab05*gtm1_in(i,j,bi,bj)
gtm1_in(i,j,bi,bj) = gt_in(i,j,bi,bj)
gt_in(i,j,bi,bj) = gt_tmp
c time step:
fld_in(i,j,bi,bj)=fld_in(i,j,bi,bj)
& -gt_in(i,j,bi,bj)*recip_rA(i,j,bi,bj)*smooth2DdelTime
gt_in(i,j,bi,bj)=0
ENDDO
ENDDO
ENDDO
ENDDO
_EXCH_XY_RL ( gt_in , myThid )
_EXCH_XY_RL ( fld_in , myThid )
_EXCH_XY_RL ( gtm1_in , myThid )
#ifdef ALLOW_TAMC_CHECKPOINTING
endif
enddo
endif
enddo
endif
enddo
#else /* ALLOW_TAMC_CHECKPOINTING */
enddo
#endif
end