/[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.1 by adcroft, Thu Sep 13 17:43:56 2001 UTC revision 1.5 by edhill, Mon Sep 13 16:13:42 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    #ifdef USE_FLT_ALT_NOISE
59          Real*8 PORT_RAND_NORM
60    #else
61        Real*8 PORT_RAND        Real*8 PORT_RAND
62    #undef _USE_INTEGERS
63    #ifdef _USE_INTEGERS
64          integer seed
65    #else
66          Real*8 seed
67    #endif
68    #endif
69    
70  c     == end of interface ==  c     == end of interface ==
71    
72        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
73        DO bi=myBxLo(myThid),myBxHi(myThid)           DO bi=myBxLo(myThid),myBxHi(myThid)
74                do ip=1,npart_tile(bi,bj)
       do ip=1,npart_tile(bi,bj)  
75    
76  c     If float has died move to level 0  c     If float has died move to level 0
77  c  c
78           if(                 if(
79       &  (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)))
80       &      ) then       & then
81                      kpart(ip,bi,bj) = 0.
82              kpart(ip,bi,bj) = 0.                 else
   
          else  
83  c     Start integration between tstart and tend (individual for each float)  c     Start integration between tstart and tend (individual for each float)
84  c  c
85           if(                    if(
86       &  (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))
87       &   .and.       &  .and.
88       &  (  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))
89       &   .and.       & .and.
90       &  (   iup(ip,bi,bj).ne. -3.)       & (   iup(ip,bi,bj).ne. -3.)
91       &      ) then       & ) then
92    
93  c     Convert to local indices  c     Convert to local indices
94  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))  
95    
96              scalex=recip_dxF(INT(xx),INT(yy),bi,bj)  C Note: global2local_i and global2local_j use delX and delY.
97              scaley=recip_dyF(INT(xx),INT(yy),bi,bj)  C This may be a problem, especially if you are using a curvilinear
98              iG = myXGlobalLo + (bi-1)*sNx  C grid. More information below.
99              jG = myYGlobalLo + (bj-1)*sNy                       xx=global2local_i(xpart(ip,bi,bj),bi,bj,mythid)
100                         yy=global2local_j(ypart(ip,bi,bj),bi,bj,mythid)
101                         kp=INT(kpart(ip,bi,bj))
102    
103                         scalex=recip_dxF(INT(xx),INT(yy),bi,bj)
104                         scaley=recip_dyF(INT(xx),INT(yy),bi,bj)
105                         iG = myXGlobalLo + (bi-1)*sNx
106                         jG = myYGlobalLo + (bj-1)*sNy
107    
108    
109  #ifdef ALLOW_3D_FLT  #ifdef ALLOW_3D_FLT
110              if (iup(ip,bi,bj).eq.-1.) then                       if (iup(ip,bi,bj).eq.-1.) then
111                 scalez=recip_drF(kp)  c                        zz=global2local_k(kpart(ip,bi,bj),bi,bj,mythid)
112                 zt=global2local_j(kpart(ip,bi,bj),bi,bj,mythid)  
113                 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)
114                 call flt_bilinear3D(xx,yy,vv,zp,vVel,3,bi,bj)                          scalez=recip_drF(kp)
115                 call flt_bilinear3D(zz,yy,ww,zp,wVel,4,bi,bj)  c We should not do any special conversions for zz, since flt_trilinear
116                 zt=zz+0.5*deltaTmom*zz*scalez  c expects it to be just a normal kpart type variable.
117              else                          zz=kpart(ip,bi,bj)
118                            call flt_trilinear(xx,yy,zz,uu,uVel,2,bi,bj)
119                            call flt_trilinear(xx,yy,zz,vv,vVel,3,bi,bj)
120                            call flt_trilinear(zz,yy,zz,ww,wVel,4,bi,bj)
121                            zt=zz+0.5*deltaTmom*ww*scalez
122                         else
123    #endif
124                            call flt_bilinear(xx,yy,uu,kp,uVel,2,bi,bj)
125                            call flt_bilinear(xx,yy,vv,kp,vVel,3,bi,bj)
126    #ifdef ALLOW_3D_FLT
127                         endif
128    #endif
129    
130    #ifdef USE_FLT_ALT_NOISE
131    c When using this alternative scheme the noise probably should not be added twice.
132    #else
133                         if (iup(ip,bi,bj).ne.-2.) then
134                            uu = uu + uu*(PORT_RAND(seed)-0.5)*flt_noise
135                            vv = vv + vv*(PORT_RAND(seed)-0.5)*flt_noise
136    #ifdef ALLOW_3D_FLT
137    #ifdef ALLOW_FLT_3D_NOISE
138                            if (iup(ip,bi,bj).eq.-1.) then
139                               ww = ww + ww*(PORT_RAND(seed)-0.5)*flt_noise
140                            endif
141  #endif  #endif
                call flt_bilinear(xx,yy,uu,kp,uVel,2,bi,bj)  
                call flt_bilinear(xx,yy,vv,kp,vVel,3,bi,bj)  
 #ifdef ALLOW_3D_FLT  
             endif  
