Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:40:50 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
3394409633 Oliv*0001 #include "FLT_OPTIONS.h"
c21f2e71d0 Jean*0002 #undef _USE_INTEGERS
3394409633 Oliv*0003 
                0004       SUBROUTINE FLT_RUNGA4 (
                0005      I                        myTime, myIter, myThid )
                0006 
                0007 C     ==================================================================
                0008 C     SUBROUTINE FLT_RUNGA4
                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 
                0018 C     !USES:
                0019       IMPLICIT NONE
                0020 
                0021 C     == global variables ==
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
                0026 #include "DYNVARS.h"
730d8469b1 Oliv*0027 #include "FLT_SIZE.h"
3394409633 Oliv*0028 #include "FLT.h"
                0029 
                0030 C     == routine arguments ==
                0031       _RL myTime
                0032       INTEGER myIter, myThid
                0033 
                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 */
3394409633 Oliv*0047 
                0048 C     == local variables ==
                0049       INTEGER bi, bj
                0050       INTEGER ip
                0051       INTEGER ic, jc, kc, iG, jG
                0052       _RL u1, v1, u2, v2, u3, v3, u4, v4
                0053 #ifdef ALLOW_3D_FLT
                0054       _RL w1, w2, w3, w4, ktz1, ktz2, ktz3, kz, scalez
                0055       _RL kzlo, kzhi
                0056 #endif
                0057       _RL ix, jy, itx1, jty1, itx2, jty2, itx3, jty3
                0058       _RL scalex, scaley
c21f2e71d0 Jean*0059 
                0060 C     == end of interface ==
                0061 
                0062 #ifndef USE_FLT_ALT_NOISE
3394409633 Oliv*0063 #ifdef _USE_INTEGERS
                0064       seed = -1
                0065 #else
                0066       seed = -1.d0
                0067 #endif
c21f2e71d0 Jean*0068 #endif /* ndef USE_FLT_ALT_NOISE */
3394409633 Oliv*0069 
                0070 #ifdef ALLOW_3D_FLT
                0071       kzlo = 0.5 _d 0
                0072       kzhi = 0.5 _d 0 + DFLOAT(Nr)
                0073 #endif
                0074 
                0075       DO bj=myByLo(myThid),myByHi(myThid)
                0076        DO bi=myBxLo(myThid),myBxHi(myThid)
                0077          DO ip=1,npart_tile(bi,bj)
                0078 
                0079 C     If float has died move to level 0
                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
                0084 C     Start integration between tstart and tend (individual for each float)
                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
                0089 
                0090               ix = ipart(ip,bi,bj)
                0091               jy = jpart(ip,bi,bj)
                0092               ic=NINT(ix)
                0093               jc=NINT(jy)
                0094               kc=NINT(kpart(ip,bi,bj))
                0095 
                0096               scalex=recip_dxF(ic,jc,bi,bj)
                0097               scaley=recip_dyF(ic,jc,bi,bj)
c21f2e71d0 Jean*0098               iG = myXGlobalLo + (bi-1)*sNx + ic-1
3394409633 Oliv*0099               jG = myYGlobalLo + (bj-1)*sNy + jc-1
                0100 
                0101 C     First step
                0102 
                0103 #ifdef ALLOW_3D_FLT
                0104               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0105 c               kz=global2local_k(kpart(ip,bi,bj),bi,bj,mjtyhid)
                0106 
                0107 C recip_drF is in units 1/r (so IF r is in m this is in 1/m)
                0108                 scalez=rkSign*recip_drF(kc)
                0109 C We should not do any special conversions for kz, since flt_trilinear
                0110 C expects it to be just a normal kpart type variable.
                0111                 kz=kpart(ip,bi,bj)
                0112                 CALL FLT_TRILINEAR(ix,jy,kz,u1,uVel,1,bi,bj,myThid)
                0113                 CALL FLT_TRILINEAR(ix,jy,kz,v1,vVel,2,bi,bj,myThid)
                0114                 CALL FLT_TRILINEAR(ix,jy,kz,w1,wVel,4,bi,bj,myThid)
                0115               ELSE
                0116 #else  /* ALLOW_3D_FLT */
                0117               IF ( .TRUE. ) THEN
                0118 #endif /* ALLOW_3D_FLT */
                0119                 CALL FLT_BILINEAR(ix,jy,u1,uVel,kc,1,bi,bj,myThid)
                0120                 CALL FLT_BILINEAR(ix,jy,v1,vVel,kc,2,bi,bj,myThid)
                0121               ENDIF
                0122 
                0123 C When using this alternative scheme the noise probably should not be added twice.
                0124 #ifndef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0125               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
