Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90
- Property svn:eol-style deleted
r1152 r2528 11 11 !! flo_interp : interpolation 12 12 !!---------------------------------------------------------------------- 13 !! * Modules used14 13 USE flo_oce ! ocean drifting floats 15 14 USE oce ! ocean dynamics and tracers … … 20 19 PRIVATE 21 20 22 !! * Accessibility23 PUBLIC flo_4rk ! routine called by floats.F90 24 25 !! * Module variables26 REAL(wp), DIMENSION (4) :: & ! RK4 and Lagrange interpolation27 tcoef1 = (/ 1.0 , 0.5 , 0.5 , 0.0 /) , & ! coeffients for28 tcoef2 = (/ 0.0 , 0.5 , 0.5 , 1.0 /) , & ! lagrangian interp.29 scoef2 = (/ 1.0 , 2.0 , 2.0 , 1.0 /) , & ! RK4 coefficients30 rcoef = (/-1./6. , 1./2. ,-1./2. , 1./6. /) ! ??? 31 REAL(wp), DIMENSION (3) :: &32 scoef1 = (/ .5, .5, 1. /) ! compute position with interpolated 33 !!---------------------------------------------------------------------- 34 !! OPA 9.0 , LOCEAN-IPSL (2005)21 PUBLIC flo_4rk ! routine called by floats.F90 22 23 ! ! RK4 and Lagrange interpolation coefficients 24 REAL(wp), DIMENSION (4) :: tcoef1 = (/ 1.0 , 0.5 , 0.5 , 0.0 /) ! 25 REAL(wp), DIMENSION (4) :: tcoef2 = (/ 0.0 , 0.5 , 0.5 , 1.0 /) ! 26 REAL(wp), DIMENSION (4) :: scoef2 = (/ 1.0 , 2.0 , 2.0 , 1.0 /) ! 27 REAL(wp), DIMENSION (4) :: rcoef = (/-1./6. , 1./2. ,-1./2. , 1./6. /) ! 28 REAL(wp), DIMENSION (3) :: scoef1 = (/ 0.5 , 0.5 , 1.0 /) ! 29 30 !! * Substitutions 31 # include "domzgr_substitute.h90" 32 !!---------------------------------------------------------------------- 33 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 35 34 !! $Id$ 36 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 37 !!---------------------------------------------------------------------- 38 35 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 36 !!---------------------------------------------------------------------- 39 37 CONTAINS 40 38 … … 50 48 !! We need to know the velocity field, the old positions of the 51 49 !! floats and the grid defined on the domain. 50 !!---------------------------------------------------------------------- 51 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 52 !! 53 !!----------------------------------------------------------------------54 !! * Arguments55 INTEGER, INTENT(in) :: kt ! ocean time-step index56 57 !! * Local declarations58 53 INTEGER :: jfl, jind ! dummy loop indices 59 REAL(wp), DIMENSION ( jpnfl) :: & 60 zgifl, zgjfl, zgkfl, & ! index RK positions 61 zufl, zvfl, zwfl ! interpolated velocity at the 62 ! float position 63 REAL(wp), DIMENSION ( jpnfl, 4 ) :: & 64 zrkxfl, zrkyfl, zrkzfl ! RK coefficients 54 REAL(wp), DIMENSION(jpnfl) :: zgifl , zgjfl , zgkfl ! index RK positions 55 REAL(wp), DIMENSION(jpnfl) :: zufl , zvfl , zwfl ! interpolated velocity at the float position 56 REAL(wp), DIMENSION(jpnfl,4) :: zrkxfl, zrkyfl, zrkzfl ! RK coefficients 65 57 !!--------------------------------------------------------------------- 66 58 … … 127 119 END DO 128 120 129 DO jind = 1, 4 130 121 DO jind = 1, 4 122 131 123 ! for each step we compute the compute the velocity with Lagrange interpolation 132 CALL flo_interp( zgifl,zgjfl,zgkfl,zufl,zvfl,zwfl,jind)124 CALL flo_interp( zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind ) 133 125 134 126 ! computation of Runge-Kutta factor 135 136 127 DO jfl = 1, jpnfl 137 128 zrkxfl(jfl,jind) = rdt*zufl(jfl) … … 154 145 END DO 155 146 END DO 156 147 ! 157 148 END SUBROUTINE flo_4rk 158 149 159 150 160 151 SUBROUTINE flo_interp( pxt , pyt , pzt , & 161 & pufl, pvfl, pwfl, ki nd)152 & pufl, pvfl, pwfl, ki ) 162 153 !!---------------------------------------------------------------------- 163 154 !! *** ROUTINE flointerp *** … … 169 160 !! compute velocity at the date and the position we need to 170 161 !! integrated with RK method. 162 !!---------------------------------------------------------------------- 163 REAL(wp) , DIMENSION(jpnfl), INTENT(in ) :: pxt , pyt , pzt ! position of the float 164 REAL(wp) , DIMENSION(jpnfl), INTENT( out) :: pufl, pvfl, pwfl ! velocity at this position 165 INTEGER , INTENT(in ) :: ki ! 171 166 !! 172 !!---------------------------------------------------------------------- 173 !! * Local declarations 174 INTEGER :: & 175 kind, & 176 jfl, jind1, jind2, jind3, & 177 zsumu, zsumv, zsumw 178 INTEGER , DIMENSION ( jpnfl ) :: & 179 iilu, ijlu, iklu, & ! nearest neighbour INDEX-u 180 iilv, ijlv, iklv, & ! nearest neighbour INDEX-v 181 iilw, ijlw, iklw ! nearest neighbour INDEX-w 182 INTEGER , DIMENSION ( jpnfl, 4 ) :: & 183 iidu, ijdu, ikdu, & ! 64 nearest neighbour INDEX-u 184 iidv, ijdv, ikdv, & ! 64 nearest neighbour INDEX-v 185 iidw, ijdw, ikdw ! 64 nearest neighbour INDEX-w 186 REAL(wp) , DIMENSION ( jpnfl ) :: & 187 pxt , pyt , pzt, & ! position of the float 188 pufl, pvfl, pwfl ! velocity at this position 189 REAL(wp) , DIMENSION ( jpnfl, 4, 4, 4 ) :: & 190 ztufl, ztvfl, ztwfl ! velocity at choosen time step 191 REAL(wp) , DIMENSION ( jpnfl, 4 ) :: & 192 zlagxu, zlagyu, zlagzu, & ! Lagrange coefficients 193 zlagxv, zlagyv, zlagzv, & 194 zlagxw, zlagyw, zlagzw 167 INTEGER :: jfl, jind1, jind2, jind3 ! dummy loop indices 168 REAL(wp) :: zsumu, zsumv, zsumw ! local scalar 169 INTEGER , DIMENSION(jpnfl) :: iilu, ijlu, iklu ! nearest neighbour INDEX-u 170 INTEGER , DIMENSION(jpnfl) :: iilv, ijlv, iklv ! nearest neighbour INDEX-v 171 INTEGER , DIMENSION(jpnfl) :: iilw, ijlw, iklw ! nearest neighbour INDEX-w 172 INTEGER , DIMENSION(jpnfl,4) :: iidu, ijdu, ikdu ! 64 nearest neighbour INDEX-u 173 INTEGER , DIMENSION(jpnfl,4) :: iidv, ijdv, ikdv ! 64 nearest neighbour INDEX-v 174 INTEGER , DIMENSION(jpnfl,4) :: iidw, ijdw, ikdw ! 64 nearest neighbour INDEX-w 175 REAL(wp) , DIMENSION(jpnfl,4,4,4) :: ztufl , ztvfl , ztwfl ! velocity at choosen time step 176 REAL(wp) , DIMENSION(jpnfl,4) :: zlagxu, zlagyu, zlagzu ! Lagrange coefficients 177 REAL(wp) , DIMENSION(jpnfl,4) :: zlagxv, zlagyv, zlagzv ! - - 178 REAL(wp) , DIMENSION(jpnfl,4) :: zlagxw, zlagyw, zlagzw ! - - 195 179 !!--------------------------------------------------------------------- 196 180 … … 208 192 DO jfl = 1, jpnfl 209 193 ! i-direction 210 IF( iilu(jfl) <= 2 ) THEN 211 iidu(jfl,jind1) = jind1 212 ELSE 213 IF( iilu(jfl) >= jpi-1 ) THEN 214 iidu(jfl,jind1) = jpi + jind1 - 4 215 ELSE 216 iidu(jfl,jind1) = iilu(jfl) + jind1 - 2 194 IF( iilu(jfl) <= 2 ) THEN ; iidu(jfl,jind1) = jind1 195 ELSE 196 IF( iilu(jfl) >= jpi-1 ) THEN ; iidu(jfl,jind1) = jpi + jind1 - 4 197 ELSE ; iidu(jfl,jind1) = iilu(jfl) + jind1 - 2 217 198 ENDIF 218 199 ENDIF 219 200 ! j-direction 220 IF( ijlu(jfl) <= 2 ) THEN 221 ijdu(jfl,jind1) = jind1 222 ELSE 223 IF( ijlu(jfl) >= jpj-1 ) THEN 224 ijdu(jfl,jind1) = jpj + jind1 - 4 225 ELSE 226 ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2 201 IF( ijlu(jfl) <= 2 ) THEN ; ijdu(jfl,jind1) = jind1 202 ELSE 203 IF( ijlu(jfl) >= jpj-1 ) THEN ; ijdu(jfl,jind1) = jpj + jind1 - 4 204 ELSE ; ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2 227 205 ENDIF 228 206 ENDIF 229 207 ! k-direction 230 IF( iklu(jfl) <= 2 ) THEN 231 ikdu(jfl,jind1) = jind1 232 ELSE 233 IF( iklu(jfl) >= jpk-1 ) THEN 234 ikdu(jfl,jind1) = jpk + jind1 - 4 235 ELSE 236 ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2 208 IF( iklu(jfl) <= 2 ) THEN ; ikdu(jfl,jind1) = jind1 209 ELSE 210 IF( iklu(jfl) >= jpk-1 ) THEN ; ikdu(jfl,jind1) = jpk + jind1 - 4 211 ELSE ; ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2 237 212 ENDIF 238 213 ENDIF … … 241 216 242 217 ! Lagrange coefficients 243 244 218 DO jfl = 1, jpnfl 245 219 DO jind1 = 1, 4 … … 249 223 END DO 250 224 END DO 251 252 225 DO jind1 = 1, 4 253 226 DO jind2 = 1, 4 … … 269 242 DO jind3 = 1, 4 270 243 ztufl(jfl,jind1,jind2,jind3) = & 271 & ( tcoef1(ki nd) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) + &272 & tcoef2(ki nd) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) ) &244 & ( tcoef1(ki) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) + & 245 & tcoef2(ki) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) ) & 273 246 & / e1u(iidu(jfl,jind1),ijdu(jfl,jind2)) 274 247 END DO … … 301 274 DO jfl = 1, jpnfl 302 275 ! i-direction 303 IF( iilv(jfl) <= 2 ) THEN 304 iidv(jfl,jind1) = jind1 305 ELSE 306 IF( iilv(jfl) >= jpi-1 ) THEN 307 iidv(jfl,jind1) = jpi + jind1 - 4 308 ELSE 309 iidv(jfl,jind1) = iilv(jfl) + jind1 - 2 276 IF( iilv(jfl) <= 2 ) THEN ; iidv(jfl,jind1) = jind1 277 ELSE 278 IF( iilv(jfl) >= jpi-1 ) THEN ; iidv(jfl,jind1) = jpi + jind1 - 4 279 ELSE ; iidv(jfl,jind1) = iilv(jfl) + jind1 - 2 310 280 ENDIF 311 281 ENDIF 312 282 ! j-direction 313 IF( ijlv(jfl) <= 2 ) THEN 314 ijdv(jfl,jind1) = jind1 315 ELSE 316 IF( ijlv(jfl) >= jpj-1 ) THEN 317 ijdv(jfl,jind1) = jpj + jind1 - 4 318 ELSE 319 ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2 283 IF( ijlv(jfl) <= 2 ) THEN ; ijdv(jfl,jind1) = jind1 284 ELSE 285 IF( ijlv(jfl) >= jpj-1 ) THEN ; ijdv(jfl,jind1) = jpj + jind1 - 4 286 ELSE ; ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2 320 287 ENDIF 321 288 ENDIF 322 289 ! k-direction 323 IF( iklv(jfl) <= 2 ) THEN 324 ikdv(jfl,jind1) = jind1 325 ELSE 326 IF( iklv(jfl) >= jpk-1 ) THEN 327 ikdv(jfl,jind1) = jpk + jind1 - 4 328 ELSE 329 ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2 290 IF( iklv(jfl) <= 2 ) THEN ; ikdv(jfl,jind1) = jind1 291 ELSE 292 IF( iklv(jfl) >= jpk-1 ) THEN ; ikdv(jfl,jind1) = jpk + jind1 - 4 293 ELSE ; ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2 330 294 ENDIF 331 295 ENDIF … … 347 311 DO jfl = 1, jpnfl 348 312 IF( jind1 /= jind2 ) THEN 349 zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2)) ) )313 zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2)) ) ) 350 314 zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) ) 351 zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2)) ) )315 zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2)) ) ) 352 316 ENDIF 353 317 END DO … … 362 326 DO jind3 = 1 ,4 363 327 ztvfl(jfl,jind1,jind2,jind3)= & 364 & ( tcoef1(ki nd) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) + &365 & tcoef2(ki nd) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) ) &328 & ( tcoef1(ki) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) + & 329 & tcoef2(ki) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) ) & 366 330 & / e2v(iidv(jfl,jind1),ijdv(jfl,jind2)) 367 331 END DO … … 385 349 ! nearest neightboring point for computation of w 386 350 DO jfl = 1, jpnfl 387 iilw(jfl) = INT( pxt(jfl))388 ijlw(jfl) = INT( pyt(jfl))389 iklw(jfl) = INT( pzt(jfl)+.5)351 iilw(jfl) = INT( pxt(jfl) ) 352 ijlw(jfl) = INT( pyt(jfl) ) 353 iklw(jfl) = INT( pzt(jfl)+.5) 390 354 END DO 391 355 … … 394 358 DO jfl = 1, jpnfl 395 359 ! i-direction 396 IF( iilw(jfl) <= 2 ) THEN 397 iidw(jfl,jind1) = jind1 398 ELSE 399 IF( iilw(jfl) >= jpi-1 ) THEN 400 iidw(jfl,jind1) = jpi + jind1 - 4 401 ELSE 402 iidw(jfl,jind1) = iilw(jfl) + jind1 - 2 360 IF( iilw(jfl) <= 2 ) THEN ; iidw(jfl,jind1) = jind1 361 ELSE 362 IF( iilw(jfl) >= jpi-1 ) THEN ; iidw(jfl,jind1) = jpi + jind1 - 4 363 ELSE ; iidw(jfl,jind1) = iilw(jfl) + jind1 - 2 403 364 ENDIF 404 365 ENDIF 405 366 ! j-direction 406 IF( ijlw(jfl) <= 2 ) THEN 407 ijdw(jfl,jind1) = jind1 408 ELSE 409 IF( ijlw(jfl) >= jpj-1 ) THEN 410 ijdw(jfl,jind1) = jpj + jind1 - 4 411 ELSE 412 ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2 367 IF( ijlw(jfl) <= 2 ) THEN ; ijdw(jfl,jind1) = jind1 368 ELSE 369 IF( ijlw(jfl) >= jpj-1 ) THEN ; ijdw(jfl,jind1) = jpj + jind1 - 4 370 ELSE ; ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2 413 371 ENDIF 414 372 ENDIF 415 373 ! k-direction 416 IF( iklw(jfl) <= 2 ) THEN 417 ikdw(jfl,jind1) = jind1 418 ELSE 419 IF( iklw(jfl) >= jpk-1 ) THEN 420 ikdw(jfl,jind1) = jpk + jind1 - 4 421 ELSE 422 ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2 423 ENDIF 424 ENDIF 425 END DO 426 END DO 427 DO jind1 = 1, 4 428 DO jfl = 1, jpnfl 429 IF( iklw(jfl) <= 2 ) THEN 430 ikdw(jfl,jind1) = jind1 431 ELSE 432 IF( iklw(jfl) >= jpk-1 ) THEN 433 ikdw(jfl,jind1) = jpk + jind1 - 4 434 ELSE 435 ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2 374 IF( iklw(jfl) <= 2 ) THEN ; ikdw(jfl,jind1) = jind1 375 ELSE 376 IF( iklw(jfl) >= jpk-1 ) THEN ; ikdw(jfl,jind1) = jpk + jind1 - 4 377 ELSE ; ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2 378 ENDIF 379 ENDIF 380 END DO 381 END DO 382 DO jind1 = 1, 4 383 DO jfl = 1, jpnfl 384 IF( iklw(jfl) <= 2 ) THEN ; ikdw(jfl,jind1) = jind1 385 ELSE 386 IF( iklw(jfl) >= jpk-1 ) THEN ; ikdw(jfl,jind1) = jpk + jind1 - 4 387 ELSE ; ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2 436 388 ENDIF 437 389 ENDIF … … 440 392 441 393 ! Lagrange coefficients for w interpolation 442 443 394 DO jfl = 1, jpnfl 444 395 DO jind1 = 1, 4 … … 448 399 END DO 449 400 END DO 450 451 401 DO jind1 = 1, 4 452 402 DO jind2 = 1, 4 453 403 DO jfl = 1, jpnfl 454 404 IF( jind1 /= jind2 ) THEN 455 zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2)) ) )456 zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2)) ) )405 zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2)) ) ) 406 zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2)) ) ) 457 407 zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) ) 458 408 ENDIF … … 462 412 463 413 ! velocity w when we compute at middle time step 464 465 414 DO jfl = 1, jpnfl 466 415 DO jind1 = 1, 4 … … 468 417 DO jind3 = 1, 4 469 418 ztwfl(jfl,jind1,jind2,jind3)= & 470 & ( tcoef1(kind) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+ & 471 & tcoef2(kind) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) ) & 472 !!bug e3w instead of fse3 473 & / e3w(ikdw(jfl,jind3)) 419 & ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+ & 420 & tcoef2(ki) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) ) & 421 & / fse3w(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) 474 422 END DO 475 423 END DO 476 424 END DO 477 425 478 zsumw =0.426 zsumw = 0.e0 479 427 DO jind1 = 1, 4 480 428 DO jind2 = 1, 4 … … 487 435 pwfl(jfl) = zsumw 488 436 END DO 489 437 ! 490 438 END SUBROUTINE flo_interp 491 439 492 440 # else 493 441 !!---------------------------------------------------------------------- 494 !! Default option Empty module 495 !!---------------------------------------------------------------------- 496 CONTAINS 497 SUBROUTINE flo_4rk ! Empty routine 498 END SUBROUTINE flo_4rk 442 !! No floats Dummy module 443 !!---------------------------------------------------------------------- 499 444 #endif 500 445
Note: See TracChangeset
for help on using the changeset viewer.