Changeset 4990 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
- Timestamp:
- 2014-12-15T17:42:49+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r4499 r4990 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 ; zwz(:,:,:) = 0.e0 106 109 ! 107 110 ! ! =========== … … 132 135 ! upstream tracer flux in the k direction 133 136 ! Surface value 134 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable 135 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 137 IF( lk_vvl ) THEN 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 zwz(ji,jj, mikt(ji,jj) ) = 0.e0 ! volume variable 141 END DO 142 END DO 143 ELSE 144 DO jj = 1, jpj 145 DO ji = 1, jpi 146 zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface 147 END DO 148 END DO 136 149 ENDIF 137 150 ! Interior value 138 DO j k = 2, jpkm1139 DO j j = 1, jpj140 DO j i = 1, jpi151 DO jj = 1, jpj 152 DO ji = 1, jpi 153 DO jk = mikt(ji,jj)+1, jpkm1 141 154 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 142 155 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) … … 157 170 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 158 171 ! update and guess with monotonic sheme 159 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 172 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 160 173 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 161 174 END DO … … 191 204 zwz(:,:,1) = 0.e0 ! Surface value 192 205 ! 193 DO jk = 2, jpkm1 ! Interior value 206 DO jj = 1, jpj 207 DO ji = 1, jpi 208 ik=mikt(ji,jj) 209 ! surface value 210 zwz(ji,jj,1:ik) = 0.e0 211 ! Interior value 212 DO jk = mikt(ji,jj)+1, jpkm1 213 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) 214 END DO 215 END DO 216 END DO 217 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions 218 CALL lbc_lnk( zwz, 'W', 1. ) 219 220 ! 4. monotonicity algorithm 221 ! ------------------------- 222 CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt ) 223 224 225 ! 5. final trend with corrected fluxes 226 ! ------------------------------------ 227 DO jk = 1, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 231 ! total advective trends 232 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 233 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 234 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 235 ! add them to the general tracer trends 236 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra * tmask(ji,jj,jk) 237 END DO 238 END DO 239 END DO 240 241 ! ! trend diagnostics (contribution of upstream fluxes) 242 IF( l_trd ) THEN 243 ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed 244 ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed 245 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 246 247 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 248 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 249 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 250 END IF 251 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 252 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 253 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:) 254 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:) 255 ENDIF 256 ! 257 END DO 258 ! 259 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 260 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 261 ! 262 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd') 263 ! 264 END SUBROUTINE tra_adv_tvd 265 266 SUBROUTINE tra_adv_tvd_zts ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, & 267 & ptb, ptn, pta, kjpt ) 268 !!---------------------------------------------------------------------- 269 !! *** ROUTINE tra_adv_tvd_zts *** 270 !! 271 !! ** Purpose : Compute the now trend due to total advection of 272 !! tracers and add it to the general trend of tracer equations 273 !! 274 !! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with 275 !! corrected flux (monotonic correction). This version use sub- 276 !! timestepping for the vertical advection which increases stability 277 !! when vertical metrics are small. 278 !! note: - this advection scheme needs a leap-frog time scheme 279 !! 280 !! ** Action : - update (pta) with the now advective tracer trends 281 !! - save the trends 282 !!---------------------------------------------------------------------- 283 USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace 284 ! 285 INTEGER , INTENT(in ) :: kt ! ocean time-step index 286 INTEGER , INTENT(in ) :: kit000 ! first time step index 287 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 288 INTEGER , INTENT(in ) :: kjpt ! number of tracers 289 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 290 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components 291 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields 292 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 293 ! 294 REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection 295 REAL(wp), DIMENSION( jpk ) :: zr_p2dt ! reciprocal of tracer timestep 296 INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices 297 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 298 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 299 INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps 300 REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection 301 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 302 REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - - 303 REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - - 304 REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav 305 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwz, zhdiv, zwz_sav, zwzts 306 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz 307 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs 308 !!---------------------------------------------------------------------- 309 ! 310 IF( nn_timing == 1 ) CALL timing_start('tra_adv_tvd_zts') 311 ! 312 CALL wrk_alloc( jpi, jpj, zwx_sav, zwy_sav ) 313 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwz , zhdiv, zwz_sav, zwzts ) 314 CALL wrk_alloc( jpi, jpj, jpk, 3, ztrs ) 315 ! 316 IF( kt == kit000 ) THEN 317 IF(lwp) WRITE(numout,*) 318 IF(lwp) WRITE(numout,*) 'tra_adv_tvd_zts : TVD ZTS advection scheme on ', cdtype 319 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 320 ENDIF 321 ! 322 l_trd = .FALSE. 323 IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE. 324 ! 325 IF( l_trd ) THEN 326 CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 327 ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp 328 ENDIF 329 ! 330 zwi(:,:,:) = 0._wp 331 z_rzts = 1._wp / REAL( jnzts, wp ) 332 zr_p2dt(:) = 1._wp / p2dt(:) 333 ! 334 ! ! =========== 335 DO jn = 1, kjpt ! tracer loop 336 ! ! =========== 337 ! 1. Bottom value : flux set to zero 338 ! ---------------------------------- 339 zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp 340 zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp 341 342 ! 2. upstream advection with initial mass fluxes & intermediate update 343 ! -------------------------------------------------------------------- 344 ! upstream tracer flux in the i and j direction 345 DO jk = 1, jpkm1 346 DO jj = 1, jpjm1 347 DO ji = 1, fs_jpim1 ! vector opt. 348 ! upstream scheme 349 zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) ) 350 zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) ) 351 zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) ) 352 zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) ) 353 zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) ) 354 zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) ) 355 END DO 356 END DO 357 END DO 358 359 ! upstream tracer flux in the k direction 360 ! Surface value 361 IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0._wp ! volume variable 362 ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface 363 ENDIF 364 ! Interior value 365 DO jk = 2, jpkm1 194 366 DO jj = 1, jpj 195 367 DO ji = 1, jpi 196 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) 368 zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) ) 369 zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) ) 370 zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) 371 END DO 372 END DO 373 END DO 374 375 ! total advective trend 376 DO jk = 1, jpkm1 377 z2dtt = p2dt(jk) 378 DO jj = 2, jpjm1 379 DO ji = fs_2, fs_jpim1 ! vector opt. 380 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 381 ! total intermediate advective trends 382 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 383 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) & 384 & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 385 ! update and guess with monotonic sheme 386 pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra 387 zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk) 388 END DO 389 END DO 390 END DO 391 ! ! Lateral boundary conditions on zwi (unchanged sign) 392 CALL lbc_lnk( zwi, 'T', 1. ) 393 394 ! ! trend diagnostics (contribution of upstream fluxes) 395 IF( l_trd ) THEN 396 ! store intermediate advective trends 397 ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:) 398 END IF 399 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) 400 IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN 401 IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) 402 IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) 403 ENDIF 404 405 ! 3. antidiffusive flux : high order minus low order 406 ! -------------------------------------------------- 407 ! antidiffusive flux on i and j 408 409 410 DO jk = 1, jpkm1 411 412 DO jj = 1, jpjm1 413 DO ji = 1, fs_jpim1 ! vector opt. 414 zwx_sav(ji,jj) = zwx(ji,jj,jk) 415 zwy_sav(ji,jj) = zwy(ji,jj,jk) 416 417 zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) 418 zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) 419 END DO 420 END DO 421 422 DO jj = 2, jpjm1 ! partial horizontal divergence 423 DO ji = fs_2, fs_jpim1 424 zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) & 425 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) ) 426 END DO 427 END DO 428 429 DO jj = 1, jpjm1 430 DO ji = 1, fs_jpim1 ! vector opt. 431 zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj) 432 zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj) 433 END DO 434 END DO 435 END DO 436 437 ! antidiffusive flux on k 438 zwz(:,:,1) = 0._wp ! Surface value 439 zwz_sav(:,:,:) = zwz(:,:,:) 440 ! 441 ztrs(:,:,:,1) = ptb(:,:,:,jn) 442 zwzts(:,:,:) = 0._wp 443 444 DO jl = 1, jnzts ! Start of sub timestepping loop 445 446 IF( jl == 1 ) THEN ! Euler forward to kick things off 447 jtb = 1 ; jtn = 1 ; jta = 2 448 zts(:) = p2dt(:) * z_rzts 449 jtaken = MOD( jnzts + 1 , 2) ! Toggle to collect every second flux 450 ! starting at jl =1 if jnzts is odd; 451 ! starting at jl =2 otherwise 452 ELSEIF( jl == 2 ) THEN ! First leapfrog step 453 jtb = 1 ; jtn = 2 ; jta = 3 454 zts(:) = 2._wp * p2dt(:) * z_rzts 455 ELSE ! Shuffle pointers for subsequent leapfrog steps 456 jtb = MOD(jtb,3) + 1 457 jtn = MOD(jtn,3) + 1 458 jta = MOD(jta,3) + 1 459 ENDIF 460 DO jk = 2, jpkm1 ! Interior value 461 DO jj = 2, jpjm1 462 DO ji = fs_2, fs_jpim1 463 zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) 464 IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk)*zts(jk) ! Accumulate time-weighted vertcal flux 465 END DO 466 END DO 467 END DO 468 469 jtaken = MOD( jtaken + 1 , 2 ) 470 471 DO jk = 2, jpkm1 ! Interior value 472 DO jj = 2, jpjm1 473 DO ji = fs_2, fs_jpim1 474 zbtr = 1._wp / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) ) 475 ! total advective trends 476 ztra = - zbtr * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) 477 ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) + zts(jk) * ztra 478 END DO 479 END DO 480 END DO 481 482 END DO 483 484 DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping 485 DO jj = 2, jpjm1 486 DO ji = fs_2, fs_jpim1 487 zwz(ji,jj,jk) = zwzts(ji,jj,jk) * zr_p2dt(jk) - zwz_sav(ji,jj,jk) 197 488 END DO 198 489 END DO … … 228 519 ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed 229 520 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) )521 CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) ) 522 CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) ) 523 CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) ) 233 524 END IF 234 525 ! ! "Poleward" heat and salt transports (contribution of upstream fluxes) … … 240 531 END DO 241 532 ! 242 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz ) 533 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwz, zhdiv, zwz_sav, zwzts ) 534 CALL wrk_dealloc( jpi, jpj, jpk, 3, ztrs ) 535 CALL wrk_dealloc( jpi, jpj, zwx_sav, zwy_sav ) 243 536 IF( l_trd ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz ) 244 537 ! 245 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd') 246 ! 247 END SUBROUTINE tra_adv_tvd 248 538 IF( nn_timing == 1 ) CALL timing_stop('tra_adv_tvd_zts') 539 ! 540 END SUBROUTINE tra_adv_tvd_zts 249 541 250 542 SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt ) … … 261 553 !! in-space based differencing for fluid 262 554 !!---------------------------------------------------------------------- 263 !264 !!----------------------------------------------------------------------265 555 REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 266 556 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 267 557 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions 268 558 ! 269 INTEGER :: ji, jj, jk ! dummy loop indices270 INTEGER :: ikm1 ! local integer559 INTEGER :: ji, jj, jk ! dummy loop indices 560 INTEGER :: ikm1 ! local integer 271 561 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars 272 562 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - … … 278 568 CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 279 569 ! 280 281 570 zbig = 1.e+40_wp 282 571 zrtrn = 1.e-15_wp 283 zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp 284 572 zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp 285 573 286 574 ! Search local extrema 287 575 ! -------------------- 288 576 ! 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 ) )293 294 DO j k = 1, jpkm1295 ikm1 = MAX(jk-1,1)296 z2dtt = p2dt(jk)297 DO jj = 2, jpjm1298 DO ji = fs_2, fs_jpim1 ! vector opt.299 577 zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), & 578 & paft * tmask - zbig * ( 1._wp - tmask ) ) 579 zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), & 580 & paft * tmask + zbig * ( 1._wp - tmask ) ) 581 582 DO jj = 2, jpjm1 583 DO ji = fs_2, fs_jpim1 ! vector opt. 584 DO jk = mikt(ji,jj), jpkm1 585 ikm1 = MAX(jk-1,mikt(ji,jj)) 586 z2dtt = p2dt(jk) 587 300 588 ! search maximum in neighbourhood 301 589 zup = MAX( zbup(ji ,jj ,jk ), & … … 334 622 DO jj = 2, jpjm1 335 623 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) )624 zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 625 zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 338 626 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) )627 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu ) 628 629 zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 630 zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 343 631 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 )632 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv ) 345 633 346 634 ! monotonic flux in the k direction, i.e. pcc … … 349 637 zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) ) 350 638 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 )639 pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb ) 352 640 END DO 353 641 END DO
Note: See TracChangeset
for help on using the changeset viewer.