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
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 IMPLICIT NONE
0020
0021
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
0031 _RL myTime
0032 INTEGER myIter, myThid
0033
0034
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
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
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
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
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
0102
0103 #ifdef ALLOW_3D_FLT
0104 IF (iup(ip,bi,bj).EQ.-1.) THEN
0105
0106
0107
0108 scalez=rkSign*recip_drF(kc)
0109
0110
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
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
0139
0140
0141 itx1=ix+0.5*flt_deltaT*u1*scalex
0142 jty1=jy+0.5*flt_deltaT*v1*scaley
0143
0144
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
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
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
0292
0293
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
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
0311 IF (iup(ip,bi,bj).EQ.-1.) THEN
0312
0313 IF (kpart(ip,bi,bj).LT.kzlo) kpart(ip,bi,bj)=kzlo
0314 & +kzlo-kpart(ip,bi,bj)
0315
0316 IF (kpart(ip,bi,bj).GT.kzhi) kpart(ip,bi,bj)=kzhi
0317
0318
0319
0320 ENDIF
0321 #endif /* ALLOW_3D_FLT */
0322
0323 #ifdef ALLOW_OBCS
0324 IF ( useOBCS ) THEN
0325
0326 IF ( maskInC(ic,jc,bi,bj).EQ.0. .AND.
0327 & maskC(ic,jc,1,bi,bj).EQ.1. ) THEN
0328
0329
0330 tend(ip,bi,bj) = myTime - flt_deltaT
0331 ENDIF
0332 ENDIF
0333 #endif /* ALLOW_OBCS */
0334
0335 ENDIF
0336 ENDIF
0337
0338
0339 ENDDO
0340
0341 ENDDO
0342 ENDDO
0343
0344 RETURN
0345 END