Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r4499 r5260 22 22 USE oce ! ocean dynamics and active tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE trdmod_oce ! tracers trends 24 USE trc_oce ! share passive tracers/Ocean variables 25 USE trd_oce ! trends: ocean variables 25 26 USE trdtra ! tracers trends 26 USE in_out_manager ! I/O manager27 27 USE dynspg_oce ! choice/control of key cpp for surface pressure gradient 28 USE diaptr ! poleward transport diagnostics 29 ! 28 30 USE lib_mpp ! MPP library 29 31 USE lbclnk ! ocean lateral boundary condition (or mpp link) 30 USE diaptr ! poleward transport diagnostics 31 USE trc_oce ! share passive tracers/Ocean variables 32 USE in_out_manager ! I/O manager 32 33 USE wrk_nemo ! Memory Allocation 33 34 USE timing ! Timing … … 37 38 PRIVATE 38 39 39 PUBLIC tra_adv_tvd ! routine called by step.F90 40 PUBLIC tra_adv_tvd ! routine called by traadv.F90 41 PUBLIC tra_adv_tvd_zts ! routine called by traadv.F90 40 42 41 43 LOGICAL :: l_trd ! flag to compute trends … … 77 79 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 78 80 ! 79 INTEGER :: ji, jj, jk, jn ! dummy loop indices 81 INTEGER :: ji, jj, jk, jn ! dummy loop indices 82 INTEGER :: ik 80 83 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 81 84 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - … … 93 96 IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype 94 97 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 98 ! 99 l_trd = .FALSE. 100 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 95 101 ENDIF 96 !97 l_trd = .FALSE.98 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.99 102 ! 100 103 IF( l_trd ) THEN … … 103 106 ENDIF 104 107 ! 105 zwi(:,:,:) = 0.e0 108 zwi(:,:,:) = 0.e0 ; 106 109 ! 107 110 ! ! =========== 108 111 DO jn = 1, kjpt ! tracer loop 109 112 ! ! =========== 110 ! 1. Bottom value : flux set to zero113 ! 1. Bottom and k=1 value : flux set to zero 111 114 ! ---------------------------------- 112 115 zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0 113 116 zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0 114 117 118 zwz(:,:,1 ) = 0._wp 115 119 ! 2. upstream advection with initial mass fluxes & intermediate update 116 120 ! -------------------------------------------------------------------- … … 131 135 132 136 ! upstream tracer flux in the k direction 133 ! Surface value134 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable135 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface136 ENDIF137 137 ! Interior value 138 138 DO jk = 2, jpkm1 … … 141 141 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 142 142 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 143 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 144 END DO 145 END DO 146 END DO 143 zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk) 144 END DO 145 END DO 146 END DO 147 ! Surface value 148 IF( lk_vvl ) THEN 149 IF ( ln_isfcav ) THEN 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 153 END DO 154 END DO 155 ELSE 156 zwz(:,:,1) = 0.e0 ! volume variable 157 END IF 158 ELSE 159 IF ( ln_isfcav ) THEN 160 DO jj = 1, jpj 161 DO ji = 1, jpi 162 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 163 END DO 164 END DO 165 ELSE 166 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 167 END IF 168 ENDIF 147 169 148 170 ! total advective trend … … 157 179 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 158 180 ! update and guess with monotonic sheme 159 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 181 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 160 182 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 161 183 END DO … … 171 193 END IF 172 194 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 173 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN174 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) )175 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) )195 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 196 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 197 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 176 198 ENDIF 177 199 … … 189 211 190 212 ! antidiffusive flux on k 191 zwz(:,:,1) = 0.e0 ! Surface value 192 ! 193 DO jk = 2, jpkm1 ! Interior value 213 ! Interior value 214 DO jk = 2, jpkm1 194 215 DO jj = 1, jpj 195 216 DO ji = 1, jpi 196 217 zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk) 218 END DO 219 END DO 220 END DO 221 ! surface value 222 IF ( ln_isfcav ) THEN 223 DO jj = 1, jpj 224 DO ji = 1, jpi 225 zwz(ji,jj,mikt(ji,jj)) = 0.e0 226 END DO 227 END DO 228 ELSE 229 zwz(:,:,1) = 0.e0 230 END IF 231 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 232 CALL lbc_lnk( zwz, 'W', 1. ) 233 234 ! 4. monotonicity algorithm 235 ! ------------------------- 236 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 237 238 239 ! 5. final trend with corrected fluxes 240 ! ------------------------------------ 241 DO jk = 1, jpkm1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 245 ! total advective trends 246 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 247 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 248 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 249 ! add them to the general tracer trends 250 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 251 END DO 252 END DO 253 END DO 254 255 ! ! trend diagnostics (contribution of upstream fluxes) 256 IF( l_trd ) THEN 257 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 258 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 259 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 260 261 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 262 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 263 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 264 END IF 265 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 266 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 267 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 268 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 269 ENDIF 270 ! 271 END DO 272 ! 273 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 274 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 275 ! 276 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd') 277 ! 278 END SUBROUTINE tra_adv_tvd 279 280 SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 281 & ptb, ptn, pta, kjpt ) 282 !!---------------------------------------------------------------------- 283 !! *** ROUTINE tra_adv_tvd_zts *** 284 !! 285 !! ** Purpose : Compute the now trend due to total advection of 286 !! tracers and add it to the general trend of tracer equations 287 !! 288 !! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with 289 !! corrected flux (monotonic correction). This version use sub- 290 !! timestepping for the vertical advection which increases stability 291 !! when vertical metrics are small. 292 !! note: - this advection scheme needs a leap-frog time scheme 293 !! 294 !! ** Action : - update (pta) with the now advective tracer trends 295 !! - save the trends 296 !!---------------------------------------------------------------------- 297 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 298 ! 299 INTEGER , INTENT(in ) :: kt ! ocean time-step index 300 INTEGER , INTENT(in ) :: kit000 ! first time step index 301 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 302 INTEGER , INTENT(in ) :: kjpt ! number of tracers 303 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 304 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 305 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 306 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 307 ! 308 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 309 REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep 310 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 311 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 312 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 313 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 314 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 315 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 316 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 317 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 318 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav 319 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 320 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 321 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 322 !!---------------------------------------------------------------------- 323 ! 324 IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd_zts') 325 ! 326 CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 327 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 328 CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 329 ! 330 IF( kt == kit000 ) THEN 331 IF(lwp) WRITE(numout,*) 332 IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 333 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 334 ENDIF 335 ! 336 l_trd = .FALSE. 337 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 338 ! 339 IF( l_trd ) THEN 340 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 341 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 342 ENDIF 343 ! 344 zwi(:,:,:) = 0._wp 345 z_rzts = 1._wp / REAL( jnzts, wp ) 346 zr_p2dt(:) = 1._wp / p2dt(:) 347 ! 348 ! ! =========== 349 DO jn = 1, kjpt ! tracer loop 350 ! ! =========== 351 ! 1. Bottom value : flux set to zero 352 ! ---------------------------------- 353 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 354 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 355 356 ! 2. upstream advection with initial mass fluxes & intermediate update 357 ! -------------------------------------------------------------------- 358 ! upstream tracer flux in the i and j direction 359 DO jk = 1, jpkm1 360 DO jj = 1, jpjm1 361 DO ji = 1, fs_jpim1 ! vector opt. 362 ! upstream scheme 363 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 364 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 365 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 366 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 367 zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 368 zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 369 END DO 370 END DO 371 END DO 372 373 ! upstream tracer flux in the k direction 374 ! Interior value 375 DO jk = 2, jpkm1 376 DO jj = 1, jpj 377 DO ji = 1, jpi 378 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 379 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 380 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 381 END DO 382 END DO 383 END DO 384 ! Surface value 385 IF( lk_vvl ) THEN 386 IF ( ln_isfcav ) THEN 387 DO jj = 1, jpj 388 DO ji = 1, jpi 389 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable + isf 390 END DO 391 END DO 392 ELSE 393 zwz(:,:,1) = 0.e0 ! volume variable + no isf 394 END IF 395 ELSE 396 IF ( ln_isfcav ) THEN 397 DO jj = 1, jpj 398 DO ji = 1, jpi 399 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface + isf 400 END DO 401 END DO 402 ELSE 403 zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface + no isf 404 END IF 405 ENDIF 406 407 ! total advective trend 408 DO jk = 1, jpkm1 409 z2dtt = p2dt(jk) 410 DO jj = 2, jpjm1 411 DO ji = fs_2, fs_jpim1 ! vector opt. 412 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 413 ! total intermediate advective trends 414 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 415 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 416 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 417 ! update and guess with monotonic sheme 418 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 419 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 420 END DO 421 END DO 422 END DO 423 ! ! Lateral boundary conditions on zwi (unchanged sign) 424 CALL lbc_lnk( zwi, 'T', 1. ) 425 426 ! ! trend diagnostics (contribution of upstream fluxes) 427 IF( l_trd ) THEN 428 ! store intermediate advective trends 429 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 430 END IF 431 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 432 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 433 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) 434 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) 435 ENDIF 436 437 ! 3. antidiffusive flux : high order minus low order 438 ! -------------------------------------------------- 439 ! antidiffusive flux on i and j 440 441 442 DO jk = 1, jpkm1 443 444 DO jj = 1, jpjm1 445 DO ji = 1, fs_jpim1 ! vector opt. 446 zwx_sav(ji,jj) = zwx(ji,jj,jk) 447 zwy_sav(ji,jj) = zwy(ji,jj,jk) 448 449 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 450 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 451 END DO 452 END DO 453 454 DO jj = 2, jpjm1 ! partial horizontal divergence 455 DO ji = fs_2, fs_jpim1 456 zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 457 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 458 END DO 459 END DO 460 461 DO jj = 1, jpjm1 462 DO ji = 1, fs_jpim1 ! vector opt. 463 zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 464 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 465 END DO 466 END DO 467 END DO 468 469 ! antidiffusive flux on k 470 zwz(:,:,1) = 0._wp ! Surface value 471 zwz_sav(:,:,:) = zwz(:,:,:) 472 ! 473 ztrs(:,:,:,1) = ptb(:,:,:,jn) 474 zwzts(:,:,:) = 0._wp 475 476 DO jl = 1, jnzts ! Start of sub timestepping loop 477 478 IF( jl == 1 ) THEN ! Euler forward to kick things off 479 jtb = 1 ; jtn = 1 ; jta = 2 480 zts(:) = p2dt(:) * z_rzts 481 jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux 482 ! starting at jl =1 if jnzts is odd; 483 ! starting at jl =2 otherwise 484 ELSEIF( jl == 2 ) THEN ! First leapfrog step 485 jtb = 1 ; jtn = 2 ; jta = 3 486 zts(:) = 2._wp * p2dt(:) * z_rzts 487 ELSE ! Shuffle pointers for subsequent leapfrog steps 488 jtb = MOD(jtb,3) + 1 489 jtn = MOD(jtn,3) + 1 490 jta = MOD(jta,3) + 1 491 ENDIF 492 DO jk = 2, jpkm1 ! Interior value 493 DO jj = 2, jpjm1 494 DO ji = fs_2, fs_jpim1 495 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 496 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux 497 END DO 498 END DO 499 END DO 500 501 jtaken = MOD( jtaken + 1 , 2 ) 502 503 DO jk = 2, jpkm1 ! Interior value 504 DO jj = 2, jpjm1 505 DO ji = fs_2, fs_jpim1 506 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 507 ! total advective trends 508 ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 509 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 510 END DO 511 END DO 512 END DO 513 514 END DO 515 516 DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping 517 DO jj = 2, jpjm1 518 DO ji = fs_2, fs_jpim1 519 zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 197 520 END DO 198 521 END DO … … 228 551 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 229 552 230 CALL trd_tra( kt, cdtype, jn, jptra_ trd_xad, ztrdx, pun, ptn(:,:,:,jn) )231 CALL trd_tra( kt, cdtype, jn, jptra_ trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )232 CALL trd_tra( kt, cdtype, jn, jptra_ trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )553 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 554 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 555 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 233 556 END IF 234 557 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 235 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 )) THEN236 IF( jn == jp_tem ) htr_adv(:) = ptr_ vj( zwy(:,:,:) ) + htr_adv(:)237 IF( jn == jp_sal ) str_adv(:) = ptr_ vj( zwy(:,:,:) ) + str_adv(:)558 IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN 559 IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:) 560 IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:) 238 561 ENDIF 239 562 ! 240 563 END DO 241 564 ! 242 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 565 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 566 CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 567 CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 243 568 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 244 569 ! 245 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd') 246 ! 247 END SUBROUTINE tra_adv_tvd 248 570 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts') 571 ! 572 END SUBROUTINE tra_adv_tvd_zts 249 573 250 574 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) … … 261 585 !! in-space based differencing for fluid 262 586 !!---------------------------------------------------------------------- 263 !264 !!----------------------------------------------------------------------265 587 REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 266 588 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 267 589 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 268 590 ! 269 INTEGER :: ji, jj, jk ! dummy loop indices270 INTEGER :: ikm1 ! local integer591 INTEGER :: ji, jj, jk ! dummy loop indices 592 INTEGER :: ikm1 ! local integer 271 593 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars 272 594 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - … … 278 600 CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 279 601 ! 280 281 602 zbig = 1.e+40_wp 282 603 zrtrn = 1.e-15_wp 283 zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp 284 604 zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp 285 605 286 606 ! Search local extrema 287 607 ! -------------------- 288 608 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 289 zbup = MAX( pbef * tmask - zbig * ( 1. e0- tmask ), &290 & paft * tmask - zbig * ( 1. e0- tmask ) )291 zbdo = MIN( pbef * tmask + zbig * ( 1. e0- tmask ), &292 & paft * tmask + zbig * ( 1. e0- tmask ) )609 zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), & 610 & paft * tmask - zbig * ( 1._wp - tmask ) ) 611 zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), & 612 & paft * tmask + zbig * ( 1._wp - tmask ) ) 293 613 294 614 DO jk = 1, jpkm1 … … 334 654 DO jj = 2, jpjm1 335 655 DO ji = fs_2, fs_jpim1 ! vector opt. 336 zau = MIN( 1. e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )337 zbu = MIN( 1. e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )656 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 657 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 338 658 zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) ) 339 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1. e0- zcu) * zbu )340 341 zav = MIN( 1. e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )342 zbv = MIN( 1. e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )659 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 660 661 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 662 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 343 663 zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 344 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1. e0- zcv) * zbv )664 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 345 665 346 666 ! monotonic flux in the k direction, i.e. pcc … … 349 669 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 350 670 zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) ) 351 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1. e0- zc) * zb )671 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 352 672 END DO 353 673 END DO
Note: See TracChangeset
for help on using the changeset viewer.