/[MITgcm]/MITgcm/pkg/flt/flt_runga4.F
ViewVC logotype

Annotation of /MITgcm/pkg/flt/flt_runga4.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (hide annotations) (download)
Thu Dec 22 19:04:45 2011 UTC (12 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint64, checkpoint65, checkpoint66g, checkpoint66f, checkpoint66e, checkpoint66d, checkpoint66c, checkpoint66b, checkpoint66a, checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, checkpoint66i, checkpoint66h, checkpoint65z, checkpoint65x, checkpoint65y, checkpoint65r, checkpoint65s, checkpoint65p, checkpoint65q, checkpoint65v, checkpoint65w, checkpoint65t, checkpoint65u, checkpoint65j, checkpoint65k, checkpoint65h, checkpoint65i, checkpoint65n, checkpoint65o, checkpoint65l, checkpoint65m, checkpoint65b, checkpoint65c, checkpoint65a, checkpoint65f, checkpoint65g, checkpoint65d, checkpoint65e, HEAD
Changes since 1.3: +3 -4 lines
remove/avoid un-used variables

1 jmc 1.4 C $Header: /u/gcmpack/MITgcm/pkg/flt/flt_runga4.F,v 1.3 2011/08/31 21:32:46 jmc Exp $
2 jahn 1.2 C $Name: $
3 jahn 1.1
4     #include "FLT_OPTIONS.h"
5 jmc 1.3 #undef _USE_INTEGERS
6 jahn 1.1
7     SUBROUTINE FLT_RUNGA4 (
8     I myTime, myIter, myThid )
9    
10     C ==================================================================
11     C SUBROUTINE FLT_RUNGA4
12     C ==================================================================
13     C o This routine steps floats forward with second order Runge-Kutta
14     C
15     C started: Arne Biastoch
16     C
17     C changed: 2004.06.10 Antti Westerlund (antti.westerlund@helsinki.fi)
18     C and Sergio Jaramillo (sju@eos.ubc.ca)
19     C ==================================================================
20    
21     C !USES:
22     IMPLICIT NONE
23    
24     C == global variables ==
25     #include "SIZE.h"
26     #include "EEPARAMS.h"
27     #include "PARAMS.h"
28     #include "GRID.h"
29     #include "DYNVARS.h"
30 jahn 1.2 #include "FLT_SIZE.h"
31 jahn 1.1 #include "FLT.h"
32    
33     C == routine arguments ==
34     _RL myTime
35     INTEGER myIter, myThid
36    
37     C == Functions ==
38 jmc 1.3 #ifdef USE_FLT_ALT_NOISE
39     Real*8 PORT_RAND_NORM
40     EXTERNAL PORT_RAND_NORM
41     #else
42     Real*8 PORT_RAND
43     EXTERNAL PORT_RAND
44     #ifdef _USE_INTEGERS
45     INTEGER seed
46     #else
47     Real*8 seed
48     #endif
49     #endif /* USE_FLT_ALT_NOISE */
50 jahn 1.1
51     C == local variables ==
52     INTEGER bi, bj
53     INTEGER ip
54     INTEGER ic, jc, kc, iG, jG
55     _RL u1, v1, u2, v2, u3, v3, u4, v4
56     #ifdef ALLOW_3D_FLT
57     _RL w1, w2, w3, w4, ktz1, ktz2, ktz3, kz, scalez
58     _RL kzlo, kzhi
59     #endif
60     _RL ix, jy, itx1, jty1, itx2, jty2, itx3, jty3
61     _RL scalex, scaley
62 jmc 1.3
63     C == end of interface ==
64    
65     #ifndef USE_FLT_ALT_NOISE
66 jahn 1.1 #ifdef _USE_INTEGERS
67     seed = -1
68     #else
69     seed = -1.d0
70     #endif
71 jmc 1.3 #endif /* ndef USE_FLT_ALT_NOISE */
72 jahn 1.1
73     #ifdef ALLOW_3D_FLT
74     kzlo = 0.5 _d 0
75     kzhi = 0.5 _d 0 + DFLOAT(Nr)
76     #endif
77    
78     DO bj=myByLo(myThid),myByHi(myThid)
79     DO bi=myBxLo(myThid),myBxHi(myThid)
80     DO ip=1,npart_tile(bi,bj)
81    
82     C If float has died move to level 0
83     IF ( tend(ip,bi,bj).NE.-1. .AND. myTime.GT.tend(ip,bi,bj)
84     & ) THEN
85     kpart(ip,bi,bj) = 0.
86     ELSE
87     C Start integration between tstart and tend (individual for each float)
88     IF ( (tstart(ip,bi,bj).EQ.-1..OR.myTime.GE.tstart(ip,bi,bj))
89     & .AND.( tend(ip,bi,bj).EQ.-1..OR.myTime.LE. tend(ip,bi,bj))
90     & .AND.( iup(ip,bi,bj).NE.-3.)
91     & ) THEN
92    
93     ix = ipart(ip,bi,bj)
94     jy = jpart(ip,bi,bj)
95     ic=NINT(ix)
96     jc=NINT(jy)
97     kc=NINT(kpart(ip,bi,bj))
98    
99     scalex=recip_dxF(ic,jc,bi,bj)
100     scaley=recip_dyF(ic,jc,bi,bj)
101 jmc 1.3 iG = myXGlobalLo + (bi-1)*sNx + ic-1
102 jahn 1.1 jG = myYGlobalLo + (bj-1)*sNy + jc-1
103    
104     C First step
105    
106     #ifdef ALLOW_3D_FLT
107     IF (iup(ip,bi,bj).EQ.-1.) THEN
108     c kz=global2local_k(kpart(ip,bi,bj),bi,bj,mjtyhid)
109    
110     C recip_drF is in units 1/r (so IF r is in m this is in 1/m)
111     scalez=rkSign*recip_drF(kc)
112     C We should not do any special conversions for kz, since flt_trilinear
113     C expects it to be just a normal kpart type variable.
114     kz=kpart(ip,bi,bj)
115     CALL FLT_TRILINEAR(ix,jy,kz,u1,uVel,1,bi,bj,myThid)
116     CALL FLT_TRILINEAR(ix,jy,kz,v1,vVel,2,bi,bj,myThid)
117     CALL FLT_TRILINEAR(ix,jy,kz,w1,wVel,4,bi,bj,myThid)
118     ELSE
119     #else /* ALLOW_3D_FLT */
120     IF ( .TRUE. ) THEN
121     #endif /* ALLOW_3D_FLT */
122     CALL FLT_BILINEAR(ix,jy,u1,uVel,kc,1,bi,bj,myThid)
123     CALL FLT_BILINEAR(ix,jy,v1,vVel,kc,2,bi,bj,myThid)
124     ENDIF
125    
126     C When using this alternative scheme the noise probably should not be added twice.
127     #ifndef USE_FLT_ALT_NOISE
128 jmc 1.3 IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
129 jahn 1.1 u1 = u1 + u1*(PORT_RAND(seed)-0.5)*flt_noise
130     v1 = v1 + v1*(PORT_RAND(seed)-0.5)*flt_noise
131     #ifdef ALLOW_3D_FLT
132     #ifdef ALLOW_FLT_3D_NOISE
133     IF (iup(ip,bi,bj).EQ.-1.) THEN
134     w1 = w1 + w1*(PORT_RAND(seed)-0.5)*flt_noise
135     ENDIF
136     #endif
137     #endif /* ALLOW_3D_FLT */
138     ENDIF
139 jmc 1.3 #endif /* ndef USE_FLT_ALT_NOISE */
140 jahn 1.1
141     C ix and itx are in indices. Therefore it is necessary to multiply
142     C with a grid scale factor.
143    
144     itx1=ix+0.5*flt_deltaT*u1*scalex
145     jty1=jy+0.5*flt_deltaT*v1*scaley
146    
147     C Second step
148    
149     #ifdef ALLOW_3D_FLT
150     IF (iup(ip,bi,bj).EQ.-1.) THEN
151     ktz1=kz+0.5*flt_deltaT*w1*scalez
152     CALL FLT_TRILINEAR(itx1,jty1,ktz1,u2,uVel,
153     & 1,bi,bj,myThid)
154     CALL FLT_TRILINEAR(itx1,jty1,ktz1,v2,vVel,
155     & 2,bi,bj,myThid)
156     CALL FLT_TRILINEAR(itx1,jty1,ktz1,w2,wVel,
157     & 4,bi,bj,myThid)
158     ELSE
159     #else /* ALLOW_3D_FLT */
160     IF ( .TRUE. ) THEN
161     #endif /* ALLOW_3D_FLT */
162     CALL FLT_BILINEAR(itx1,jty1,u2,uVel,
163     & kc,1,bi,bj,myThid)
164     CALL FLT_BILINEAR(itx1,jty1,v2,vVel,
165     & kc,2,bi,bj,myThid)
166     ENDIF
167    
168 jmc 1.3 IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
169 jahn 1.1 #ifdef USE_FLT_ALT_NOISE
170 jmc 1.3 u2 = u2 + PORT_RAND_NORM()*flt_noise
171     v2 = v2 + PORT_RAND_NORM()*flt_noise
172 jahn 1.1 #ifdef ALLOW_3D_FLT
173     #ifdef ALLOW_FLT_3D_NOISE
174     IF (iup(ip,bi,bj).EQ.-1.) THEN
175 jmc 1.3 w2 = w2 + PORT_RAND_NORM()*flt_noise
176 jahn 1.1 ENDIF
177     #endif
178     #endif /* ALLOW_3D_FLT */
179    
180     #else /* USE_FLT_ALT_NOISE */
181     u2 = u2 + u2*(PORT_RAND(seed)-0.5)*flt_noise
182     v2 = v2 + v2*(PORT_RAND(seed)-0.5)*flt_noise
183     #ifdef ALLOW_3D_FLT
184     #ifdef ALLOW_FLT_3D_NOISE
185     IF (iup(ip,bi,bj).EQ.-1.) THEN
186     w2 = w2 + w2*(PORT_RAND(seed)-0.5)*flt_noise
187     ENDIF
188     #endif
189     #endif /* ALLOW_3D_FLT */
190    
191     #endif /* USE_FLT_ALT_NOISE */
192     ENDIF
193    
194     itx2=ix+0.5*flt_deltaT*u2*scalex
195     jty2=jy+0.5*flt_deltaT*v2*scaley
196    
197     C Third step
198    
199     #ifdef ALLOW_3D_FLT
200     IF (iup(ip,bi,bj).EQ.-1.) THEN
201     ktz2=kz+0.5*flt_deltaT*w2*scalez
202     CALL FLT_TRILINEAR(itx2,jty2,ktz2,u3,uVel,
203     & 1,bi,bj,myThid)
204     CALL FLT_TRILINEAR(itx2,jty2,ktz2,v3,vVel,
205     & 2,bi,bj,myThid)
206     CALL FLT_TRILINEAR(itx2,jty2,ktz2,w3,wVel,
207     & 4,bi,bj,myThid)
208     ELSE
209     #else /* ALLOW_3D_FLT */
210     IF ( .TRUE. ) THEN
211     #endif /* ALLOW_3D_FLT */
212     CALL FLT_BILINEAR(itx2,jty2,u3,uVel,
213     & kc,1,bi,bj,myThid)
214     CALL FLT_BILINEAR(itx2,jty2,v3,vVel,
215     & kc,2,bi,bj,myThid)
216     ENDIF
217    
218 jmc 1.3 IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
219 jahn 1.1 #ifdef USE_FLT_ALT_NOISE
220 jmc 1.3 u3 = u3 + PORT_RAND_NORM()*flt_noise
221     v3 = v3 + PORT_RAND_NORM()*flt_noise
222 jahn 1.1 #ifdef ALLOW_3D_FLT
223     #ifdef ALLOW_FLT_3D_NOISE
224     IF (iup(ip,bi,bj).EQ.-1.) THEN
225 jmc 1.3 w3 = w3 + PORT_RAND_NORM()*flt_noise
226 jahn 1.1 ENDIF
227     #endif
228     #endif /* ALLOW_3D_FLT */
229    
230     #else /* USE_FLT_ALT_NOISE */
231     u3 = u3 + u3*(PORT_RAND(seed)-0.5)*flt_noise
232     v3 = v3 + v3*(PORT_RAND(seed)-0.5)*flt_noise
233     #ifdef ALLOW_3D_FLT
234     #ifdef ALLOW_FLT_3D_NOISE
235     IF (iup(ip,bi,bj).EQ.-1.) THEN
236     w3 = w3 + w3*(PORT_RAND(seed)-0.5)*flt_noise
237     ENDIF
238     #endif
239     #endif /* ALLOW_3D_FLT */
240    
241     #endif /* USE_FLT_ALT_NOISE */
242     ENDIF
243    
244     itx3=ix+flt_deltaT*u3*scalex
245     jty3=jy+flt_deltaT*v3*scaley
246    
247     C Fourth step
248    
249     #ifdef ALLOW_3D_FLT
250     IF (iup(ip,bi,bj).EQ.-1.) THEN
251     ktz3=kz+0.5*flt_deltaT*w2*scalez
252     CALL FLT_TRILINEAR(itx3,jty3,ktz3,u4,uVel,
253     & 1,bi,bj,myThid)
254     CALL FLT_TRILINEAR(itx3,jty3,ktz3,v4,vVel,
255     & 2,bi,bj,myThid)
256     CALL FLT_TRILINEAR(itx3,jty3,ktz3,w4,wVel,
257     & 4,bi,bj,myThid)
258     ELSE
259     #else /* ALLOW_3D_FLT */
260     IF ( .TRUE. ) THEN
261     #endif /* ALLOW_3D_FLT */
262     CALL FLT_BILINEAR(itx3,jty3,u4,uVel,
263     & kc,1,bi,bj,myThid)
264     CALL FLT_BILINEAR(itx3,jty3,v4,vVel,
265     & kc,2,bi,bj,myThid)
266     ENDIF
267    
268 jmc 1.3 IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
269 jahn 1.1 #ifdef USE_FLT_ALT_NOISE
270 jmc 1.3 u4 = u4 + PORT_RAND_NORM()*flt_noise
271     v4 = v4 + PORT_RAND_NORM()*flt_noise
272 jahn 1.1 #ifdef ALLOW_3D_FLT
273     #ifdef ALLOW_FLT_3D_NOISE
274     IF (iup(ip,bi,bj).EQ.-1.) THEN
275 jmc 1.3 w4 = w4 + PORT_RAND_NORM()*flt_noise
276 jahn 1.1 ENDIF
277     #endif
278     #endif /* ALLOW_3D_FLT */
279    
280     #else /* USE_FLT_ALT_NOISE */
281     u4 = u4 + u4*(PORT_RAND(seed)-0.5)*flt_noise
282     v4 = v4 + v4*(PORT_RAND(seed)-0.5)*flt_noise
283     #ifdef ALLOW_3D_FLT
284     #ifdef ALLOW_FLT_3D_NOISE
285     IF (iup(ip,bi,bj).EQ.-1.) THEN
286     w4 = w4 + w4*(PORT_RAND(seed)-0.5)*flt_noise
287     ENDIF
288     #endif
289     #endif /* ALLOW_3D_FLT */
290    
291     #endif /* USE_FLT_ALT_NOISE */
292     ENDIF
293 jmc 1.3
294 jahn 1.1 C ipart is in coordinates. Therefore it is necessary to multiply
295     C with a grid scale factor divided by the number grid points per
296     C geographical coordinate.
297 jmc 1.3 ipart(ip,bi,bj) = ipart(ip,bi,bj) + flt_deltaT/6
298 jahn 1.1 & * (u1 + 2*u2 + 2*u3 + u4) * scalex
299 jmc 1.3 jpart(ip,bi,bj) = jpart(ip,bi,bj) + flt_deltaT/6
300 jahn 1.1 & * (v1 + 2*v2 + 2*v3 + v4) * scaley
301     #ifdef ALLOW_3D_FLT
302     IF (iup(ip,bi,bj).EQ.-1.) THEN
303     kpart(ip,bi,bj) = kpart(ip,bi,bj) + flt_deltaT/6
304     & * (w1 + 2*w2 + 2*w3 + w4) * scalez
305     ENDIF
306     #endif /* ALLOW_3D_FLT */
307    
308     C-- new horizontal grid indices
309 jmc 1.4 ic = MAX( 1-OLx, MIN( NINT(ipart(ip,bi,bj)), sNx+OLx ) )
310     jc = MAX( 1-OLy, MIN( NINT(jpart(ip,bi,bj)), sNy+OLy ) )
311 jahn 1.1
312     #ifdef ALLOW_3D_FLT
313     C If float is 3D, make sure that it remains in water
314     IF (iup(ip,bi,bj).EQ.-1.) THEN
315     C reflect on surface
316     IF (kpart(ip,bi,bj).LT.kzlo) kpart(ip,bi,bj)=kzlo
317     & +kzlo-kpart(ip,bi,bj)
318     C stop at bottom
319     IF (kpart(ip,bi,bj).GT.kzhi) kpart(ip,bi,bj)=kzhi
320     C-to work also with non flat-bottom set-up:
321     c IF ( kpart(ip,bi,bj).GT.kLowC(ic,jc,bi,bj)+0.5 )
322     c & kpart(ip,bi,bj) = kLowC(ic,jc,bi,bj)+0.5
323     ENDIF
324     #endif /* ALLOW_3D_FLT */
325    
326     #ifdef ALLOW_OBCS
327     IF ( useOBCS ) THEN
328     C-- stop floats which enter the OB region:
329     IF ( maskInC(ic,jc,bi,bj).EQ.0. .AND.
330     & maskC(ic,jc,1,bi,bj).EQ.1. ) THEN
331     C for now, just reset "tend" to myTime-deltaT
332     C (a better way would be to remove this one & re-order the list of floats)
333     tend(ip,bi,bj) = myTime - flt_deltaT
334     ENDIF
335     ENDIF
336     #endif /* ALLOW_OBCS */
337    
338     ENDIF
339     ENDIF
340    
341     C- end ip loop
342     ENDDO
343     C- end bi,bj loops
344     ENDDO
345     ENDDO
346    
347     RETURN
348     END

  ViewVC Help
Powered by ViewVC 1.1.22