142  #endif  #endif
143                         endif
144              if (iup(ip,bi,bj).ne.-2.) then  #endif
                uu = uu + uu*(PORT_RAND()-0.5)*flt_noise  
                vv = vv + vv*(PORT_RAND()-0.5)*flt_noise  
             endif  
145    
146  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
147  c with a grid scale factor.  c with a grid scale factor.
148  c  c
149              xt=xx+0.5*deltaTmom*uu*scalex                       xt=xx+0.5*deltaTmom*uu*scalex
150              yt=yy+0.5*deltaTmom*vv*scaley                       yt=yy+0.5*deltaTmom*vv*scaley
151    
152  c     Second step  c     Second step
153  c  c
154                            
155  #ifdef ALLOW_3D_FLT  #ifdef ALLOW_3D_FLT
156              if (iup(ip,bi,bj).eq.-1.) then                       if (iup(ip,bi,bj).eq.-1.) then
157                 call flt_bilinear3D(xt,yt,u1,zt,uVel,2,bi,bj)                          call flt_trilinear(xt,yt,zt,u1,uVel,2,bi,bj)
158                 call flt_bilinear3D(xt,yt,v1,zt,vVel,3,bi,bj)                          call flt_trilinear(xt,yt,zt,v1,vVel,3,bi,bj)
159                 call flt_bilinear3D(xx,yy,w1,zt,wVel,4,bi,bj)                          call flt_trilinear(xt,yt,zt,w1,wVel,4,bi,bj)
160                 kpart(ip,bi,bj) = kpart(ip,bi,bj) + deltaTmom*w1*scalez                       else
161              else  #endif
162  #endif                          call flt_bilinear(xt,yt,u1,kp,uVel,2,bi,bj)
163              call flt_bilinear(xt,yt,u1,kp,uVel,2,bi,bj)                          call flt_bilinear(xt,yt,v1,kp,vVel,3,bi,bj)
164              call flt_bilinear(xt,yt,v1,kp,vVel,3,bi,bj)  #ifdef ALLOW_3D_FLT
165  #ifdef ALLOW_3D_FLT                       endif
166              endif  #endif
167  #endif  
168                         if (iup(ip,bi,bj).ne.-2.) then
169              if (iup(ip,bi,bj).ne.-2.) then  #ifdef USE_FLT_ALT_NOISE
170                 u1 = u1 + u1*(PORT_RAND()-0.5)*flt_noise                          u1 = u1 + port_rand_norm()*flt_noise
171                 v1 = v1 + v1*(PORT_RAND()-0.5)*flt_noise                          v1 = v1 + port_rand_norm()*flt_noise
172              endif  #ifdef ALLOW_3D_FLT
173    #ifdef ALLOW_FLT_3D_NOISE
174                            if (iup(ip,bi,bj).eq.-1.) then
175                               w1 = w1 + port_rand_norm()*flt_noise
176                            endif
177    #endif
178    #endif
179    
180    #else
181                            u1 = u1 + u1*(PORT_RAND(seed)-0.5)*flt_noise
182                            v1 = v1 + v1*(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                               w1 = w1 + w1*(PORT_RAND(seed)-0.5)*flt_noise
187                            endif
188    #endif
189    #endif
190    
191    #endif
192                         endif
193    
194  c xpart is in coordinates. Therefore it is necessary to multiply  c xpart is in coordinates. Therefore it is necessary to multiply
195  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
196  c geographical coordinate.  c geographical coordinate.
197  c  c
198              xpart(ip,bi,bj) = xpart(ip,bi,bj)  C This will only work if delX & delY are available.
199       &                        + deltaTmom*u1*scalex*delX(iG)  C This may be a problem, especially if you are using a curvilinear
200              ypart(ip,bi,bj) = ypart(ip,bi,bj)  C grid. In that case you have to replace them for the values of
201       &                        + deltaTmom*v1*scaley*delY(jG)  C your grid, which can be troublesome.
202                         xpart(ip,bi,bj) = xpart(ip,bi,bj)
203           endif       &                    + deltaTmom*u1*scalex*delX(iG)
204        endif                       ypart(ip,bi,bj) = ypart(ip,bi,bj)
205         &                    + deltaTmom*v1*scaley*delY(jG)
206        enddo  #ifdef ALLOW_3D_FLT
207                         if (iup(ip,bi,bj).eq.-1.) then
208                            kpart(ip,bi,bj) = kpart(ip,bi,bj)
209         &                                  + deltaTmom*w1*scalez
210                         endif
211    #endif
212    
213        ENDDO  #ifdef ALLOW_3D_FLT
214    c If float is 3D, make sure that it remains in water
215                         if (iup(ip,bi,bj).eq.-1.) then
216    c reflect on surface
217                            if(kpart(ip,bi,bj).lt.1.0)
218         &                       kpart(ip,bi,bj)=1.0
219         &                       +abs(1.0-kpart(ip,bi,bj))
220    c stop at bottom
221                            if(kpart(ip,bi,bj).gt.Nr)
222         &                       kpart(ip,bi,bj)=Nr
223                         endif
224    #endif
225                      endif
226                   endif
227                
228                enddo
229             ENDDO
230        ENDDO        ENDDO
231  c  c
232        return        return
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.5
  ViewVC Help
Powered by ViewVC 1.1.22