pkg/flt/flt_runga2.F

0a3ae49bfc Jean*0001 #include "FLT_OPTIONS.h"
c21f2e71d0 Jean*0002 #undef _USE_INTEGERS
c806179eb4 Alis*0003
eacecc7041 Jean*0004       SUBROUTINE FLT_RUNGA2 (
                0005      I                        myTime, myIter, myThid )
                0006
                0007 C     ==================================================================
51ec3c32fe Jean*0008 C     SUBROUTINE FLT_RUNGA2
eacecc7041 Jean*0009 C     ==================================================================
                0010 C     o This routine steps floats forward with second order Runge-Kutta
                0011 C
                0012 C     started: Arne Biastoch
                0013 C
                0014 C     changed: 2004.06.10 Antti Westerlund (antti.westerlund@helsinki.fi)
                0015 C              and Sergio Jaramillo (sju@eos.ubc.ca)
                0016 C     ==================================================================
                0017
51ec3c32fe Jean*0018 C     !USES:
                0019       IMPLICIT NONE
c806179eb4 Alis*0020
51ec3c32fe Jean*0021 C     == global variables ==
c806179eb4 Alis*0022 #include "SIZE.h"
51ec3c32fe Jean*0023 #include "EEPARAMS.h"
c806179eb4 Alis*0024 #include "PARAMS.h"
                0025 #include "GRID.h"
51ec3c32fe Jean*0026 #include "DYNVARS.h"
730d8469b1 Oliv*0027 #include "FLT_SIZE.h"
c806179eb4 Alis*0028 #include "FLT.h"
                0029
eacecc7041 Jean*0030 C     == routine arguments ==
                0031       _RL myTime
                0032       INTEGER myIter, myThid
                0033
51ec3c32fe Jean*0034 C     == Functions ==
c21f2e71d0 Jean*0035 #ifdef USE_FLT_ALT_NOISE
                0036       Real*8   PORT_RAND_NORM
                0037       EXTERNAL PORT_RAND_NORM
                0038 #else
                0039       Real*8   PORT_RAND
                0040       EXTERNAL PORT_RAND
                0041 #ifdef _USE_INTEGERS
                0042       INTEGER seed
                0043 #else
                0044       Real*8 seed
                0045 #endif
                0046 #endif /* USE_FLT_ALT_NOISE */
c806179eb4 Alis*0047
51ec3c32fe Jean*0048 C     == local variables ==
                0049       INTEGER bi, bj
                0050       INTEGER ip
                0051       INTEGER ic, jc, kc, iG, jG
eacecc7041 Jean*0052       _RL uu, vv, u1, v1
c806179eb4 Alis*0053 #ifdef ALLOW_3D_FLT
d5477ff298 Jean*0054       _RL ww, w1, ktz, kz, scalez
51ec3c32fe Jean*0055       _RL kzlo, kzhi
c806179eb4 Alis*0056 #endif
d5477ff298 Jean*0057       _RL ix, jy, itx, jty
c806179eb4 Alis*0058       _RL scalex, scaley
c21f2e71d0 Jean*0059
                0060 C     == end of interface ==
                0061
                0062 #ifndef USE_FLT_ALT_NOISE
e3c94c7d93 Ed H*0063 #ifdef _USE_INTEGERS
c25f76287c Jean*0064       seed = -1
e3c94c7d93 Ed H*0065 #else
c25f76287c Jean*0066       seed = -1.d0
e3c94c7d93 Ed H*0067 #endif
c21f2e71d0 Jean*0068 #endif /* ndef USE_FLT_ALT_NOISE */
c806179eb4 Alis*0069
51ec3c32fe Jean*0070 #ifdef ALLOW_3D_FLT
                0071       kzlo = 0.5 _d 0
                0072       kzhi = 0.5 _d 0 + DFLOAT(Nr)
                0073 #endif
                0074