3394409633 Oliv*0126                 u1 = u1 + u1*(PORT_RAND(seed)-0.5)*flt_noise
                0127                 v1 = v1 + v1*(PORT_RAND(seed)-0.5)*flt_noise
                0128 #ifdef ALLOW_3D_FLT
                0129 #ifdef ALLOW_FLT_3D_NOISE
                0130                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0131                   w1 = w1 + w1*(PORT_RAND(seed)-0.5)*flt_noise
                0132                 ENDIF
                0133 #endif
                0134 #endif /* ALLOW_3D_FLT */
                0135               ENDIF
c21f2e71d0 Jean*0136 #endif /* ndef USE_FLT_ALT_NOISE */
3394409633 Oliv*0137 
                0138 C ix and itx are in indices. Therefore it is necessary to multiply
                0139 C with a grid scale factor.
                0140 
                0141               itx1=ix+0.5*flt_deltaT*u1*scalex
                0142               jty1=jy+0.5*flt_deltaT*v1*scaley
                0143 
                0144 C     Second step
                0145 
                0146 #ifdef ALLOW_3D_FLT
                0147               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0148                 ktz1=kz+0.5*flt_deltaT*w1*scalez
                0149                 CALL FLT_TRILINEAR(itx1,jty1,ktz1,u2,uVel,
                0150      &               1,bi,bj,myThid)
                0151                 CALL FLT_TRILINEAR(itx1,jty1,ktz1,v2,vVel,
                0152      &               2,bi,bj,myThid)
                0153                 CALL FLT_TRILINEAR(itx1,jty1,ktz1,w2,wVel,
                0154      &               4,bi,bj,myThid)
                0155               ELSE
                0156 #else  /* ALLOW_3D_FLT */
                0157               IF ( .TRUE. ) THEN
                0158 #endif /* ALLOW_3D_FLT */
                0159                 CALL FLT_BILINEAR(itx1,jty1,u2,uVel,
                0160      &               kc,1,bi,bj,myThid)
                0161                 CALL FLT_BILINEAR(itx1,jty1,v2,vVel,
                0162      &               kc,2,bi,bj,myThid)
                0163               ENDIF
                0164 
c21f2e71d0 Jean*0165               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
3394409633 Oliv*0166 #ifdef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0167                 u2 = u2 + PORT_RAND_NORM()*flt_noise
                0168                 v2 = v2 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0169 #ifdef ALLOW_3D_FLT
                0170 #ifdef ALLOW_FLT_3D_NOISE
                0171                 IF (iup(ip,bi,bj).EQ.-1.) THEN
c21f2e71d0 Jean*0172                   w2 = w2 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0173                 ENDIF
                0174 #endif
                0175 #endif /* ALLOW_3D_FLT */
                0176 
                0177 #else /* USE_FLT_ALT_NOISE */
                0178                 u2 = u2 + u2*(PORT_RAND(seed)-0.5)*flt_noise
                0179                 v2 = v2 + v2*(PORT_RAND(seed)-0.5)*flt_noise
                0180 #ifdef ALLOW_3D_FLT
                0181 #ifdef ALLOW_FLT_3D_NOISE
                0182                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0183                   w2 = w2 + w2*(PORT_RAND(seed)-0.5)*flt_noise
                0184                 ENDIF
                0185 #endif
                0186 #endif /* ALLOW_3D_FLT */
                0187 
                0188 #endif /* USE_FLT_ALT_NOISE */
                0189               ENDIF
                0190 
                0191               itx2=ix+0.5*flt_deltaT*u2*scalex
                0192               jty2=jy+0.5*flt_deltaT*v2*scaley
                0193 
                0194 C     Third step
                0195 
                0196 #ifdef ALLOW_3D_FLT
                0197               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0198                 ktz2=kz+0.5*flt_deltaT*w2*scalez
                0199                 CALL FLT_TRILINEAR(itx2,jty2,ktz2,u3,uVel,
                0200      &               1,bi,bj,myThid)
                0201                 CALL FLT_TRILINEAR(itx2,jty2,ktz2,v3,vVel,
                0202      &               2,bi,bj,myThid)
                0203                 CALL FLT_TRILINEAR(itx2,jty2,ktz2,w3,wVel,
                0204      &               4,bi,bj,myThid)
                0205               ELSE
                0206 #else  /* ALLOW_3D_FLT */
                0207               IF ( .TRUE. ) THEN
                0208 #endif /* ALLOW_3D_FLT */
                0209                 CALL FLT_BILINEAR(itx2,jty2,u3,uVel,
                0210      &               kc,1,bi,bj,myThid)
                0211                 CALL FLT_BILINEAR(itx2,jty2,v3,vVel,
                0212      &               kc,2,bi,bj,myThid)
                0213               ENDIF
                0214 
