Changeset 13469 for NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN
- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynspg_ts.F90
r13466 r13469 253 253 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 254 254 CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 255 DO jj = 2, jpjm1 256 DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters 257 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 258 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 259 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 260 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 261 END DO 262 END DO 255 DO_2D_00_00 256 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 257 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 258 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 259 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 260 END_2D 263 261 ELSE ! now suface pressure gradient 264 DO jj = 2, jpjm1 265 DO ji = fs_2, fs_jpim1 ! vector opt. 266 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 267 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 268 END DO 269 END DO 270 ENDIF 271 ! 272 ENDIF 273 ! 274 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 275 DO ji = fs_2, fs_jpim1 276 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 277 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 278 END DO 279 END DO 262 DO_2D_00_00 263 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 264 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 265 END_2D 266 ENDIF 267 ! 268 ENDIF 269 ! 270 DO_2D_00_00 271 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 272 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 273 END_2D 280 274 ! 281 275 ! != Add bottom stress contribution from baroclinic velocities =! … … 287 281 IF( ln_apr_dyn ) THEN 288 282 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 292 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 293 END DO 294 END DO 283 DO_2D_00_00 284 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 285 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 286 END_2D 295 287 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 296 288 zztmp = grav * r1_2 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 ! vector opt. 299 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 300 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 301 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 302 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 303 END DO 304 END DO 289 DO_2D_00_00 290 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 291 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 292 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 293 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 294 END_2D 305 295 ENDIF 306 296 ENDIF … … 309 299 ! ! ---------------------------------- ! 310 300 IF( ln_bt_fw ) THEN ! Add wind forcing 311 DO jj = 2, jpjm1 312 DO ji = fs_2, fs_jpim1 ! vector opt. 313 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj) 314 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj) 315 END DO 316 END DO 301 DO_2D_00_00 302 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj) 303 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj) 304 END_2D 317 305 ELSE 318 306 zztmp = r1_rau0 * r1_2 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 ! vector opt. 321 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 322 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 323 END DO 324 END DO 307 DO_2D_00_00 308 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj) 309 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj) 310 END_2D 325 311 ENDIF 326 312 ! … … 457 443 ! 458 444 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 459 DO jj = 1, jpj 460 DO ji = 1, jpim1 ! not jpi-column 461 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 462 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 463 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 464 END DO 465 END DO 466 DO jj = 1, jpjm1 ! not jpj-row 467 DO ji = 1, jpi 468 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 469 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 470 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 471 END DO 472 END DO 445 DO_2D_11_10 446 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 447 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 448 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 449 END_2D 450 DO_2D_10_11 451 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 452 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 453 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 454 END_2D 473 455 ! 474 456 ENDIF … … 526 508 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 527 509 !-------------------------------------------------------------------------! 528 DO jj = 2, jpjm1 ! INNER domain 529 DO ji = 2, jpim1 530 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 531 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 532 END DO 533 END DO 510 DO_2D_00_00 511 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 512 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rdtbt * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 513 END_2D 534 514 ! 535 515 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) … … 553 533 ! Sea Surface Height at u-,v-points (vvl case only) 554 534 IF( .NOT.ln_linssh ) THEN 555 DO jj = 2, jpjm1 ! INNER domain, will be extended to whole domain later 556 DO ji = 2, jpim1 ! NO Vector Opt. 557 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 558 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 559 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 560 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 561 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 562 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 563 END DO 564 END DO 535 DO_2D_00_00 536 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 537 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 538 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 539 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 540 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 541 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 542 END_2D 565 543 ENDIF 566 544 ! … … 575 553 ! ! Surface pressure gradient 576 554 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 577 DO jj = 2, jpjm1 578 DO ji = 2, jpim1 579 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 580 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 581 END DO 582 END DO 555 DO_2D_00_00 556 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 557 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 558 END_2D 583 559 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 584 560 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters … … 595 571 ! Add tidal astronomical forcing if defined 596 572 IF ( ln_tide .AND. ln_tide_pot ) THEN 597 DO jj = 2, jpjm1 598 DO ji = fs_2, fs_jpim1 ! vector opt. 599 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 600 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 601 END DO 602 END DO 573 DO_2D_00_00 574 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 575 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 576 END_2D 603 577 ENDIF 604 578 ! … … 606 580 !jth do implicitly instead 607 581 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 608 DO jj = 2, jpjm1 609 DO ji = fs_2, fs_jpim1 ! vector opt. 610 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 611 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 612 END DO 613 END DO 582 DO_2D_00_00 583 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 584 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 585 END_2D 614 586 ENDIF 615 587 ! … … 626 598 !------------------------------------------------------------------------------------------------------------------------! 627 599 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 628 DO jj = 2, jpjm1 629 DO ji = fs_2, fs_jpim1 ! vector opt. 630 ua_e(ji,jj) = ( un_e(ji,jj) & 631 & + rdtbt * ( zu_spg(ji,jj) & 632 & + zu_trd(ji,jj) & 633 & + zu_frc(ji,jj) ) & 634 & ) * ssumask(ji,jj) 635 636 va_e(ji,jj) = ( vn_e(ji,jj) & 637 & + rdtbt * ( zv_spg(ji,jj) & 638 & + zv_trd(ji,jj) & 639 & + zv_frc(ji,jj) ) & 640 & ) * ssvmask(ji,jj) 641 END DO 642 END DO 600 DO_2D_00_00 601 ua_e(ji,jj) = ( un_e(ji,jj) & 602 & + rdtbt * ( zu_spg(ji,jj) & 603 & + zu_trd(ji,jj) & 604 & + zu_frc(ji,jj) ) & 605 & ) * ssumask(ji,jj) 606 607 va_e(ji,jj) = ( vn_e(ji,jj) & 608 & + rdtbt * ( zv_spg(ji,jj) & 609 & + zv_trd(ji,jj) & 610 & + zv_frc(ji,jj) ) & 611 & ) * ssvmask(ji,jj) 612 END_2D 643 613 ! 644 614 ELSE !* Flux form 645 DO jj = 2, jpjm1 646 DO ji = 2, jpim1 647 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 648 ! ! backward interpolated depth used in spg terms at jn+1/2 649 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 650 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 651 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 652 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 653 ! ! inverse depth at jn+1 654 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 655 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 656 ! 657 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 658 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 659 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 660 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 661 ! 662 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 663 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 664 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 665 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv 666 END DO 667 END DO 615 DO_2D_00_00 616 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 617 ! ! backward interpolated depth used in spg terms at jn+1/2 618 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 619 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 620 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 621 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 622 ! ! inverse depth at jn+1 623 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 624 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 625 ! 626 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 627 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 628 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 629 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 630 ! 631 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 632 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 633 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 634 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hv 635 END_2D 668 636 ENDIF 669 637 !jth implicit bottom friction: 670 638 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 671 DO jj = 2, jpjm1 672 DO ji = fs_2, fs_jpim1 ! vector opt. 673 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj)) 674 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj)) 675 END DO 676 END DO 639 DO_2D_00_00 640 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj)) 641 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj)) 642 END_2D 677 643 ENDIF 678 644 … … 737 703 IF (ln_bt_fw) THEN 738 704 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 739 DO jj = 1, jpj 740 DO ji = 1, jpi 741 zun_save = un_adv(ji,jj) 742 zvn_save = vn_adv(ji,jj) 743 ! ! apply the previously computed correction 744 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 745 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 746 ! ! Update corrective fluxes for next time step 747 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 748 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 749 ! ! Save integrated transport for next computation 750 ub2_b(ji,jj) = zun_save 751 vb2_b(ji,jj) = zvn_save 752 END DO 753 END DO 705 DO_2D_11_11 706 zun_save = un_adv(ji,jj) 707 zvn_save = vn_adv(ji,jj) 708 ! ! apply the previously computed correction 709 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 710 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 711 ! ! Update corrective fluxes for next time step 712 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 713 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 714 ! ! Save integrated transport for next computation 715 ub2_b(ji,jj) = zun_save 716 vb2_b(ji,jj) = zvn_save 717 END_2D 754 718 ELSE 755 719 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero … … 770 734 ELSE 771 735 ! At this stage, ssha has been corrected: compute new depths at velocity points 772 DO jj = 1, jpjm1 773 DO ji = 1, jpim1 ! NO Vector Opt. 774 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 775 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 776 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 777 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 778 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 779 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 780 END DO 781 END DO 736 DO_2D_10_10 737 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 738 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 739 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 740 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 741 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 742 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 743 END_2D 782 744 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 783 745 ! … … 794 756 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 795 757 DO jk = 1, jpkm1 796 u n(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)797 v n(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)758 uu(:,:,jk,Nii) = ( uu(:,:,jk,Nii) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk) 759 vv(:,:,jk,Nii) = ( vv(:,:,jk,Nii) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 798 760 END DO 799 761 … … 802 764 CALL lbc_lnk_multi( 'dynspg_ts', zuwdav2, 'U', 1._wp, zvwdav2, 'V', 1._wp) 803 765 DO jk = 1, jpkm1 804 u n(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &805 & + zuwdav2(:,:)*(u n(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk)806 v n(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &807 & + zvwdav2(:,:)*(v n(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk)766 uu(:,:,jk,Nii) = ( un_adv(:,:)*r1_hu_n(:,:) & 767 & + zuwdav2(:,:)*(uu(:,:,jk,Nii) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk) 768 vv(:,:,jk,Nii) = ( vn_adv(:,:)*r1_hv_n(:,:) & 769 & + zvwdav2(:,:)*(vv(:,:,jk,Nii) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk) 808 770 END DO 809 771 END IF … … 1007 969 ! Max courant number for ext. grav. waves 1008 970 ! 1009 DO jj = 1, jpj 1010 DO ji =1, jpi 1011 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1012 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1013 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 1014 END DO 1015 END DO 971 DO_2D_11_11 972 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 973 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 974 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 975 END_2D 1016 976 ! 1017 977 zcmax = MAXVAL( zcu(:,:) ) … … 1133 1093 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1134 1094 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1135 DO jj = 1, jpjm1 1136 DO ji = 1, jpim1 1137 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1138 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1139 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1140 END DO 1141 END DO 1095 DO_2D_10_10 1096 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1097 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1098 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1099 END_2D 1142 1100 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1143 DO jj = 1, jpjm1 1144 DO ji = 1, jpim1 1145 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1146 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1147 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1148 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1149 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1150 END DO 1151 END DO 1101 DO_2D_10_10 1102 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1103 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1104 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1105 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1106 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1107 END_2D 1152 1108 END SELECT 1153 1109 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1154 1110 ! 1155 1111 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1156 DO jj = 2, jpj 1157 DO ji = 2, jpi 1158 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1159 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1160 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1161 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1162 END DO 1163 END DO 1112 DO_2D_01_01 1113 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1114 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1115 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1116 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1117 END_2D 1164 1118 ! 1165 1119 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1166 1120 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1167 DO jj = 2, jpj 1168 DO ji = 2, jpi 1169 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1170 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1171 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1172 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1173 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1174 END DO 1175 END DO 1121 DO_2D_01_01 1122 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1123 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1124 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1125 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1126 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1127 END_2D 1176 1128 ! 1177 1129 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! … … 1200 1152 ! 1201 1153 !zhf(:,:) = hbatf(:,:) 1202 DO jj = 1, jpjm1 1203 DO ji = 1, jpim1 1204 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1205 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1206 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1207 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1208 END DO 1209 END DO 1154 DO_2D_10_10 1155 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1156 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1157 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1158 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1159 END_2D 1210 1160 ENDIF 1211 1161 ! … … 1221 1171 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1222 1172 ! JC: TBC. hf should be greater than 0 1223 DO jj = 1, jpj 1224 DO ji = 1, jpi 1225 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1226 END DO 1227 END DO 1173 DO_2D_11_11 1174 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1175 END_2D 1228 1176 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1229 1177 END SELECT … … 1246 1194 SELECT CASE( nvor_scheme ) 1247 1195 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1248 DO jj = 2, jpjm1 1249 DO ji = 2, jpim1 1250 z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 1251 z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1252 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1253 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 1254 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1255 ! 1256 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1257 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 1258 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1259 END DO 1260 END DO 1196 DO_2D_00_00 1197 z1_hu = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 1198 z1_hv = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1199 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1200 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 1201 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1202 ! 1203 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1204 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 1205 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1206 END_2D 1261 1207 ! 1262 1208 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1263 DO jj = 2, jpjm1 1264 DO ji = fs_2, fs_jpim1 ! vector opt. 1265 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1266 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1267 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1268 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1269 ! energy conserving formulation for planetary vorticity term 1270 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1271 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1272 END DO 1273 END DO 1209 DO_2D_00_00 1210 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1211 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1212 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1213 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1214 ! energy conserving formulation for planetary vorticity term 1215 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1216 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1217 END_2D 1274 1218 ! 1275 1219 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1276 DO jj = 2, jpjm1 1277 DO ji = fs_2, fs_jpim1 ! vector opt. 1278 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1279 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1280 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1281 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1282 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1283 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1284 END DO 1285 END DO 1220 DO_2D_00_00 1221 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1222 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1223 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1224 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1225 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1226 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1227 END_2D 1286 1228 ! 1287 1229 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1288 DO jj = 2, jpjm1 1289 DO ji = fs_2, fs_jpim1 ! vector opt. 1290 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1291 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1292 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1293 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1294 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1295 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1296 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1297 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1298 END DO 1299 END DO 1230 DO_2D_00_00 1231 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1232 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1233 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1234 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1235 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1236 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1237 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1238 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1239 END_2D 1300 1240 ! 1301 1241 END SELECT … … 1322 1262 ! 1323 1263 IF( ln_wd_dl_rmp ) THEN 1324 DO jj = 1, jpj 1325 DO ji = 1, jpi 1326 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1327 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1328 ptmsk(ji,jj) = 1._wp 1329 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1330 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1331 ELSE 1332 ptmsk(ji,jj) = 0._wp 1333 ENDIF 1334 END DO 1335 END DO 1264 DO_2D_11_11 1265 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1266 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1267 ptmsk(ji,jj) = 1._wp 1268 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1269 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1270 ELSE 1271 ptmsk(ji,jj) = 0._wp 1272 ENDIF 1273 END_2D 1336 1274 ELSE 1337 DO jj = 1, jpj 1338 DO ji = 1, jpi 1339 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1340 ELSE ; ptmsk(ji,jj) = 0._wp 1341 ENDIF 1342 END DO 1343 END DO 1275 DO_2D_11_11 1276 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1277 ELSE ; ptmsk(ji,jj) = 0._wp 1278 ENDIF 1279 END_2D 1344 1280 ENDIF 1345 1281 ! … … 1365 1301 !!---------------------------------------------------------------------- 1366 1302 ! 1367 DO jj = 1, jpj 1368 DO ji = 1, jpim1 ! not jpi-column 1369 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1370 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1371 ENDIF 1372 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1373 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1374 END DO 1375 END DO 1376 ! 1377 DO jj = 1, jpjm1 ! not jpj-row 1378 DO ji = 1, jpi 1379 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1380 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1381 ENDIF 1382 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1383 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1384 END DO 1385 END DO 1303 DO_2D_11_10 1304 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1305 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1306 ENDIF 1307 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1308 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1309 END_2D 1310 ! 1311 DO_2D_10_11 1312 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1313 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1314 ENDIF 1315 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1316 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1317 END_2D 1386 1318 ! 1387 1319 END SUBROUTINE wad_Umsk … … 1399 1331 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1400 1332 !!---------------------------------------------------------------------- 1401 DO jj = 2, jpjm1 1402 DO ji = 2, jpim1 1403 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1404 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1405 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1406 & > rn_wdmin1 + rn_wdmin2 1407 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1408 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1409 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1410 IF(ll_tmp1) THEN 1411 zcpx(ji,jj) = 1.0_wp 1412 ELSEIF(ll_tmp2) THEN 1413 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1414 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1415 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1416 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1417 ELSE 1418 zcpx(ji,jj) = 0._wp 1419 ENDIF 1420 ! 1421 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1422 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1423 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1424 & > rn_wdmin1 + rn_wdmin2 1425 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1426 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1427 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1428 1429 IF(ll_tmp1) THEN 1430 zcpy(ji,jj) = 1.0_wp 1431 ELSE IF(ll_tmp2) THEN 1432 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1433 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1434 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1435 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1436 ELSE 1437 zcpy(ji,jj) = 0._wp 1438 ENDIF 1439 END DO 1440 END DO 1333 DO_2D_00_00 1334 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1335 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1336 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1337 & > rn_wdmin1 + rn_wdmin2 1338 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1339 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1340 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1341 IF(ll_tmp1) THEN 1342 zcpx(ji,jj) = 1.0_wp 1343 ELSEIF(ll_tmp2) THEN 1344 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1345 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1346 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1347 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1348 ELSE 1349 zcpx(ji,jj) = 0._wp 1350 ENDIF 1351 ! 1352 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1353 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1354 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1355 & > rn_wdmin1 + rn_wdmin2 1356 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1357 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1358 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1359 1360 IF(ll_tmp1) THEN 1361 zcpy(ji,jj) = 1.0_wp 1362 ELSE IF(ll_tmp2) THEN 1363 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1364 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1365 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1366 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1367 ELSE 1368 zcpy(ji,jj) = 0._wp 1369 ENDIF 1370 END_2D 1441 1371 1442 1372 END SUBROUTINE wad_spg … … 1467 1397 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1468 1398 1469 DO jj = 2, jpjm1 1470 DO ji = 2, jpim1 ! INNER domain 1471 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1472 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1473 END DO 1474 END DO 1399 DO_2D_00_00 1400 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1401 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1402 END_2D 1475 1403 ELSE ! bottom friction only 1476 DO jj = 2, jpjm1 1477 DO ji = 2, jpim1 ! INNER domain 1478 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1479 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1480 END DO 1481 END DO 1404 DO_2D_00_00 1405 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1406 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1407 END_2D 1482 1408 ENDIF 1483 1409 ! … … 1486 1412 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1487 1413 1488 DO jj = 2, jpjm1 1489 DO ji = 2, jpim1 ! INNER domain 1490 ikbu = mbku(ji,jj) 1491 ikbv = mbkv(ji,jj) 1492 zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 1493 zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 1494 END DO 1495 END DO 1414 DO_2D_00_00 1415 ikbu = mbku(ji,jj) 1416 ikbv = mbkv(ji,jj) 1417 zu_i(ji,jj) = uu(ji,jj,ikbu,Nii) - un_b(ji,jj) 1418 zv_i(ji,jj) = vv(ji,jj,ikbv,Nii) - vn_b(ji,jj) 1419 END_2D 1496 1420 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1497 1421 1498 DO jj = 2, jpjm1 1499 DO ji = 2, jpim1 ! INNER domain 1500 ikbu = mbku(ji,jj) 1501 ikbv = mbkv(ji,jj) 1502 zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 1503 zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 1504 END DO 1505 END DO 1422 DO_2D_00_00 1423 ikbu = mbku(ji,jj) 1424 ikbv = mbkv(ji,jj) 1425 zu_i(ji,jj) = uu(ji,jj,ikbu,Nnn) - ub_b(ji,jj) 1426 zv_i(ji,jj) = vv(ji,jj,ikbv,Nnn) - vb_b(ji,jj) 1427 END_2D 1506 1428 ENDIF 1507 1429 ! 1508 1430 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1509 1431 zztmp = -1._wp / rdtbt 1510 DO jj = 2, jpjm1 1511 DO ji = 2, jpim1 ! INNER domain 1512 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1513 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1514 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1515 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1516 END DO 1517 END DO 1432 DO_2D_00_00 1433 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1434 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1435 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1436 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1437 END_2D 1518 1438 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1519 1439 1520 DO jj = 2, jpjm1 1521 DO ji = 2, jpim1 ! INNER domain 1522 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1523 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1524 END DO 1525 END DO 1440 DO_2D_00_00 1441 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1442 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1443 END_2D 1526 1444 END IF 1527 1445 ! … … 1532 1450 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1533 1451 1534 DO jj = 2, jpjm1 1535 DO ji = 2, jpim1 ! INNER domain 1536 iktu = miku(ji,jj) 1537 iktv = mikv(ji,jj) 1538 zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 1539 zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 1540 END DO 1541 END DO 1452 DO_2D_00_00 1453 iktu = miku(ji,jj) 1454 iktv = mikv(ji,jj) 1455 zu_i(ji,jj) = uu(ji,jj,iktu,Nii) - un_b(ji,jj) 1456 zv_i(ji,jj) = vv(ji,jj,iktv,Nii) - vn_b(ji,jj) 1457 END_2D 1542 1458 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1543 1459 1544 DO jj = 2, jpjm1 1545 DO ji = 2, jpim1 ! INNER domain 1546 iktu = miku(ji,jj) 1547 iktv = mikv(ji,jj) 1548 zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 1549 zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 1550 END DO 1551 END DO 1460 DO_2D_00_00 1461 iktu = miku(ji,jj) 1462 iktv = mikv(ji,jj) 1463 zu_i(ji,jj) = uu(ji,jj,iktu,Nnn) - ub_b(ji,jj) 1464 zv_i(ji,jj) = vv(ji,jj,iktv,Nnn) - vb_b(ji,jj) 1465 END_2D 1552 1466 ENDIF 1553 1467 ! 1554 1468 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1555 1469 1556 DO jj = 2, jpjm1 1557 DO ji = 2, jpim1 ! INNER domain 1558 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1559 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1560 END DO 1561 END DO 1470 DO_2D_00_00 1471 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1472 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1473 END_2D 1562 1474 ! 1563 1475 ENDIF -
NEMO/branches/2020/temporary_r4_trunk/src/OCE/DYN/dynzdf.F90
r13466 r13469 110 110 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 111 111 DO jk = 1, jpkm1 112 ua(:,:,jk) = ( u b(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)113 va(:,:,jk) = ( v b(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)112 ua(:,:,jk) = ( uu(:,:,jk,Nnn) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 113 va(:,:,jk) = ( vv(:,:,jk,Nnn) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 114 114 END DO 115 115 ELSE ! applied on thickness weighted velocity 116 116 DO jk = 1, jpkm1 117 ua(:,:,jk) = ( e3u_b(:,:,jk) * u b(:,:,jk) &117 ua(:,:,jk) = ( e3u_b(:,:,jk) * uu(:,:,jk,Nnn) & 118 118 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 119 va(:,:,jk) = ( e3v_b(:,:,jk) * v b(:,:,jk) &119 va(:,:,jk) = ( e3v_b(:,:,jk) * vv(:,:,jk,Nnn) & 120 120 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 121 121 END DO … … 131 131 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 132 132 END DO 133 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 134 DO ji = fs_2, fs_jpim1 ! vector opt. 135 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 136 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 133 DO_2D_00_00 134 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 135 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 136 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 137 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 138 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 139 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 140 END_2D 141 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF) 142 DO_2D_00_00 143 iku = miku(ji,jj) ! top ocean level at u- and v-points 144 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 137 145 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 138 146 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 139 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 140 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 141 END DO 142 END DO 143 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF) 144 DO jj = 2, jpjm1 145 DO ji = fs_2, fs_jpim1 ! vector opt. 146 iku = miku(ji,jj) ! top ocean level at u- and v-points 147 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 148 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 149 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 150 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 151 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 152 END DO 153 END DO 147 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 148 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 149 END_2D 154 150 END IF 155 151 ENDIF … … 162 158 SELECT CASE( nldf_dyn ) 163 159 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 164 DO jk = 1, jpkm1 165 DO jj = 2, jpjm1 166 DO ji = fs_2, fs_jpim1 ! vector opt. 167 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 168 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 169 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 170 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 171 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 172 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 173 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 174 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 175 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 176 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 177 END DO 178 END DO 179 END DO 160 DO_3D_00_00( 1, jpkm1 ) 161 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 162 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 163 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 164 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 165 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 166 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 167 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 168 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 169 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 170 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 171 END_3D 180 172 CASE DEFAULT ! iso-level lateral mixing 181 DO jk = 1, jpkm1 182 DO jj = 2, jpjm1 183 DO ji = fs_2, fs_jpim1 ! vector opt. 184 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 185 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 186 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 187 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 188 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 189 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 190 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 191 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 192 END DO 193 END DO 194 END DO 173 DO_3D_00_00( 1, jpkm1 ) 174 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 175 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 176 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 177 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 178 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua 179 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp ) 180 zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp ) 181 zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) ) 182 END_3D 195 183 END SELECT 196 DO jj = 2, jpjm1 !* Surface boundary conditions 197 DO ji = fs_2, fs_jpim1 ! vector opt. 198 zwi(ji,jj,1) = 0._wp 199 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 200 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 201 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 202 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 203 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 204 END DO 205 END DO 184 DO_2D_00_00 185 zwi(ji,jj,1) = 0._wp 186 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 187 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw_n(ji,jj,2) ) * wumask(ji,jj,2) 188 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 189 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) 190 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 0._wp ) ) 191 END_2D 206 192 ELSE 207 193 SELECT CASE( nldf_dyn ) 208 194 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 209 DO jk = 1, jpkm1 210 DO jj = 2, jpjm1 211 DO ji = fs_2, fs_jpim1 ! vector opt. 212 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 213 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 214 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 215 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 216 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 217 zwi(ji,jj,jk) = zzwi 218 zws(ji,jj,jk) = zzws 219 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 220 END DO 221 END DO 222 END DO 195 DO_3D_00_00( 1, jpkm1 ) 196 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 197 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 198 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 199 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 200 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 201 zwi(ji,jj,jk) = zzwi 202 zws(ji,jj,jk) = zzws 203 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 204 END_3D 223 205 CASE DEFAULT ! iso-level lateral mixing 224 DO jk = 1, jpkm1 225 DO jj = 2, jpjm1 226 DO ji = fs_2, fs_jpim1 ! vector opt. 227 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 228 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 229 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 230 zwi(ji,jj,jk) = zzwi 231 zws(ji,jj,jk) = zzws 232 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 233 END DO 234 END DO 235 END DO 206 DO_3D_00_00( 1, jpkm1 ) 207 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 208 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 209 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 210 zwi(ji,jj,jk) = zzwi 211 zws(ji,jj,jk) = zzws 212 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 213 END_3D 236 214 END SELECT 237 DO jj = 2, jpjm1 !* Surface boundary conditions 238 DO ji = fs_2, fs_jpim1 ! vector opt. 239 zwi(ji,jj,1) = 0._wp 240 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 241 END DO 242 END DO 215 DO_2D_00_00 216 zwi(ji,jj,1) = 0._wp 217 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 218 END_2D 243 219 ENDIF 244 220 ! … … 251 227 ! 252 228 IF ( ln_drgimp ) THEN ! implicit bottom friction 253 DO jj = 2, jpjm1 254 DO ji = 2, jpim1 255 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 229 DO_2D_00_00 230 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 231 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 232 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 233 END_2D 234 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit) 235 DO_2D_00_00 236 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 237 iku = miku(ji,jj) ! ocean top level at u- and v-points 256 238 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 257 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 258 END DO 259 END DO 260 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit) 261 DO jj = 2, jpjm1 262 DO ji = 2, jpim1 263 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 264 iku = miku(ji,jj) ! ocean top level at u- and v-points 265 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 266 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 267 END DO 268 END DO 239 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 240 END_2D 269 241 END IF 270 242 ENDIF … … 285 257 !----------------------------------------------------------------------- 286 258 ! 287 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 288 DO jj = 2, jpjm1 289 DO ji = fs_2, fs_jpim1 ! vector opt. 290 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 291 END DO 292 END DO 293 END DO 294 ! 295 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 296 DO ji = fs_2, fs_jpim1 ! vector opt. 297 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 298 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 299 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 300 END DO 301 END DO 302 DO jk = 2, jpkm1 303 DO jj = 2, jpjm1 304 DO ji = fs_2, fs_jpim1 305 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 306 END DO 307 END DO 308 END DO 309 ! 310 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 311 DO ji = fs_2, fs_jpim1 ! vector opt. 312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 313 END DO 314 END DO 315 DO jk = jpk-2, 1, -1 316 DO jj = 2, jpjm1 317 DO ji = fs_2, fs_jpim1 318 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 319 END DO 320 END DO 321 END DO 259 DO_3D_00_00( 2, jpkm1 ) 260 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 261 END_3D 262 ! 263 DO_2D_00_00 264 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 265 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 266 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 267 END_2D 268 DO_3D_00_00( 2, jpkm1 ) 269 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 270 END_3D 271 ! 272 DO_2D_00_00 273 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 274 END_2D 275 DO_3D_00_00( jpk-2, 1, -1 ) 276 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 277 END_3D 322 278 ! 323 279 ! !== Vertical diffusion on v ==! … … 328 284 SELECT CASE( nldf_dyn ) 329 285 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 330 DO jk = 1, jpkm1 331 DO jj = 2, jpjm1 332 DO ji = fs_2, fs_jpim1 ! vector opt. 333 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 334 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 335 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 336 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 337 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 338 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 339 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 340 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 341 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 342 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 343 END DO 344 END DO 345 END DO 286 DO_3D_00_00( 1, jpkm1 ) 287 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 288 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 289 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 290 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 291 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 292 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 293 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 294 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 295 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 296 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 297 END_3D 346 298 CASE DEFAULT ! iso-level lateral mixing 347 DO jk = 1, jpkm1 348 DO jj = 2, jpjm1 349 DO ji = fs_2, fs_jpim1 ! vector opt. 350 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 351 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 352 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 353 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 354 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 355 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 356 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 357 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 358 END DO 359 END DO 360 END DO 299 DO_3D_00_00( 1, jpkm1 ) 300 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 301 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 302 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 303 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 304 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va 305 zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp ) 306 zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp ) 307 zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 308 END_3D 361 309 END SELECT 362 DO jj = 2, jpjm1 !* Surface boundary conditions 363 DO ji = fs_2, fs_jpim1 ! vector opt. 364 zwi(ji,jj,1) = 0._wp 365 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 366 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 367 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 368 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 369 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 370 END DO 371 END DO 310 DO_2D_00_00 311 zwi(ji,jj,1) = 0._wp 312 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 313 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw_n(ji,jj,2) ) * wvmask(ji,jj,2) 314 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 315 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) 316 zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 0._wp ) ) 317 END_2D 372 318 ELSE 373 319 SELECT CASE( nldf_dyn ) 374 320 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 375 DO jk = 1, jpkm1 376 DO jj = 2, jpjm1 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 379 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 380 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 381 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 382 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 383 zwi(ji,jj,jk) = zzwi 384 zws(ji,jj,jk) = zzws 385 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 386 END DO 387 END DO 388 END DO 321 DO_3D_00_00( 1, jpkm1 ) 322 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 323 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 324 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 325 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 326 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 327 zwi(ji,jj,jk) = zzwi 328 zws(ji,jj,jk) = zzws 329 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 330 END_3D 389 331 CASE DEFAULT ! iso-level lateral mixing 390 DO jk = 1, jpkm1 391 DO jj = 2, jpjm1 392 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 394 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 395 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 396 zwi(ji,jj,jk) = zzwi 397 zws(ji,jj,jk) = zzws 398 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 399 END DO 400 END DO 401 END DO 332 DO_3D_00_00( 1, jpkm1 ) 333 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 334 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 335 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 336 zwi(ji,jj,jk) = zzwi 337 zws(ji,jj,jk) = zzws 338 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 339 END_3D 402 340 END SELECT 403 DO jj = 2, jpjm1 !* Surface boundary conditions 404 DO ji = fs_2, fs_jpim1 ! vector opt. 405 zwi(ji,jj,1) = 0._wp 406 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 407 END DO 408 END DO 341 DO_2D_00_00 342 zwi(ji,jj,1) = 0._wp 343 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 344 END_2D 409 345 ENDIF 410 346 ! … … 416 352 ! 417 353 IF( ln_drgimp ) THEN 418 DO jj = 2, jpjm1 419 DO ji = 2, jpim1 420 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 354 DO_2D_00_00 355 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 356 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 357 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 358 END_2D 359 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 360 DO_2D_00_00 361 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 421 362 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 422 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 423 END DO 424 END DO 425 IF ( ln_isfcav.OR.ln_drgice_imp ) THEN 426 DO jj = 2, jpjm1 427 DO ji = 2, jpim1 428 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 429 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 430 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 431 END DO 432 END DO 363 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 364 END_2D 433 365 ENDIF 434 366 ENDIF … … 449 381 !----------------------------------------------------------------------- 450 382 ! 451 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 452 DO jj = 2, jpjm1 453 DO ji = fs_2, fs_jpim1 ! vector opt. 454 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 455 END DO 456 END DO 457 END DO 458 ! 459 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 460 DO ji = fs_2, fs_jpim1 ! vector opt. 461 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 462 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 463 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 464 END DO 465 END DO 466 DO jk = 2, jpkm1 467 DO jj = 2, jpjm1 468 DO ji = fs_2, fs_jpim1 ! vector opt. 469 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 470 END DO 471 END DO 472 END DO 473 ! 474 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 475 DO ji = fs_2, fs_jpim1 ! vector opt. 476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 477 END DO 478 END DO 479 DO jk = jpk-2, 1, -1 480 DO jj = 2, jpjm1 481 DO ji = fs_2, fs_jpim1 482 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 483 END DO 484 END DO 485 END DO 383 DO_3D_00_00( 2, jpkm1 ) 384 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 385 END_3D 386 ! 387 DO_2D_00_00 388 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 389 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 390 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 391 END_2D 392 DO_3D_00_00( 2, jpkm1 ) 393 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 394 END_3D 395 ! 396 DO_2D_00_00 397 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 398 END_2D 399 DO_3D_00_00( jpk-2, 1, -1 ) 400 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 401 END_3D 486 402 ! 487 403 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 488 ztrdu(:,:,:) = ( ua(:,:,:) - u b(:,:,:) ) / r2dt - ztrdu(:,:,:)489 ztrdv(:,:,:) = ( va(:,:,:) - v b(:,:,:) ) / r2dt - ztrdv(:,:,:)404 ztrdu(:,:,:) = ( ua(:,:,:) - uu(:,:,:,Nnn) ) / r2dt - ztrdu(:,:,:) 405 ztrdv(:,:,:) = ( va(:,:,:) - vv(:,:,:,Nnn) ) / r2dt - ztrdv(:,:,:) 490 406 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 491 407 DEALLOCATE( ztrdu, ztrdv )
Note: See TracChangeset
for help on using the changeset viewer.