c806179eb4 Alis*0075       DO bj=myByLo(myThid),myByHi(myThid)
51ec3c32fe Jean*0076        DO bi=myBxLo(myThid),myBxHi(myThid)
                0077          DO ip=1,npart_tile(bi,bj)
eacecc7041 Jean*0078
                0079 C     If float has died move to level 0
51ec3c32fe Jean*0080            IF ( tend(ip,bi,bj).NE.-1. .AND. myTime.GT.tend(ip,bi,bj)
                0081      &        ) THEN
                0082             kpart(ip,bi,bj) = 0.
                0083            ELSE
eacecc7041 Jean*0084 C     Start integration between tstart and tend (individual for each float)
51ec3c32fe Jean*0085             IF ( (tstart(ip,bi,bj).EQ.-1..OR.myTime.GE.tstart(ip,bi,bj))
                0086      &      .AND.(  tend(ip,bi,bj).EQ.-1..OR.myTime.LE.  tend(ip,bi,bj))
                0087      &      .AND.(   iup(ip,bi,bj).NE.-3.)
                0088      &         ) THEN
eacecc7041 Jean*0089
d5477ff298 Jean*0090               ix = ipart(ip,bi,bj)
                0091               jy = jpart(ip,bi,bj)
                0092               ic=NINT(ix)
                0093               jc=NINT(jy)
51ec3c32fe Jean*0094               kc=NINT(kpart(ip,bi,bj))
c806179eb4 Alis*0095
51ec3c32fe Jean*0096               scalex=recip_dxF(ic,jc,bi,bj)
                0097               scaley=recip_dyF(ic,jc,bi,bj)
                0098               iG = myXGlobalLo + (bi-1)*sNx + ic-1
                0099               jG = myYGlobalLo + (bj-1)*sNy + jc-1
c806179eb4 Alis*0100
                0101 #ifdef ALLOW_3D_FLT
51ec3c32fe Jean*0102               IF (iup(ip,bi,bj).EQ.-1.) THEN
d5477ff298 Jean*0103 c               kz=global2local_k(kpart(ip,bi,bj),bi,bj,mjtyhid)
e47631944b Ed H*0104
eacecc7041 Jean*0105 C recip_drF is in units 1/r (so IF r is in m this is in 1/m)
51ec3c32fe Jean*0106                 scalez=rkSign*recip_drF(kc)
d5477ff298 Jean*0107 C We should not do any special conversions for kz, since flt_trilinear
eacecc7041 Jean*0108 C expects it to be just a normal kpart type variable.
d5477ff298 Jean*0109                 kz=kpart(ip,bi,bj)
                0110                 CALL FLT_TRILINEAR(ix,jy,kz,uu,uVel,1,bi,bj,myThid)
                0111                 CALL FLT_TRILINEAR(ix,jy,kz,vv,vVel,2,bi,bj,myThid)
                0112                 CALL FLT_TRILINEAR(ix,jy,kz,ww,wVel,4,bi,bj,myThid)
51ec3c32fe Jean*0113               ELSE
                0114 #else  /* ALLOW_3D_FLT */
                0115               IF ( .TRUE. ) THEN
                0116 #endif /* ALLOW_3D_FLT */
d5477ff298 Jean*0117                 CALL FLT_BILINEAR(ix,jy,uu,uVel,kc,1,bi,bj,myThid)
                0118                 CALL FLT_BILINEAR(ix,jy,vv,vVel,kc,2,bi,bj,myThid)
51ec3c32fe Jean*0119               ENDIF
c806179eb4 Alis*0120
eacecc7041 Jean*0121 C When using this alternative scheme the noise probably should not be added twice.
51ec3c32fe Jean*0122 #ifndef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0123               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
51ec3c32fe Jean*0124                 uu = uu + uu*(PORT_RAND(seed)-0.5)*flt_noise
                0125                 vv = vv + vv*(PORT_RAND(seed)-0.5)*flt_noise