c21f2e71d0 Jean*0215               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
3394409633 Oliv*0216 #ifdef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0217                 u3 = u3 + PORT_RAND_NORM()*flt_noise
                0218                 v3 = v3 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0219 #ifdef ALLOW_3D_FLT
                0220 #ifdef ALLOW_FLT_3D_NOISE
                0221                 IF (iup(ip,bi,bj).EQ.-1.) THEN
c21f2e71d0 Jean*0222                   w3 = w3 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0223                 ENDIF
                0224 #endif
                0225 #endif /* ALLOW_3D_FLT */
                0226 
                0227 #else /* USE_FLT_ALT_NOISE */
                0228                 u3 = u3 + u3*(PORT_RAND(seed)-0.5)*flt_noise
                0229                 v3 = v3 + v3*(PORT_RAND(seed)-0.5)*flt_noise
                0230 #ifdef ALLOW_3D_FLT
                0231 #ifdef ALLOW_FLT_3D_NOISE
                0232                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0233                   w3 = w3 + w3*(PORT_RAND(seed)-0.5)*flt_noise
                0234                 ENDIF
                0235 #endif
                0236 #endif /* ALLOW_3D_FLT */
                0237 
                0238 #endif /* USE_FLT_ALT_NOISE */
                0239               ENDIF
                0240 
                0241               itx3=ix+flt_deltaT*u3*scalex
                0242               jty3=jy+flt_deltaT*v3*scaley
                0243 
                0244 C     Fourth step
                0245 
                0246 #ifdef ALLOW_3D_FLT
                0247               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0248                 ktz3=kz+0.5*flt_deltaT*w2*scalez
                0249                 CALL FLT_TRILINEAR(itx3,jty3,ktz3,u4,uVel,
                0250      &               1,bi,bj,myThid)
                0251                 CALL FLT_TRILINEAR(itx3,jty3,ktz3,v4,vVel,
                0252      &               2,bi,bj,myThid)
                0253                 CALL FLT_TRILINEAR(itx3,jty3,ktz3,w4,wVel,
                0254      &               4,bi,bj,myThid)
                0255               ELSE
                0256 #else  /* ALLOW_3D_FLT */
                0257               IF ( .TRUE. ) THEN
                0258 #endif /* ALLOW_3D_FLT */
                0259                 CALL FLT_BILINEAR(itx3,jty3,u4,uVel,
                0260      &               kc,1,bi,bj,myThid)
                0261                 CALL FLT_BILINEAR(itx3,jty3,v4,vVel,
                0262      &               kc,2,bi,bj,myThid)
                0263               ENDIF
                0264 
c21f2e71d0 Jean*0265               IF ( flt_noise.NE.0. .AND. iup(ip,bi,bj).NE.-2. ) THEN
3394409633 Oliv*0266 #ifdef USE_FLT_ALT_NOISE
c21f2e71d0 Jean*0267                 u4 = u4 + PORT_RAND_NORM()*flt_noise
                0268                 v4 = v4 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0269 #ifdef ALLOW_3D_FLT
                0270 #ifdef ALLOW_FLT_3D_NOISE
                0271                 IF (iup(ip,bi,bj).EQ.-1.) THEN
