- Timestamp:
- 2014-12-01T11:36:36+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4879_UKMO_NOC_MERGE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r4666 r4934 37 37 PRIVATE 38 38 39 PUBLIC tra_adv_tvd ! routine called by step.F90 39 PUBLIC tra_adv_tvd ! routine called by traadv.F90 40 PUBLIC tra_adv_tvd_zts ! routine called by traadv.F90 40 41 41 42 LOGICAL :: l_trd ! flag to compute trends … … 262 263 END SUBROUTINE tra_adv_tvd 263 264 265 SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 266 & ptb, ptn, pta, kjpt ) 267 !!---------------------------------------------------------------------- 268 !! *** ROUTINE tra_adv_tvd_zts *** 269 !! 270 !! ** Purpose : Compute the now trend due to total advection of 271 !! tracers and add it to the general trend of tracer equations 272 !! 273 !! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with 274 !! corrected flux (monotonic correction). This version use sub- 275 !! timestepping for the vertical advection which increases stability 276 !! when vertical metrics are small. 277 !! note: - this advection scheme needs a leap-frog time scheme 278 !! 279 !! ** Action : - update (pta) with the now advective tracer trends 280 !! - save the trends 281 !!---------------------------------------------------------------------- 282 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 283 ! 284 INTEGER , INTENT(in ) :: kt ! ocean time-step index 285 INTEGER , INTENT(in ) :: kit000 ! first time step index 286 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 287 INTEGER , INTENT(in ) :: kjpt ! number of tracers 288 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 289 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 290 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 291 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 292 ! 293 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 294 REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep 295 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 296 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 297 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 298 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 299 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 300 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 301 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 302 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 303 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav 304 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 305 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 306 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 307 !!---------------------------------------------------------------------- 308 ! 309 IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd_zts') 310 ! 311 CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 312 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 313 CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 314 ! 315 IF( kt == kit000 ) THEN 316 IF(lwp) WRITE(numout,*) 317 IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 318 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 319 ENDIF 320 ! 321 l_trd = .FALSE. 322 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 323 ! 324 IF( l_trd ) THEN 325 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 326 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 327 ENDIF 328 ! 329 zwi(:,:,:) = 0._wp 330 z_rzts = 1._wp / REAL( jnzts, wp ) 331 zr_p2dt(:) = 1._wp / p2dt(:) 332 ! 333 ! ! =========== 334 DO jn = 1, kjpt ! tracer loop 335 ! ! =========== 336 ! 1. Bottom value : flux set to zero 337 ! ---------------------------------- 338 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 339 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 340 341 ! 2. upstream advection with initial mass fluxes & intermediate update 342 ! -------------------------------------------------------------------- 343 ! upstream tracer flux in the i and j direction 344 DO jk = 1, jpkm1 345 DO jj = 1, jpjm1 346 DO ji = 1, fs_jpim1 ! vector opt. 347 ! upstream scheme 348 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 349 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 350 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 351 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 352 zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 353 zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 354 END DO 355 END DO 356 END DO 357 358 ! upstream tracer flux in the k direction 359 ! Surface value 360 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0._wp ! volume variable 361 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 362 ENDIF 363 ! Interior value 364 DO jk = 2, jpkm1 365 DO jj = 1, jpj 366 DO ji = 1, jpi 367 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 368 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 369 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 370 END DO 371 END DO 372 END DO 373 374 ! total advective trend 375 DO jk = 1, jpkm1 376 z2dtt = p2dt(jk) 377 DO jj = 2, jpjm1 378 DO ji = fs_2, fs_jpim1 ! vector opt. 379 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 380 ! total intermediate advective trends 381 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 382 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 383 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 384 ! update and guess with monotonic sheme 385 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 386 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 387 END DO 388 END DO 389 END DO 390 ! ! Lateral boundary conditions on zwi (unchanged sign) 391 CALL lbc_lnk( zwi, 'T', 1. ) 392 393 ! ! trend diagnostics (contribution of upstream fluxes) 394 IF( l_trd ) THEN 395 ! store intermediate advective trends 396 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 397 END IF 398 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 399 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 400 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 401 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 402 ENDIF 403 404 ! 3. antidiffusive flux : high order minus low order 405 ! -------------------------------------------------- 406 ! antidiffusive flux on i and j 407 408 409 DO jk = 1, jpkm1 410 411 DO jj = 1, jpjm1 412 DO ji = 1, fs_jpim1 ! vector opt. 413 zwx_sav(ji,jj) = zwx(ji,jj,jk) 414 zwy_sav(ji,jj) = zwy(ji,jj,jk) 415 416 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 417 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 418 END DO 419 END DO 420 421 DO jj = 2, jpjm1 ! partial horizontal divergence 422 DO ji = fs_2, fs_jpim1 423 zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 424 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 425 END DO 426 END DO 427 428 DO jj = 1, jpjm1 429 DO ji = 1, fs_jpim1 ! vector opt. 430 zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 431 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 432 END DO 433 END DO 434 END DO 435 436 ! antidiffusive flux on k 437 zwz(:,:,1) = 0._wp ! Surface value 438 zwz_sav(:,:,:) = zwz(:,:,:) 439 ! 440 ztrs(:,:,:,1) = ptb(:,:,:,jn) 441 zwzts(:,:,:) = 0._wp 442 443 DO jl = 1, jnzts ! Start of sub timestepping loop 444 445 IF( jl == 1 ) THEN ! Euler forward to kick things off 446 jtb = 1 ; jtn = 1 ; jta = 2 447 zts(:) = p2dt(:) * z_rzts 448 jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux 449 ! starting at jl =1 if jnzts is odd; 450 ! starting at jl =2 otherwise 451 ELSEIF( jl == 2 ) THEN ! First leapfrog step 452 jtb = 1 ; jtn = 2 ; jta = 3 453 zts(:) = 2._wp * p2dt(:) * z_rzts 454 ELSE ! Shuffle pointers for subsequent leapfrog steps 455 jtb = MOD(jtb,3) + 1 456 jtn = MOD(jtn,3) + 1 457 jta = MOD(jta,3) + 1 458 ENDIF 459 DO jk = 2, jpkm1 ! Interior value 460 DO jj = 2, jpjm1 461 DO ji = fs_2, fs_jpim1 462 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 463 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux 464 END DO 465 END DO 466 END DO 467 468 jtaken = MOD( jtaken + 1 , 2 ) 469 470 DO jk = 2, jpkm1 ! Interior value 471 DO jj = 2, jpjm1 472 DO ji = fs_2, fs_jpim1 473 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 474 ! total advective trends 475 ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 476 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 477 END DO 478 END DO 479 END DO 480 481 END DO 482 483 DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping 484 DO jj = 2, jpjm1 485 DO ji = fs_2, fs_jpim1 486 zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 487 END DO 488 END DO 489 END DO 490 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 491 CALL lbc_lnk( zwz, 'W', 1. ) 492 493 ! 4. monotonicity algorithm 494 ! ------------------------- 495 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 496 497 498 ! 5. final trend with corrected fluxes 499 ! ------------------------------------ 500 DO jk = 1, jpkm1 501 DO jj = 2, jpjm1 502 DO ji = fs_2, fs_jpim1 ! vector opt. 503 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 504 ! total advective trends 505 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 506 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 507 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 508 ! add them to the general tracer trends 509 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 510 END DO 511 END DO 512 END DO 513 514 ! ! trend diagnostics (contribution of upstream fluxes) 515 IF( l_trd ) THEN 516 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 517 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 518 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 519 520 CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) ) 521 CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 522 CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 523 END IF 524 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 525 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 526 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 527 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 528 ENDIF 529 ! 530 END DO 531 ! 532 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 533 CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 534 CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 535 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 536 ! 537 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts') 538 ! 539 END SUBROUTINE tra_adv_tvd_zts 264 540 265 541 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
Note: See TracChangeset
for help on using the changeset viewer.