e47631944b Ed H*0126 #ifdef ALLOW_3D_FLT
                0127 #ifdef ALLOW_FLT_3D_NOISE
51ec3c32fe Jean*0128                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0129                   ww = ww + ww*(PORT_RAND(seed)-0.5)*flt_noise
                0130                 ENDIF
e47631944b Ed H*0131 #endif
51ec3c32fe Jean*0132 #endif /* ALLOW_3D_FLT */
                0133               ENDIF
c21f2e71d0 Jean*0134 #endif /* ndef USE_FLT_ALT_NOISE */
c806179eb4 Alis*0135
d5477ff298 Jean*0136 C ix and itx are in indices. Therefore it is necessary to multiply
eacecc7041 Jean*0137 C with a grid scale factor.
                0138
c84b484b2a Davi*0139               itx=ix+0.5*flt_deltaT*uu*scalex
                0140               jty=jy+0.5*flt_deltaT*vv*scaley
c806179eb4 Alis*0141
eacecc7041 Jean*0142 C     Second step
                0143
c806179eb4 Alis*0144 #ifdef ALLOW_3D_FLT
51ec3c32fe Jean*0145               IF (iup(ip,bi,bj).EQ.-1.) THEN
c84b484b2a Davi*0146                 ktz=kz+0.5*flt_deltaT*ww*scalez
d5477ff298 Jean*0147                 CALL FLT_TRILINEAR(itx,jty,ktz,u1,uVel,1,bi,bj,myThid)
                0148                 CALL FLT_TRILINEAR(itx,jty,ktz,v1,vVel,2,bi,bj,myThid)
                0149                 CALL FLT_TRILINEAR(itx,jty,ktz,w1,wVel,4,bi,bj,myThid)
51ec3c32fe Jean*0150               ELSE
                0151 #else  /* ALLOW_3D_FLT */
                0152               IF ( .TRUE. ) THEN
                0153 #endif /* ALLOW_3D_FLT */
d5477ff298 Jean*0154                 CALL FLT_BILINEAR(itx,jty,u1,uVel,kc,1,bi,bj,myThid)
                0155                 CALL FLT_BILINEAR(itx,jty,v1,vVel,kc,2,bi,bj,myThid)
51ec3c32fe Jean*0156               ENDIF
                0157
c21f2e71d0 Jean*0158               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
e47631944b Ed H*0159 #ifdef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0160                 u1 = u1 + PORT_RAND_NORM()*flt_noise
                0161                 v1 = v1 + PORT_RAND_NORM()*flt_noise
e47631944b Ed H*0162 #ifdef ALLOW_3D_FLT
                0163 #ifdef ALLOW_FLT_3D_NOISE
51ec3c32fe Jean*0164                 IF (iup(ip,bi,bj).EQ.-1.) THEN
c21f2e71d0 Jean*0165                   w1 = w1 + PORT_RAND_NORM()*flt_noise
51ec3c32fe Jean*0166                 ENDIF
e47631944b Ed H*0167 #endif
51ec3c32fe Jean*0168 #endif /* ALLOW_3D_FLT */
e47631944b Ed H*0169
51ec3c32fe Jean*0170 #else /* USE_FLT_ALT_NOISE */
                0171                 u1 = u1 + u1*(PORT_RAND(seed)-0.5)*flt_noise
                0172                 v1 = v1 + v1*(PORT_RAND(seed)-0.5)*flt_noise
e47631944b Ed H*0173 #ifdef ALLOW_3D_FLT
                0174 #ifdef ALLOW_FLT_3D_NOISE
51ec3c32fe Jean*0175                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0176                   w1 = w1 + w1*(PORT_RAND(seed)-0.5)*flt_noise
                0177                 ENDIF
e47631944b Ed H*0178 #endif
51ec3c32fe Jean*0179 #endif /* ALLOW_3D_FLT */
e47631944b Ed H*0180
51ec3c32fe Jean*0181 #endif /* USE_FLT_ALT_NOISE */
                0182               ENDIF