c21f2e71d0 Jean*0272                   w4 = w4 + PORT_RAND_NORM()*flt_noise
3394409633 Oliv*0273                 ENDIF
                0274 #endif
                0275 #endif /* ALLOW_3D_FLT */
                0276 
                0277 #else /* USE_FLT_ALT_NOISE */
                0278                 u4 = u4 + u4*(PORT_RAND(seed)-0.5)*flt_noise
                0279                 v4 = v4 + v4*(PORT_RAND(seed)-0.5)*flt_noise
                0280 #ifdef ALLOW_3D_FLT
                0281 #ifdef ALLOW_FLT_3D_NOISE
                0282                 IF (iup(ip,bi,bj).EQ.-1.) THEN
                0283                   w4 = w4 + w4*(PORT_RAND(seed)-0.5)*flt_noise
                0284                 ENDIF
                0285 #endif
                0286 #endif /* ALLOW_3D_FLT */
                0287 
                0288 #endif /* USE_FLT_ALT_NOISE */
                0289               ENDIF
c21f2e71d0 Jean*0290 
3394409633 Oliv*0291 C ipart is in coordinates. Therefore it is necessary to multiply
                0292 C with a grid scale factor divided by the number grid points per
                0293 C geographical coordinate.
c21f2e71d0 Jean*0294               ipart(ip,bi,bj) = ipart(ip,bi,bj) + flt_deltaT/6
3394409633 Oliv*0295      &                        * (u1 + 2*u2 + 2*u3 + u4) * scalex
c21f2e71d0 Jean*0296               jpart(ip,bi,bj) = jpart(ip,bi,bj) + flt_deltaT/6
3394409633 Oliv*0297      &                        * (v1 + 2*v2 + 2*v3 + v4) * scaley
                0298 #ifdef ALLOW_3D_FLT
                0299               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0300                 kpart(ip,bi,bj) = kpart(ip,bi,bj) + flt_deltaT/6
                0301      &                          * (w1 + 2*w2 + 2*w3 + w4) * scalez
                0302               ENDIF
                0303 #endif /* ALLOW_3D_FLT */
                0304 
                0305 C--  new horizontal grid indices
e1fb02e8f0 Jean*0306               ic = MAX( 1-OLx, MIN( NINT(ipart(ip,bi,bj)), sNx+OLx ) )
                0307               jc = MAX( 1-OLy, MIN( NINT(jpart(ip,bi,bj)), sNy+OLy ) )
3394409633 Oliv*0308 
                0309 #ifdef ALLOW_3D_FLT
                0310 C If float is 3D, make sure that it remains in water
                0311               IF (iup(ip,bi,bj).EQ.-1.) THEN
                0312 C reflect on surface
                0313                 IF (kpart(ip,bi,bj).LT.kzlo) kpart(ip,bi,bj)=kzlo
                0314      &                                      +kzlo-kpart(ip,bi,bj)
                0315 C stop at bottom
                0316                 IF (kpart(ip,bi,bj).GT.kzhi) kpart(ip,bi,bj)=kzhi
                0317 C-to work also with non flat-bottom set-up:
                0318 c               IF ( kpart(ip,bi,bj).GT.kLowC(ic,jc,bi,bj)+0.5 )
                0319 c    &               kpart(ip,bi,bj) = kLowC(ic,jc,bi,bj)+0.5
                0320               ENDIF
                0321 #endif /* ALLOW_3D_FLT */
                0322 
                0323 #ifdef ALLOW_OBCS
                0324               IF ( useOBCS ) THEN
                0325 C--  stop floats which enter the OB region:
                0326                IF ( maskInC(ic,jc,bi,bj).EQ.0. .AND.
                0327      &              maskC(ic,jc,1,bi,bj).EQ.1. ) THEN
                0328 C    for now, just reset "tend" to myTime-deltaT
                0329 C    (a better way would be to remove this one & re-order the list of floats)
                0330                   tend(ip,bi,bj) = myTime - flt_deltaT
                0331                ENDIF
                0332               ENDIF
                0333 #endif /* ALLOW_OBCS */
                0334 
                0335             ENDIF
                0336            ENDIF
                0337 
                0338 C-    end ip loop
                0339          ENDDO
                0340 C-    end bi,bj loops
                0341        ENDDO
                0342       ENDDO
                0343 
                0344       RETURN
                0345       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!