Changeset 10446
- Timestamp:
- 2018-12-21T16:07:19+01:00 (4 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r10439 r10446 35 35 REAL(wp) :: z1_120 = 1._wp / 120._wp ! =1/120 36 36 37 ! limiter: 1=nonosc, 2=superbee, 3=h3(rachid) 38 INTEGER :: kn_limiter = 1 39 37 40 ! if T interpolated at u/v points is negative, then interpolate T at u/v points using the upstream scheme 38 LOGICAL :: ll_neg = .TRUE.41 LOGICAL :: ll_neg = .TRUE. 39 42 40 43 ! alternate directions for upstream 41 LOGICAL :: ll_upsxy = .TRUE.44 LOGICAL :: ll_upsxy = .TRUE. 42 45 43 46 ! alternate directions for high order 44 LOGICAL :: ll_hoxy = .TRUE.47 LOGICAL :: ll_hoxy = .TRUE. 45 48 46 49 ! prelimiter: use it to avoid overshoot in H 47 LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D50 LOGICAL :: ll_prelimiter_zalesak = .TRUE. ! from: Zalesak(1979) eq. 14 => better for 1D. Not well defined in 2D 48 51 49 52 … … 213 216 !! 214 217 !! ** Purpose : Compute the now trend due to total advection of 215 !! tracers and add it to the general trend of tracer equations218 !! tracers and add it to the general trend of tracer equations 216 219 !! 217 !! ** Method : TVD scheme, i.e. 2nd order centered scheme with 218 !! corrected flux (monotonic correction) 219 !! note: - this advection scheme needs a leap-frog time scheme 220 !! ** Method : - calculate upstream fluxes and upstream solution for tracer H 221 !! - calculate tracer H at u and v points (Ultimate) 222 !! - calculate the high order fluxes using alterning directions (Macho?) 223 !! - apply a limiter on the fluxes (nonosc) 224 !! - convert this tracer flux to a tracer content flux (uH -> uV) 225 !! - calculate the high order solution for tracer content V 220 226 !! 221 !! ** Action : - pt the after advective tracer 227 !! ** Action : solve 2 equations => a) da/dt = -div(ua) 228 !! b) dV/dt = -div(uV) using dH/dt = -u.grad(H) 229 !! in eq. b), - fluxes uH are evaluated (with UMx) and limited (with nonosc). This step is necessary to get a good H. 230 !! - then we convert this flux to a "volume" flux this way => uH*ua/u 231 !! where ua is the flux from eq. a) 232 !! - at last we estimate dV/dt = -div(uH*ua/u) 233 !! 234 !! ** Note : - this method can lead to small negative V (since we only limit H) => corrected in icedyn_adv.F90 conserving mass etc. 235 !! - negative tracers at u-v points can also occur from the Ultimate scheme (usually at the ice edge) and the solution for now 236 !! is to apply an upstream scheme when it occurs. A better solution would be to degrade the order of 237 !! the scheme automatically by applying a mask of the ice cover inside Ultimate (not done). 238 !! - Eventhough 1D tests give very good results (typically the one from Schar & Smolarkiewiecz), the 2D is less good. 239 !! Large values of H can appear for very small ice concentration, and when it does it messes the things up since we 240 !! work on H (and not V). It probably comes from the prelimiter of zalesak which is coded for 1D and not 2D. 241 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 242 !! concentration is small). 243 !! To-do: expand the prelimiter from zalesak to make it work in 2D 222 244 !!---------------------------------------------------------------------- 223 245 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 235 257 INTEGER :: ji, jj, jl ! dummy loop indices 236 258 REAL(wp) :: ztra ! local scalar 237 INTEGER :: kn_limiter = 1 ! 1=nonosc ; 2=superbee ; 3=h3(rachid) 238 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zt_u, zt_v, zpt 259 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ho , zfv_ho , zpt 239 260 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zfu_ups, zfv_ups, zt_ups 240 261 !!---------------------------------------------------------------------- 241 262 ! 242 ! upstream (_ups) advection with initial mass fluxes 243 ! --------------------------------------------------- 244 ! 245 IF( .NOT. ll_upsxy ) THEN !** no alternate directions **! 246 DO jl = 1, jpl 247 DO jj = 1, jpjm1 248 DO ji = 1, fs_jpim1 249 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 250 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 251 END DO 252 END DO 253 END DO 254 255 ELSE !** alternate directions **! 256 ! 257 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 258 ! 259 DO jl = 1, jpl !-- flux in x-direction 260 DO jj = 1, jpjm1 261 DO ji = 1, fs_jpim1 262 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 263 END DO 264 END DO 265 END DO 266 ! 267 DO jl = 1, jpl !-- first guess of tracer from u-flux 268 DO jj = 2, jpjm1 269 DO ji = fs_2, fs_jpim1 270 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 271 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) & 272 & ) * tmask(ji,jj,1) 273 END DO 274 END DO 275 END DO 276 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 277 ! 278 DO jl = 1, jpl !-- flux in y-direction 279 DO jj = 1, jpjm1 280 DO ji = 1, fs_jpim1 281 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 282 END DO 283 END DO 284 END DO 285 ! 286 ELSE !== even ice time step: adv_y then adv_x ==! 287 ! 288 DO jl = 1, jpl !-- flux in y-direction 289 DO jj = 1, jpjm1 290 DO ji = 1, fs_jpim1 291 zfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 292 END DO 293 END DO 294 END DO 295 ! 296 DO jl = 1, jpl !-- first guess of tracer from v-flux 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 299 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( zfv_ups(ji,jj,jl) - zfv_ups(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 300 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 301 & * tmask(ji,jj,1) 302 END DO 303 END DO 304 END DO 305 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 306 ! 307 DO jl = 1, jpl !-- flux in x-direction 308 DO jj = 1, jpjm1 309 DO ji = 1, fs_jpim1 310 zfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 311 END DO 312 END DO 313 END DO 314 ! 315 ENDIF 316 317 ENDIF 318 ! 319 DO jl = 1, jpl !-- after tracer with upstream scheme 320 DO jj = 2, jpjm1 321 DO ji = fs_2, fs_jpim1 322 ztra = - ( zfu_ups(ji,jj,jl) - zfu_ups(ji-1,jj ,jl) & 323 & + zfv_ups(ji,jj,jl) - zfv_ups(ji ,jj-1,jl) ) * r1_e1e2t(ji,jj) 324 zt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + pdt * ztra + ( pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) & 325 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) ) & 326 & * r1_e1e2t(ji,jj) * (1.-pamsk) ) * tmask(ji,jj,1) 327 END DO 328 END DO 329 END DO 330 CALL lbc_lnk( 'icedyn_adv_umx', zt_ups, 'T', 1. ) 331 263 ! Upstream (_ups) fluxes 264 ! ----------------------- 265 CALL upstream( pamsk, jt, kt, pdt, pt, pu, pv, zt_ups, zfu_ups, zfv_ups ) 266 332 267 ! High order (_ho) fluxes 333 268 ! ----------------------- … … 336 271 CASE ( 20 ) !== centered second order ==! 337 272 ! 338 CALL cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zfu_ho, zfv_ho, & 339 & zt_ups, zfu_ups, zfv_ups ) 273 CALL cen2( pamsk, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 340 274 ! 341 275 CASE ( 1:5 ) !== 1st to 5th order ULTIMATE-MACHO scheme ==! 342 276 ! 343 CALL macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_u, zt_v, zfu_ho, zfv_ho, & 344 & zt_ups, zfu_ups, zfv_ups ) 277 CALL macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, zt_ups, zfu_ups, zfv_ups, zfu_ho, zfv_ho ) 345 278 ! 346 279 END SELECT … … 372 305 END DO 373 306 ENDIF 374 ! 375 ! in case of advection of A: output u*a (ho)376 ! ------------------------------------- ----307 ! --ho 308 ! in case of advection of A: output u*a 309 ! ------------------------------------- 377 310 IF( PRESENT( pua_ho ) ) THEN 378 311 DO jl = 1, jpl … … 391 324 DO jj = 2, jpjm1 392 325 DO ji = fs_2, fs_jpim1 393 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) * r1_e1e2t(ji,jj) * pdt394 395 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra ) * tmask(ji,jj,1)326 ztra = - ( zfu_ho(ji,jj,jl) - zfu_ho(ji-1,jj,jl) + zfv_ho(ji,jj,jl) - zfv_ho(ji,jj-1,jl) ) 327 ! 328 ptc(ji,jj,jl) = ( ptc(ji,jj,jl) + ztra * r1_e1e2t(ji,jj) * pdt ) * tmask(ji,jj,1) 396 329 END DO 397 330 END DO … … 402 335 403 336 404 SUBROUTINE cen2( pamsk, kn_limiter, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pfu_ho, pfv_ho, & 405 & pt_ups, pfu_ups, pfv_ups ) 337 SUBROUTINE upstream( pamsk, jt, kt, pdt, pt, pu, pv, pt_ups, pfu_ups, pfv_ups ) 406 338 !!--------------------------------------------------------------------- 407 !! *** ROUTINE macho***339 !! *** ROUTINE upstream *** 408 340 !! 409 !! ** Purpose : compute 410 !! 411 !! ** Method : ... ??? 412 !! TIM = transient interpolation Modeling 413 !! 414 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 341 !! ** Purpose : compute the upstream fluxes and upstream guess of tracer 415 342 !!---------------------------------------------------------------------- 416 343 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 417 INTEGER , INTENT(in ) :: kn_limiter ! limiter 344 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 345 INTEGER , INTENT(in ) :: kt ! number of iteration 346 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 347 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt ! tracer fields 348 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pu, pv ! 2 ice velocity components 349 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_ups ! upstream guess of tracer 350 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ups, pfv_ups ! upstream fluxes 351 ! 352 INTEGER :: ji, jj, jl ! dummy loop indices 353 REAL(wp) :: ztra ! local scalar 354 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt 355 !!---------------------------------------------------------------------- 356 357 IF( .NOT. ll_upsxy ) THEN !** no alternate directions **! 358 ! 359 DO jl = 1, jpl 360 DO jj = 1, jpjm1 361 DO ji = 1, fs_jpim1 362 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 363 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 364 END DO 365 END DO 366 END DO 367 ! 368 ELSE !** alternate directions **! 369 ! 370 IF( MOD( (kt - 1) / nn_fsbc , 2 ) == MOD( (jt - 1) , 2 ) ) THEN !== odd ice time step: adv_x then adv_y ==! 371 ! 372 DO jl = 1, jpl !-- flux in x-direction 373 DO jj = 1, jpjm1 374 DO ji = 1, fs_jpim1 375 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 376 END DO 377 END DO 378 END DO 379 ! 380 DO jl = 1, jpl !-- first guess of tracer from u-flux 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 383 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 384 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 385 ! 386 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 387 END DO 388 END DO 389 END DO 390 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 391 ! 392 DO jl = 1, jpl !-- flux in y-direction 393 DO jj = 1, jpjm1 394 DO ji = 1, fs_jpim1 395 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 396 END DO 397 END DO 398 END DO 399 ! 400 ELSE !== even ice time step: adv_y then adv_x ==! 401 ! 402 DO jl = 1, jpl !-- flux in y-direction 403 DO jj = 1, jpjm1 404 DO ji = 1, fs_jpim1 405 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 406 END DO 407 END DO 408 END DO 409 ! 410 DO jl = 1, jpl !-- first guess of tracer from v-flux 411 DO jj = 2, jpjm1 412 DO ji = fs_2, fs_jpim1 413 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 414 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 415 ! 416 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 417 END DO 418 END DO 419 END DO 420 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 421 ! 422 DO jl = 1, jpl !-- flux in x-direction 423 DO jj = 1, jpjm1 424 DO ji = 1, fs_jpim1 425 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 426 END DO 427 END DO 428 END DO 429 ! 430 ENDIF 431 432 ENDIF 433 ! 434 DO jl = 1, jpl !-- after tracer with upstream scheme 435 DO jj = 2, jpjm1 436 DO ji = fs_2, fs_jpim1 437 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 438 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & 439 & + ( pu (ji,jj ) - pu (ji-1,jj ) & 440 & + pv (ji,jj ) - pv (ji ,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 441 ! 442 pt_ups(ji,jj,jl) = ( pt (ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 443 END DO 444 END DO 445 END DO 446 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1. ) 447 448 END SUBROUTINE upstream 449 450 451 SUBROUTINE cen2( pamsk, jt, kt, pdt, pt, pu, pv, puc, pvc, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 452 !!--------------------------------------------------------------------- 453 !! *** ROUTINE cen2 *** 454 !! 455 !! ** Purpose : compute the high order fluxes using a centered 456 !! second order scheme 457 !!---------------------------------------------------------------------- 458 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 418 459 INTEGER , INTENT(in ) :: jt ! number of sub-iteration 419 460 INTEGER , INTENT(in ) :: kt ! number of iteration … … 423 464 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: puc, pvc ! 2 ice velocity * A components 424 465 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step 466 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer 467 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes 425 468 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 426 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content427 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes428 469 ! 429 470 INTEGER :: ji, jj, jl ! dummy loop indices 430 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt 471 REAL(wp) :: ztra ! local scalar 472 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zpt 431 473 !!---------------------------------------------------------------------- 432 474 ! … … 442 484 END DO 443 485 IF ( kn_limiter == 1 ) THEN 444 CALL nonosc_2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 445 ELSEIF( kn_limiter == 2 ) THEN 446 CALL limiter_x( pdt, pu, pt, pfu_ho ) 447 CALL limiter_y( pdt, pv, pt, pfv_ho ) 448 ELSEIF( kn_limiter == 3 ) THEN 449 CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 450 CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 486 CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 487 ELSEIF( kn_limiter == 2 .OR. kn_limiter == 3 ) THEN 488 CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 489 CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 451 490 ENDIF 452 491 ! … … 462 501 END DO 463 502 END DO 464 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 465 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 503 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 466 504 467 505 DO jl = 1, jpl !-- first guess of tracer from u-flux 468 506 DO jj = 2, jpjm1 469 507 DO ji = fs_2, fs_jpim1 470 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) * pdt * r1_e1e2t(ji,jj) & 471 & + pt(ji,jj,jl) * pdt * ( pu(ji,jj) - pu(ji-1,jj) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 472 & * tmask(ji,jj,1) 508 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 509 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) 510 ! 511 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 473 512 END DO 474 513 END DO 475 514 END DO 476 CALL lbc_lnk( 'icedyn_adv_umx', z zt, 'T', 1. )515 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 477 516 478 517 DO jl = 1, jpl !-- flux in y-direction 479 518 DO jj = 1, jpjm1 480 519 DO ji = 1, fs_jpim1 481 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( z zt(ji,jj,jl) + zzt(ji,jj+1,jl) )520 pfv_ho(ji,jj,jl) = 0.5 * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 482 521 END DO 483 522 END DO 484 523 END DO 485 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 486 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 524 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 487 525 488 526 ELSE !== even ice time step: adv_y then adv_x ==! … … 495 533 END DO 496 534 END DO 497 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 498 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 535 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 499 536 ! 500 537 DO jl = 1, jpl !-- first guess of tracer from v-flux 501 538 DO jj = 2, jpjm1 502 539 DO ji = fs_2, fs_jpim1 503 zzt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) * pdt * r1_e1e2t(ji,jj) & 504 & + pt(ji,jj,jl) * pdt * ( pv(ji,jj) - pv(ji,jj-1) ) * r1_e1e2t(ji,jj) * (1.-pamsk) ) & 505 & * tmask(ji,jj,1) 540 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 541 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) 542 ! 543 zpt(ji,jj,jl) = ( pt(ji,jj,jl) + ztra * pdt * r1_e1e2t(ji,jj) ) * tmask(ji,jj,1) 506 544 END DO 507 545 END DO 508 546 END DO 509 CALL lbc_lnk( 'icedyn_adv_umx', z zt, 'T', 1. )547 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 510 548 ! 511 549 DO jl = 1, jpl !-- flux in x-direction 512 550 DO jj = 1, jpjm1 513 551 DO ji = 1, fs_jpim1 514 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( z zt(ji,jj,jl) + zzt(ji+1,jj,jl) )552 pfu_ho(ji,jj,jl) = 0.5 * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 515 553 END DO 516 554 END DO 517 555 END DO 518 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 519 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 556 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 520 557 521 558 ENDIF 522 IF( kn_limiter == 1 ) CALL nonosc _2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )559 IF( kn_limiter == 1 ) CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 523 560 524 561 ENDIF … … 527 564 528 565 529 SUBROUTINE macho( pamsk, kn_limiter, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_u, pt_v, pfu_ho, pfv_ho, & 530 & pt_ups, pfu_ups, pfv_ups ) 566 SUBROUTINE macho( pamsk, kn_umx, jt, kt, pdt, pt, pu, pv, puc, pvc, pubox, pvbox, ptc, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 531 567 !!--------------------------------------------------------------------- 532 568 !! *** ROUTINE macho *** 533 569 !! 534 !! ** Purpose : compute 570 !! ** Purpose : compute the high order fluxes using Ultimate-Macho scheme 535 571 !! 536 !! ** Method : ... ??? 537 !! TIM = transient interpolation Modeling 572 !! ** Method : ... 538 573 !! 539 574 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 540 575 !!---------------------------------------------------------------------- 541 576 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 542 INTEGER , INTENT(in ) :: kn_limiter ! limiter543 577 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) 544 578 INTEGER , INTENT(in ) :: jt ! number of sub-iteration … … 550 584 REAL(wp), DIMENSION(:,: ) , INTENT(in ) :: pubox, pvbox ! upstream velocity 551 585 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: ptc ! tracer content at before time step 552 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pt_u, pt_v ! tracer at u- and v-points 586 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer 587 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pfu_ups, pfv_ups ! upstream fluxes 553 588 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: pfu_ho, pfv_ho ! high order fluxes 554 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pt_ups ! upstream guess of tracer content555 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pfu_ups, pfv_ups ! upstream fluxes556 589 ! 557 590 INTEGER :: ji, jj, jl ! dummy loop indices 558 REAL(wp) :: ztra 559 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zzt, zzfu_ho, zzfv_ho 591 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zt_u, zt_v, zpt 560 592 !!---------------------------------------------------------------------- 561 593 ! … … 563 595 ! 564 596 ! !-- ultimate interpolation of pt at u-point --! 565 CALL ultimate_x( kn_umx, pdt, pt, pu, pt_u, pfu_ho )597 CALL ultimate_x( kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 566 598 ! !-- limiter in x --! 567 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 568 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 569 ! !-- advective form update in zzt --! 599 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 600 ! !-- advective form update in zpt --! 570 601 DO jl = 1, jpl 571 602 DO jj = 2, jpjm1 572 603 DO ji = fs_2, fs_jpim1 573 zzt(ji,jj,jl) = pt(ji,jj,jl) - pubox(ji,jj ) * pdt * ( pt_u(ji,jj,jl) - pt_u(ji-1,jj,jl) ) * r1_e1t(ji,jj) & 574 & - pt (ji,jj,jl) * pdt * ( pu (ji,jj) - pu (ji-1,jj) ) * r1_e1e2t(ji,jj) * pamsk 575 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 576 END DO 577 END DO 578 END DO 579 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 604 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pubox(ji,jj ) * ( zt_u(ji,jj,jl) - zt_u(ji-1,jj,jl) ) * r1_e1t (ji,jj) & 605 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & 606 & * pamsk & 607 & ) * pdt ) * tmask(ji,jj,1) 608 END DO 609 END DO 610 END DO 611 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 580 612 ! 581 613 ! !-- ultimate interpolation of pt at v-point --! 582 614 IF( ll_hoxy ) THEN 583 CALL ultimate_y( kn_umx, pdt, z zt, pv, pt_v, pfv_ho )615 CALL ultimate_y( kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 584 616 ELSE 585 CALL ultimate_y( kn_umx, pdt, pt , pv, pt_v, pfv_ho )617 CALL ultimate_y( kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 586 618 ENDIF 587 619 ! !-- limiter in y --! 588 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 589 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 620 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 590 621 ! 591 622 ! … … 593 624 ! 594 625 ! !-- ultimate interpolation of pt at v-point --! 595 CALL ultimate_y( kn_umx, pdt, pt, pv, pt_v, pfv_ho )626 CALL ultimate_y( kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 596 627 ! !-- limiter in y --! 597 IF( kn_limiter == 2 ) CALL limiter_y( pdt, pv, pt, pfv_ho ) 598 IF( kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 599 ! !-- advective form update in zzt --! 628 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 629 ! !-- advective form update in zpt --! 600 630 DO jl = 1, jpl 601 631 DO jj = 2, jpjm1 602 632 DO ji = fs_2, fs_jpim1 603 zzt(ji,jj,jl) = pt(ji,jj,jl) - pvbox(ji,jj ) * pdt * ( pt_v(ji,jj,jl) - pt_v(ji,jj-1,jl) ) * r1_e2t(ji,jj) & 604 & - pt (ji,jj,jl) * pdt * ( pv (ji,jj) - pv (ji,jj-1) ) * r1_e1e2t(ji,jj) * pamsk 605 zzt(ji,jj,jl) = zzt(ji,jj,jl) * tmask(ji,jj,1) 606 END DO 607 END DO 608 END DO 609 CALL lbc_lnk( 'icedyn_adv_umx', zzt, 'T', 1. ) 633 zpt(ji,jj,jl) = ( pt(ji,jj,jl) - ( pvbox(ji,jj ) * ( zt_v(ji,jj,jl) - zt_v(ji,jj-1,jl) ) * r1_e2t (ji,jj) & 634 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & 635 & * pamsk & 636 & ) * pdt ) * tmask(ji,jj,1) 637 END DO 638 END DO 639 END DO 640 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1. ) 610 641 ! 611 642 ! !-- ultimate interpolation of pt at u-point --! 612 643 IF( ll_hoxy ) THEN 613 CALL ultimate_x( kn_umx, pdt, z zt, pu, pt_u, pfu_ho )644 CALL ultimate_x( kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 614 645 ELSE 615 CALL ultimate_x( kn_umx, pdt, pt , pu, pt_u, pfu_ho )646 CALL ultimate_x( kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 616 647 ENDIF 617 648 ! !-- limiter in x --! 618 IF( kn_limiter == 2 ) CALL limiter_x( pdt, pu, pt, pfu_ho ) 619 IF( kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 620 ! 649 IF( kn_limiter == 2 .OR. kn_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 621 650 ! 622 651 ENDIF 623 652 624 IF( kn_limiter == 1 ) THEN 625 CALL nonosc_2d ( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 626 ENDIF 653 IF( kn_limiter == 1 ) CALL nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 627 654 ! 628 655 END SUBROUTINE macho … … 633 660 !! *** ROUTINE ultimate_x *** 634 661 !! 635 !! ** Purpose : compute 662 !! ** Purpose : compute tracer at u-points 636 663 !! 637 !! ** Method : ... ??? 638 !! TIM = transient interpolation Modeling 664 !! ** Method : ... 639 665 !! 640 666 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. … … 714 740 zdx2 = e1u(ji,jj) * e1u(ji,jj) 715 741 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 716 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) 717 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) )&718 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) 719 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ))742 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 743 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 744 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 745 & - SIGN( 1._wp, zcu ) * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 720 746 END DO 721 747 END DO … … 730 756 zdx2 = e1u(ji,jj) * e1u(ji,jj) 731 757 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 732 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) 733 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) )&734 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) 735 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ))758 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 759 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 760 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 761 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) ) 736 762 END DO 737 763 END DO … … 747 773 !!rachid zdx2 = e1u(ji,jj) * e1t(ji,jj) 748 774 zdx4 = zdx2 * zdx2 749 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) 750 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) 751 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) 752 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) 753 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) 775 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( ( pt (ji+1,jj,jl) + pt (ji,jj,jl) & 776 & - zcu * ( pt (ji+1,jj,jl) - pt (ji,jj,jl) ) ) & 777 & + z1_6 * zdx2 * ( zcu*zcu - 1._wp ) * ( ztu2(ji+1,jj,jl) + ztu2(ji,jj,jl) & 778 & - 0.5_wp * zcu * ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) ) & 779 & + z1_120 * zdx4 * ( zcu*zcu - 1._wp ) * ( zcu*zcu - 4._wp ) * ( ztu4(ji+1,jj,jl) + ztu4(ji,jj,jl) & 754 780 & - SIGN( 1._wp, zcu ) * ( ztu4(ji+1,jj,jl) - ztu4(ji,jj,jl) ) ) ) 755 781 END DO … … 790 816 !! *** ROUTINE ultimate_y *** 791 817 !! 792 !! ** Purpose : compute 818 !! ** Purpose : compute tracer at v-points 793 819 !! 794 !! ** Method : ... ??? 795 !! TIM = transient interpolation Modeling 820 !! ** Method : ... 796 821 !! 797 822 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. … … 871 896 zdy2 = e2v(ji,jj) * e2v(ji,jj) 872 897 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 873 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) 874 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) 875 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) 898 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 899 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 900 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 876 901 & - SIGN( 1._wp, zcv ) * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 877 902 END DO … … 886 911 zdy2 = e2v(ji,jj) * e2v(ji,jj) 887 912 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 888 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) 889 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) 890 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) 913 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 914 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 915 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 891 916 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) ) 892 917 END DO … … 902 927 !!rachid zdy2 = e2v(ji,jj) * e2t(ji,jj) 903 928 zdy4 = zdy2 * zdy2 904 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) 905 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) 906 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) 907 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) 908 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) 929 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt (ji,jj+1,jl) + pt (ji,jj,jl) & 930 & - zcv * ( pt (ji,jj+1,jl) - pt (ji,jj,jl) ) ) & 931 & + z1_6 * zdy2 * ( zcv*zcv - 1._wp ) * ( ztv2(ji,jj+1,jl) + ztv2(ji,jj,jl) & 932 & - 0.5_wp * zcv * ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) ) & 933 & + z1_120 * zdy4 * ( zcv*zcv - 1._wp ) * ( zcv*zcv - 4._wp ) * ( ztv4(ji,jj+1,jl) + ztv4(ji,jj,jl) & 909 934 & - SIGN( 1._wp, zcv ) * ( ztv4(ji,jj+1,jl) - ztv4(ji,jj,jl) ) ) ) 910 935 END DO … … 941 966 942 967 943 SUBROUTINE nonosc _2d( pamsk, pdt, pu, puc, pv, pvc, ptc, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho )968 SUBROUTINE nonosc( pamsk, pdt, pu, pv, pt, pt_ups, pfu_ups, pfv_ups, pfu_ho, pfv_ho ) 944 969 !!--------------------------------------------------------------------- 945 970 !! *** ROUTINE nonosc *** 946 971 !! 947 972 !! ** Purpose : compute monotonic tracer fluxes from the upstream 948 973 !! scheme and the before field by a nonoscillatory algorithm 949 974 !! 950 !! ** Method : ... ??? 951 !! warning : pt and pt_ups must be masked, but the boundaries 952 !! conditions on the fluxes are not necessary zalezak (1979) 953 !! drange (1995) multi-dimensional forward-in-time and upstream- 954 !! in-space based differencing for fluid 975 !! ** Method : ... 955 976 !!---------------------------------------------------------------------- 956 977 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 957 978 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 958 979 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 959 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: puc ! ice i-velocity *A => u*e2*a960 980 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice j-velocity => v*e1 961 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pvc ! ice j-velocity *A => v*e1*a 962 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: ptc, pt, pt_ups ! before field & upstream guess of after field 963 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ups, pfu_ups ! upstream flux 981 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt, pt_ups ! before field & upstream guess of after field 982 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups, pfu_ups ! upstream flux 964 983 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho, pfu_ho ! monotonic flux 965 984 ! … … 1004 1023 CALL lbc_lnk_multi( 'icedyn_adv_umx', zti_ups, 'T', 1., ztj_ups, 'T', 1. ) 1005 1024 1006 !! this does not work ??1007 !! DO jj = 2, jpjm11008 !! DO ji = fs_2, fs_jpim11009 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji+1,jj ) - pt_ups (ji ,jj) ) .AND. &1010 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj+1) - pt_ups (ji ,jj) ) &1011 !! & ) THEN1012 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., zti_ups(ji+1,jj ) - zti_ups(ji ,jj) ) .AND. &1013 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., ztj_ups(ji,jj+1 ) - ztj_ups(ji ,jj) ) &1014 !! & ) THEN1015 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.1016 !! ENDIF1017 !! IF( SIGN( 1., pfu_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj) - pt_ups (ji-1,jj ) ) .AND. &1018 !! & SIGN( 1., pfv_ho(ji,jj) ) /= SIGN( 1., pt_ups (ji ,jj) - pt_ups (ji ,jj-1) ) &1019 !! & ) THEN1020 !! pfu_ho(ji,jj) = 0. ; pfv_ho(ji,jj) = 0.1021 !! ENDIF1022 !! ENDIF1023 !! END DO1024 !! END DO1025 1026 1025 DO jl = 1, jpl 1027 1026 DO jj = 2, jpjm1 … … 1030 1029 & pfv_ho(ji,jj,jl) * ( pt_ups(ji,jj+1,jl) - pt_ups(ji,jj,jl) ) <= 0. ) THEN 1031 1030 ! 1032 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0 .AND. &1033 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0 ) THEN1031 IF( pfu_ho(ji,jj,jl) * ( zti_ups(ji+1,jj,jl) - zti_ups(ji,jj,jl) ) <= 0. .AND. & 1032 & pfv_ho(ji,jj,jl) * ( ztj_ups(ji,jj+1,jl) - ztj_ups(ji,jj,jl) ) <= 0. ) THEN 1034 1033 pfu_ho(ji,jj,jl)=0. 1035 1034 pfv_ho(ji,jj,jl)=0. 1036 1035 ENDIF 1037 1036 ! 1038 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0 .AND. &1039 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0 ) THEN1037 IF( pfu_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji-1,jj,jl) ) <= 0. .AND. & 1038 & pfv_ho(ji,jj,jl) * ( pt_ups(ji ,jj,jl) - pt_ups(ji,jj-1,jl) ) <= 0. ) THEN 1040 1039 pfu_ho(ji,jj,jl)=0. 1041 1040 pfv_ho(ji,jj,jl)=0. … … 1049 1048 1050 1049 ENDIF 1051 1052 1050 1053 1051 ! Search local extrema … … 1086 1084 & + MAX( 0., pfv_ho(ji,jj ,jl) ) - MIN( 0., pfv_ho(ji,jj-1,jl) ) 1087 1085 ! 1088 zpos = zpos - ( pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1)) &1086 zpos = zpos - (pt(ji,jj,jl) * MIN( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MIN( 0., pv(ji,jj) - pv(ji,jj-1)) & 1089 1087 & ) * ( 1. - pamsk ) 1090 zneg = zneg + ( pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1)) &1088 zneg = zneg + (pt(ji,jj,jl) * MAX( 0., pu(ji,jj) - pu(ji-1,jj) ) + pt(ji,jj,jl) * MAX( 0., pv(ji,jj) - pv(ji,jj-1)) & 1091 1089 & ) * ( 1. - pamsk ) 1092 1090 ! … … 1154 1152 1155 1153 END DO 1156 1157 ! 1158 ! 1159 END SUBROUTINE nonosc_2d 1160 1161 SUBROUTINE limiter_x( pdt, pu, pt, pfu_ho, pfu_ups ) 1154 ! 1155 END SUBROUTINE nonosc 1156 1157 1158 SUBROUTINE limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 1162 1159 !!--------------------------------------------------------------------- 1163 1160 !! *** ROUTINE limiter_x *** … … 1165 1162 !! ** Purpose : compute flux limiter 1166 1163 !!---------------------------------------------------------------------- 1167 REAL(wp) , INTENT(in ):: pdt ! tracer time-step1168 REAL(wp), DIMENSION (:,: ), INTENT(in ):: pu ! ice i-velocity => u*e21169 REAL(wp), DIMENSION (:,:,:), INTENT(in ):: pt ! ice tracer1170 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfu_ho ! high orderflux1171 REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfu_ups ! upstreamflux1164 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1165 REAL(wp), DIMENSION(:,: ), INTENT(in ) :: pu ! ice i-velocity => u*e2 1166 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pt ! ice tracer 1167 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pfu_ups ! upstream flux 1168 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pfu_ho ! high order flux 1172 1169 ! 1173 1170 REAL(wp) :: Cr, Rjm, Rj, Rjp, uCFL, zpsi, zh3, zlimiter, Rr … … 1194 1191 Rjp = zslpx(ji+1,jj,jl) 1195 1192 1196 IF( PRESENT(pfu_ups)) THEN1193 IF( kn_limiter == 3 ) THEN 1197 1194 1198 1195 IF( pu(ji,jj) > 0. ) THEN ; Rr = Rjm … … 1210 1207 pfu_ho(ji,jj,jl) = pfu_ups(ji,jj,jl) + zlimiter 1211 1208 1212 ELSE 1209 ELSEIF( kn_limiter == 2 ) THEN 1213 1210 IF( Rj /= 0. ) THEN 1214 1211 IF( pu(ji,jj) > 0. ) THEN ; Cr = Rjm / Rj … … 1223 1220 ! -- van albada 2 -- 1224 1221 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1225 1226 1222 ! -- sweby (with beta=1) -- 1227 1223 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) ) … … 1254 1250 END SUBROUTINE limiter_x 1255 1251 1256 SUBROUTINE limiter_y( pdt, pv, pt, pfv_ho, pfv_ups ) 1252 1253 SUBROUTINE limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 1257 1254 !!--------------------------------------------------------------------- 1258 1255 !! *** ROUTINE limiter_y *** … … 1260 1257 !! ** Purpose : compute flux limiter 1261 1258 !!---------------------------------------------------------------------- 1262 REAL(wp) , INTENT(in ) 1263 REAL(wp), DIMENSION (:,: ), INTENT(in ) 1264 REAL(wp), DIMENSION (:,:,:), INTENT(in ) 1265 REAL(wp), DIMENSION (:,:,:), INTENT(in out) :: pfv_ho ! high orderflux1266 REAL(wp), DIMENSION (:,:,:), INTENT(in ), OPTIONAL :: pfv_ups ! upstreamflux1259 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 1260 REAL(wp), DIMENSION (:,: ), INTENT(in ) :: pv ! ice i-velocity => u*e2 1261 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pt ! ice tracer 1262 REAL(wp), DIMENSION (:,:,:), INTENT(in ) :: pfv_ups ! upstream flux 1263 REAL(wp), DIMENSION (:,:,:), INTENT(inout) :: pfv_ho ! high order flux 1267 1264 ! 1268 1265 REAL(wp) :: Cr, Rjm, Rj, Rjp, vCFL, zpsi, zh3, zlimiter, Rr … … 1289 1286 Rjp = zslpy(ji,jj+1,jl) 1290 1287 1291 IF( PRESENT(pfv_ups)) THEN1288 IF( kn_limiter == 3 ) THEN 1292 1289 1293 1290 IF( pv(ji,jj) > 0. ) THEN ; Rr = Rjm … … 1305 1302 pfv_ho(ji,jj,jl) = pfv_ups(ji,jj,jl) + zlimiter 1306 1303 1307 ELSE 1304 ELSEIF( kn_limiter == 2 ) THEN 1308 1305 1309 1306 IF( Rj /= 0. ) THEN … … 1319 1316 ! -- van albada 2 -- 1320 1317 !!zpsi = 2.*Cr / (Cr*Cr+1.) 1321 1322 1318 ! -- sweby (with beta=1) -- 1323 1319 !!zpsi = MAX( 0., MAX( MIN(1.,1.*Cr), MIN(1.,Cr) ) )
Note: See TracChangeset
for help on using the changeset viewer.