c806179eb4 Alis*0183
d5477ff298 Jean*0184 C ipart is in coordinates. Therefore it is necessary to multiply
eacecc7041 Jean*0185 C with a grid scale factor divided by the number grid points per
                0186 C geographical coordinate.
d5477ff298 Jean*0187               ipart(ip,bi,bj) = ipart(ip,bi,bj)
c84b484b2a Davi*0188      &                        + flt_deltaT*u1*scalex
d5477ff298 Jean*0189               jpart(ip,bi,bj) = jpart(ip,bi,bj)
c84b484b2a Davi*0190      &                        + flt_deltaT*v1*scaley
e47631944b Ed H*0191 #ifdef ALLOW_3D_FLT
51ec3c32fe Jean*0192               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0193                 kpart(ip,bi,bj) = kpart(ip,bi,bj)
c84b484b2a Davi*0194      &                          + flt_deltaT*w1*scalez
51ec3c32fe Jean*0195               ENDIF
                0196 #endif /* ALLOW_3D_FLT */
c806179eb4 Alis*0197
ac1cd9b608 Jean*0198 C--  new horizontal grid indices
e1fb02e8f0 Jean*0199               ic = MAX( 1-OLx, MIN( NINT(ipart(ip,bi,bj)), sNx+OLx ) )
                0200               jc = MAX( 1-OLy, MIN( NINT(jpart(ip,bi,bj)), sNy+OLy ) )
ac1cd9b608 Jean*0201
e47631944b Ed H*0202 #ifdef ALLOW_3D_FLT
eacecc7041 Jean*0203 C If float is 3D, make sure that it remains in water
51ec3c32fe Jean*0204               IF (iup(ip,bi,bj).EQ.-1.) THEN
eacecc7041 Jean*0205 C reflect on surface
51ec3c32fe Jean*0206                 IF (kpart(ip,bi,bj).LT.kzlo) kpart(ip,bi,bj)=kzlo
                0207      &                                      +kzlo-kpart(ip,bi,bj)
eacecc7041 Jean*0208 C stop at bottom
51ec3c32fe Jean*0209                 IF (kpart(ip,bi,bj).GT.kzhi) kpart(ip,bi,bj)=kzhi
ac1cd9b608 Jean*0210 C-to work also with non flat-bottom set-up:
                0211 c               IF ( kpart(ip,bi,bj).GT.kLowC(ic,jc,bi,bj)+0.5 )
                0212 c    &               kpart(ip,bi,bj) = kLowC(ic,jc,bi,bj)+0.5
51ec3c32fe Jean*0213               ENDIF
                0214 #endif /* ALLOW_3D_FLT */
ac1cd9b608 Jean*0215
                0216 #ifdef ALLOW_OBCS
                0217               IF ( useOBCS ) THEN
                0218 C--  stop floats which enter the OB region:
                0219                IF ( maskInC(ic,jc,bi,bj).EQ.0. .AND.
                0220      &              maskC(ic,jc,1,bi,bj).EQ.1. ) THEN
                0221 C    for now, just reset "tend" to myTime-deltaT
                0222 C    (a better way would be to remove this one & re-order the list of floats)
                0223                   tend(ip,bi,bj) = myTime - flt_deltaT
                0224                ENDIF
                0225               ENDIF
                0226 #endif /* ALLOW_OBCS */
                0227
51ec3c32fe Jean*0228             ENDIF
                0229            ENDIF
eacecc7041 Jean*0230
51ec3c32fe Jean*0231 C-    end ip loop
e47631944b Ed H*0232          ENDDO
51ec3c32fe Jean*0233 C-    end bi,bj loops
                0234        ENDDO
c806179eb4 Alis*0235       ENDDO
eacecc7041 Jean*0236
                0237       RETURN
                0238       END