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

Diff of /MITgcm/pkg/flt/flt_runga2.F

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

revision 1.2 by molod, Mon Jul 19 15:13:07 2004 UTC revision 1.3 by edhill, Tue Sep 7 16:19:30 2004 UTC
# Line 13  c     ================================== Line 13  c     ==================================
13  c     SUBROUTINE flt_runga2  c     SUBROUTINE flt_runga2
14  c     ==================================================================  c     ==================================================================
15  c  c
16  c     o This routine steps floats forward with second order Runga-Kutta  c     o This routine steps floats forward with second order Runge-Kutta
17    c
18    c     started: Arne Biastoch
19    c
20    c     changed: 2004.06.10 Antti Westerlund (antti.westerlund@helsinki.fi)
21    c              and Sergio Jaramillo (sju@eos.ubc.ca)
22  c  c
23  c     ==================================================================  c     ==================================================================
24  c     SUBROUTINE flt_runga2  c     SUBROUTINE flt_runga2
# Line 38  c     == routine arguments == Line 43  c     == routine arguments ==
43        INTEGER bi, bj        INTEGER bi, bj
44        _RL global2local_i        _RL global2local_i
45        _RL global2local_j        _RL global2local_j
46          _RL global2local_k
47  c     == local variables ==  c     == local variables ==
48    
49        integer ip, kp, iG, jG        integer ip, kp, iG, jG
# Line 50  c     == local variables == Line 55  c     == local variables ==
55        _RL scalex, scaley        _RL scalex, scaley
56        character*(max_len_mbuf) msgbuf        character*(max_len_mbuf) msgbuf
57        _RL npart_dist        _RL npart_dist
58        Real*8 PORT_RAND  #ifdef USE_FLT_ALT_NOISE
59  #undef _USE_INTEGERS        Real*8 PORT_RAND_NORM
 #ifdef _USE_INTEGERS  
       integer seed  
60  #else  #else
61        Real*8 seed        Real*8 PORT_RAND
62  #endif  #endif
63    
64  c     == end of interface ==  c     == end of interface ==
65    
66        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
67        DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
68                do ip=1,npart_tile(bi,bj)
       do ip=1,npart_tile(bi,bj)  
69    
70  c     If float has died move to level 0  c     If float has died move to level 0
71  c  c
72           if(                 if(
73       &  (tend(ip,bi,bj).ne.-1. .and. myCurrentTime.gt. tend(ip,bi,bj))       & (tend(ip,bi,bj).ne.-1. .and. myCurrentTime.gt. tend(ip,bi,bj)))
74       &      ) then       & then
75                      kpart(ip,bi,bj) = 0.
76              kpart(ip,bi,bj) = 0.                 else
   
          else  
77  c     Start integration between tstart and tend (individual for each float)  c     Start integration between tstart and tend (individual for each float)
78  c  c
79           if(                    if(
80       &  (tstart(ip,bi,bj).eq.-1. .or. myCurrentTime.ge.tstart(ip,bi,bj))       & (tstart(ip,bi,bj).eq.-1. .or. myCurrentTime.ge.tstart(ip,bi,bj))
81       &   .and.       &  .and.
82       &  (  tend(ip,bi,bj).eq.-1. .or. myCurrentTime.le.  tend(ip,bi,bj))       & (  tend(ip,bi,bj).eq.-1. .or. myCurrentTime.le.  tend(ip,bi,bj))
83       &   .and.       & .and.
84       &  (   iup(ip,bi,bj).ne. -3.)       & (   iup(ip,bi,bj).ne. -3.)
85       &      ) then       & ) then
86    
87  c     Convert to local indices  c     Convert to local indices
88  c  c
             xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)  
             yy=global2local_j(ypart(ip,bi,bj),bi,bj,mythid)  
             kp=INT(kpart(ip,bi,bj))  
89    
90              scalex=recip_dxF(INT(xx),INT(yy),bi,bj)  C Note: global2local_i and global2local_j use delX and delY.
91              scaley=recip_dyF(INT(xx),INT(yy),bi,bj)  C This may be a problem, especially if you are using a curvilinear
92              iG = myXGlobalLo + (bi-1)*sNx  C grid. More information below.
93              jG = myYGlobalLo + (bj-1)*sNy                       xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)
94                         yy=global2local_j(ypart(ip,bi,bj),bi,bj,mythid)
95                         kp=INT(kpart(ip,bi,bj))
96    
97                         scalex=recip_dxF(INT(xx),INT(yy),bi,bj)
98                         scaley=recip_dyF(INT(xx),INT(yy),bi,bj)
99                         iG = myXGlobalLo + (bi-1)*sNx
100                         jG = myYGlobalLo + (bj-1)*sNy
101    
102    
103  #ifdef ALLOW_3D_FLT  #ifdef ALLOW_3D_FLT
104              if (iup(ip,bi,bj).eq.-1.) then                       if (iup(ip,bi,bj).eq.-1.) then
105                 scalez=recip_drF(kp)  c                        zz=global2local_k(kpart(ip,bi,bj),bi,bj,mythid)
106                 zt=global2local_j(kpart(ip,bi,bj),bi,bj,mythid)  
107                 call flt_bilinear3D(xx,yy,uu,zp,uVel,2,bi,bj)  c recip_drF is in units 1/r (so if r is in m this is in 1/m)
108                 call flt_bilinear3D(xx,yy,vv,zp,vVel,3,bi,bj)                          scalez=recip_drF(kp)
109                 call flt_bilinear3D(zz,yy,ww,zp,wVel,4,bi,bj)  c We shouldn't do any special conversions for zz, since flt_trilinear
110                 zt=zz+0.5*deltaTmom*zz*scalez  c expects it to be just a normal kpart type variable.
111              else                          zz=kpart(ip,bi,bj)
112                            call flt_trilinear(xx,yy,zz,uu,uVel,2,bi,bj)
113                            call flt_trilinear(xx,yy,zz,vv,vVel,3,bi,bj)
114                            call flt_trilinear(zz,yy,zz,ww,wVel,4,bi,bj)
115                            zt=zz+0.5*deltaTmom*ww*scalez
116                         else
117  #endif  #endif
118                 call flt_bilinear(xx,yy,uu,kp,uVel,2,bi,bj)                          call flt_bilinear(xx,yy,uu,kp,uVel,2,bi,bj)
119                 call flt_bilinear(xx,yy,vv,kp,vVel,3,bi,bj)                          call flt_bilinear(xx,yy,vv,kp,vVel,3,bi,bj)
120  #ifdef ALLOW_3D_FLT  #ifdef ALLOW_3D_FLT
121              endif                       endif
122  #endif  #endif
123    
124              if (iup(ip,bi,bj).ne.-2.) then  #ifdef USE_FLT_ALT_NOISE
125                 uu = uu + uu*(PORT_RAND(seed)-0.5)*flt_noise  c When using this alternative scheme the noise probably should not be added twice.
126                 vv = vv + vv*(PORT_RAND(seed)-0.5)*flt_noise  #else
127              endif                       if (iup(ip,bi,bj).ne.-2.) then
128                            uu = uu + uu*(PORT_RAND()-0.5)*flt_noise
129                            vv = vv + vv*(PORT_RAND()-0.5)*flt_noise
130    #ifdef ALLOW_3D_FLT
131    #ifdef ALLOW_FLT_3D_NOISE
132                            if (iup(ip,bi,bj).eq.-1.) then
133                               ww = ww + ww*(PORT_RAND()-0.5)*flt_noise
134                            endif
135    #endif
136    #endif
137                         endif
138    #endif
139    
140  c xx and xt are in indices. Therefore it is necessary to multiply  c xx and xt are in indices. Therefore it is necessary to multiply
141  c with a grid scale factor.  c with a grid scale factor.
142  c  c
143              xt=xx+0.5*deltaTmom*uu*scalex                       xt=xx+0.5*deltaTmom*uu*scalex
144              yt=yy+0.5*deltaTmom*vv*scaley                       yt=yy+0.5*deltaTmom*vv*scaley
145    
146  c     Second step  c     Second step
147  c  c
148                            
149  #ifdef ALLOW_3D_FLT  #ifdef ALLOW_3D_FLT
150              if (iup(ip,bi,bj).eq.-1.) then                       if (iup(ip,bi,bj).eq.-1.) then
151                 call flt_bilinear3D(xt,yt,u1,zt,uVel,2,bi,bj)                          call flt_trilinear(xt,yt,zt,u1,uVel,2,bi,bj)
152                 call flt_bilinear3D(xt,yt,v1,zt,vVel,3,bi,bj)                          call flt_trilinear(xt,yt,zt,v1,vVel,3,bi,bj)
153                 call flt_bilinear3D(xx,yy,w1,zt,wVel,4,bi,bj)                          call flt_trilinear(xt,yt,zt,w1,wVel,4,bi,bj)
154                 kpart(ip,bi,bj) = kpart(ip,bi,bj) + deltaTmom*w1*scalez                       else
155              else  #endif
156  #endif                          call flt_bilinear(xt,yt,u1,kp,uVel,2,bi,bj)
157              call flt_bilinear(xt,yt,u1,kp,uVel,2,bi,bj)                          call flt_bilinear(xt,yt,v1,kp,vVel,3,bi,bj)
158              call flt_bilinear(xt,yt,v1,kp,vVel,3,bi,bj)  #ifdef ALLOW_3D_FLT
159  #ifdef ALLOW_3D_FLT                       endif
160              endif  #endif
161  #endif  
162                         if (iup(ip,bi,bj).ne.-2.) then
163              if (iup(ip,bi,bj).ne.-2.) then  #ifdef USE_FLT_ALT_NOISE
164                 u1 = u1 + u1*(PORT_RAND(seed)-0.5)*flt_noise                          u1 = u1 + port_rand_norm()*flt_noise
165                 v1 = v1 + v1*(PORT_RAND(seed)-0.5)*flt_noise                          v1 = v1 + port_rand_norm()*flt_noise
166              endif  #ifdef ALLOW_3D_FLT
167    #ifdef ALLOW_FLT_3D_NOISE
168                            if (iup(ip,bi,bj).eq.-1.) then
169                               w1 = w1 + port_rand_norm()*flt_noise
170                            endif
171    #endif
172    #endif
173    
174    #else
175                            u1 = u1 + u1*(PORT_RAND()-0.5)*flt_noise
176                            v1 = v1 + v1*(PORT_RAND()-0.5)*flt_noise
177    #ifdef ALLOW_3D_FLT
178    #ifdef ALLOW_FLT_3D_NOISE
179                            if (iup(ip,bi,bj).eq.-1.) then
180                               w1 = w1 + w1*(PORT_RAND()-0.5)*flt_noise
181                            endif
182    #endif
183    #endif
184    
185    #endif
186                         endif
187    
188  c xpart is in coordinates. Therefore it is necessary to multiply  c xpart is in coordinates. Therefore it is necessary to multiply
189  c with a grid scale factor divided by the number grid points per  c with a grid scale factor divided by the number grid points per
190  c geographical coordinate.  c geographical coordinate.
191  c  c
192              xpart(ip,bi,bj) = xpart(ip,bi,bj)  C This will only work if delX & delY are available.
193       &                        + deltaTmom*u1*scalex*delX(iG)  C This may be a problem, especially if you are using a curvilinear
194              ypart(ip,bi,bj) = ypart(ip,bi,bj)  C grid. In that case you have to replace them for the values of
195       &                        + deltaTmom*v1*scaley*delY(jG)  C your grid, which can be troublesome.
196                         xpart(ip,bi,bj) = xpart(ip,bi,bj)
197           endif       &                    + deltaTmom*u1*scalex*delX(iG)
198        endif                       ypart(ip,bi,bj) = ypart(ip,bi,bj)
199         &                    + deltaTmom*v1*scaley*delY(jG)
200        enddo  #ifdef ALLOW_3D_FLT
201                         if (iup(ip,bi,bj).eq.-1.) then
202                            kpart(ip,bi,bj) = kpart(ip,bi,bj)
203         &                                  + deltaTmom*w1*scalez
204                         endif
205    #endif
206    
207        ENDDO  #ifdef ALLOW_3D_FLT
208    c If float is 3D, make sure that it remains in water
209                         if (iup(ip,bi,bj).eq.-1.) then
210    c reflect on surface
211                            if(kpart(ip,bi,bj).lt.1.0)
212         &                       kpart(ip,bi,bj)=1.0
213         &                       +abs(1.0-kpart(ip,bi,bj))
214    c stop at bottom
215                            if(kpart(ip,bi,bj).gt.Nr)
216         &                       kpart(ip,bi,bj)=Nr
217                         endif
218    #endif
219                      endif
220                   endif
221                
222                enddo
223             ENDDO
224        ENDDO        ENDDO
225  c  c
226        return        return
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22