- Timestamp:
- 2021-12-03T20:32:50+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r14318_RK3_stage1
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynvor.F90 (modified) (22 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r14318_RK3_stage1
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/dev_r14318_RK3_stage1/src/OCE/DYN/dynvor.F90
r14547 r15574 321 321 INTEGER :: ji, jj, jk ! dummy loop indices 322 322 REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars 323 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwt ! 2D workspace 324 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 325 !!---------------------------------------------------------------------- 326 ! 327 IF( kt == nit000 ) THEN 328 IF(lwp) WRITE(numout,*) 329 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 330 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 323 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwt ! 2D workspace 324 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined 325 !!---------------------------------------------------------------------- 326 ! 327 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 328 IF( kt == nit000 ) THEN 329 IF(lwp) WRITE(numout,*) 330 IF(lwp) WRITE(numout,*) 'dyn:vor_enT : vorticity term: t-point energy conserving scheme' 331 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 332 ENDIF 331 333 ENDIF 332 334 ! … … 335 337 ! 336 338 CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used 337 ALLOCATE( zwz( jpi,jpj,jpk) )339 ALLOCATE( zwz(A2D(nn_hls),jpk) ) 338 340 DO jk = 1, jpkm1 ! Horizontal slab 339 DO_2D( 1, 0, 1, 0)341 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 340 342 zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 341 343 & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj) 342 344 END_2D 343 345 IF( ln_dynvor_msk ) THEN ! mask relative vorticity 344 DO_2D( 1, 0, 1, 0)346 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 345 347 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 346 348 END_2D 347 349 ENDIF 348 350 END DO 349 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )351 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 350 352 ! 351 353 END SELECT … … 358 360 ! 359 361 CASE ( np_COR ) !* Coriolis (planetary vorticity) 360 zwt(:,:) = ff_t(:,:) * e1e2t(:,:)*e3t(:,:,jk,Kmm) 362 DO_2D( 0, 1, 0, 1 ) 363 zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm) 364 END_2D 361 365 CASE ( np_RVO ) !* relative vorticity 362 366 DO_2D( 0, 1, 0, 1 ) … … 437 441 INTEGER :: ji, jj, jk ! dummy loop indices 438 442 REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars 439 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz ! 2D workspace 440 !!---------------------------------------------------------------------- 441 ! 442 IF( kt == nit000 ) THEN 443 IF(lwp) WRITE(numout,*) 444 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 445 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 443 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace 444 !!---------------------------------------------------------------------- 445 ! 446 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 447 IF( kt == nit000 ) THEN 448 IF(lwp) WRITE(numout,*) 449 IF(lwp) WRITE(numout,*) 'dyn:vor_ene : vorticity term: energy conserving scheme' 450 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 451 ENDIF 446 452 ENDIF 447 453 ! … … 452 458 SELECT CASE( kvor ) !== vorticity considered ==! 453 459 CASE ( np_COR ) !* Coriolis (planetary vorticity) 454 zwz(:,:) = ff_f(:,:) 460 DO_2D( 1, 0, 1, 0 ) 461 zwz(ji,jj) = ff_f(ji,jj) 462 END_2D 455 463 CASE ( np_RVO ) !* relative vorticity 456 464 DO_2D( 1, 0, 1, 0 ) … … 518 526 #endif 519 527 ! !== horizontal fluxes ==! 520 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 521 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 528 DO_2D( 1, 1, 1, 1 ) 529 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 530 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 531 END_2D 522 532 ! 523 533 ! !== compute and add the vorticity term trend =! … … 564 574 INTEGER :: ji, jj, jk ! dummy loop indices 565 575 REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars 566 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zwy, zwz, zww ! 2D workspace 567 !!---------------------------------------------------------------------- 568 ! 569 IF( kt == nit000 ) THEN 570 IF(lwp) WRITE(numout,*) 571 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 572 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 576 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace 577 !!---------------------------------------------------------------------- 578 ! 579 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 580 IF( kt == nit000 ) THEN 581 IF(lwp) WRITE(numout,*) 582 IF(lwp) WRITE(numout,*) 'dyn:vor_ens : vorticity term: enstrophy conserving scheme' 583 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 584 ENDIF 573 585 ENDIF 574 586 ! ! =============== … … 578 590 SELECT CASE( kvor ) !== vorticity considered ==! 579 591 CASE ( np_COR ) !* Coriolis (planetary vorticity) 580 zwz(:,:) = ff_f(:,:) 592 DO_2D( 1, 0, 1, 0 ) 593 zwz(ji,jj) = ff_f(ji,jj) 594 END_2D 581 595 CASE ( np_RVO ) !* relative vorticity 582 596 DO_2D( 1, 0, 1, 0 ) … … 645 659 #endif 646 660 ! !== horizontal fluxes ==! 647 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 648 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 661 DO_2D( 1, 1, 1, 1 ) 662 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 663 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 664 END_2D 649 665 ! 650 666 ! !== compute and add the vorticity term trend =! … … 690 706 REAL(wp) :: zua, zva ! local scalars 691 707 REAL(wp) :: zmsk, ze3f ! local scalars 692 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy , z1_e3f 693 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 694 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 695 !!---------------------------------------------------------------------- 696 ! 697 IF( kt == nit000 ) THEN 698 IF(lwp) WRITE(numout,*) 699 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 700 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 708 REAL(wp), DIMENSION(A2D(nn_hls)) :: z1_e3f 709 #if defined key_loop_fusion 710 REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1 711 REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1 712 REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1 713 #else 714 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy 715 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse 716 #endif 717 REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined 718 !!---------------------------------------------------------------------- 719 ! 720 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 721 IF( kt == nit000 ) THEN 722 IF(lwp) WRITE(numout,*) 723 IF(lwp) WRITE(numout,*) 'dyn:vor_een : vorticity term: energy and enstrophy conserving scheme' 724 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 725 ENDIF 701 726 ENDIF 702 727 ! … … 706 731 ! 707 732 #if defined key_qco || defined key_linssh 708 DO_2D( 1, 0, 1, 0) ! == reciprocal of e3 at F-point (key_qco)733 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! == reciprocal of e3 at F-point (key_qco) 709 734 z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk) 710 735 END_2D … … 712 737 SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point 713 738 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 714 DO_2D( 1, 0, 1, 0 ) 715 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 716 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 717 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 718 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 739 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 740 ! round brackets added to fix the order of floating point operations 741 ! needed to ensure halo 1 - halo 2 compatibility 742 ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 743 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) & 744 & + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 745 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) ) 719 746 IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f 720 747 ELSE ; z1_e3f(ji,jj) = 0._wp … … 722 749 END_2D 723 750 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 724 DO_2D( 1, 0, 1, 0 ) 725 ze3f = ( e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 726 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk) & 727 & + e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 728 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk) ) 751 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 752 ! round brackets added to fix the order of floating point operations 753 ! needed to ensure halo 1 - halo 2 compatibility 754 ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) & 755 & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) & 756 & + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) & 757 & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) ) 729 758 zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) & 730 759 & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) ) … … 739 768 ! 740 769 CASE ( np_COR ) !* Coriolis (planetary vorticity) 741 DO_2D( 1, 0, 1, 0)770 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 742 771 zwz(ji,jj,jk) = ff_f(ji,jj) * z1_e3f(ji,jj) 743 772 END_2D 744 773 CASE ( np_RVO ) !* relative vorticity 745 DO_2D( 1, 0, 1, 0)774 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 746 775 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 747 776 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj) 748 777 END_2D 749 778 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 750 DO_2D( 1, 0, 1, 0)779 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 751 780 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 752 781 END_2D 753 782 ENDIF 754 783 CASE ( np_MET ) !* metric term 755 DO_2D( 1, 0, 1, 0)784 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 756 785 zwz(ji,jj,jk) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 757 786 & - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) 758 787 END_2D 759 788 CASE ( np_CRV ) !* Coriolis + relative vorticity 760 DO_2D( 1, 0, 1, 0 ) 761 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 762 & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & 763 & * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 789 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 790 ! round brackets added to fix the order of floating point operations 791 ! needed to ensure halo 1 - halo 2 compatibility 792 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 793 & ) & ! bracket for halo 1 - halo 2 compatibility 794 & - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) & 795 & ) & ! bracket for halo 1 - halo 2 compatibility 796 & ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj) 764 797 END_2D 765 798 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 766 DO_2D( 1, 0, 1, 0)799 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 767 800 zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 768 801 END_2D 769 802 ENDIF 770 803 CASE ( np_CME ) !* Coriolis + metric 771 DO_2D( 1, 0, 1, 0)804 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 772 805 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 773 806 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj) … … 780 813 ! ! =============== 781 814 ! 782 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 783 ! 784 ! ! =============== 785 DO jk = 1, jpkm1 ! Horizontal slab 786 ! ! =============== 787 ! 815 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 816 ! 817 ! ! =============== 818 ! ! Horizontal slab 819 ! ! =============== 820 #if defined key_loop_fusion 821 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 788 822 ! !== horizontal fluxes ==! 789 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 790 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 823 zwx = e2u(ji ,jj ) * e3u(ji ,jj ,jk,Kmm) * pu(ji ,jj ,jk) 824 zwx_im1 = e2u(ji-1,jj ) * e3u(ji-1,jj ,jk,Kmm) * pu(ji-1,jj ,jk) 825 zwx_jp1 = e2u(ji ,jj+1) * e3u(ji ,jj+1,jk,Kmm) * pu(ji ,jj+1,jk) 826 zwx_im1_jp1 = e2u(ji-1,jj+1) * e3u(ji-1,jj+1,jk,Kmm) * pu(ji-1,jj+1,jk) 827 zwy = e1v(ji ,jj ) * e3v(ji ,jj ,jk,Kmm) * pv(ji ,jj ,jk) 828 zwy_ip1 = e1v(ji+1,jj ) * e3v(ji+1,jj ,jk,Kmm) * pv(ji+1,jj ,jk) 829 zwy_jm1 = e1v(ji ,jj-1) * e3v(ji ,jj-1,jk,Kmm) * pv(ji ,jj-1,jk) 830 zwy_ip1_jm1 = e1v(ji+1,jj-1) * e3v(ji+1,jj-1,jk,Kmm) * pv(ji+1,jj-1,jk) 831 ! !== compute and add the vorticity term trend =! 832 ztne = zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) 833 ztnw = zwz(ji-1,jj-1,jk) + zwz(ji-1,jj ,jk) + zwz(ji ,jj ,jk) 834 ztnw_ip1 = zwz(ji ,jj-1,jk) + zwz(ji ,jj ,jk) + zwz(ji+1,jj ,jk) 835 ztse = zwz(ji ,jj ,jk) + zwz(ji ,jj-1,jk) + zwz(ji-1,jj-1,jk) 836 ztse_jp1 = zwz(ji ,jj+1,jk) + zwz(ji ,jj ,jk) + zwz(ji-1,jj ,jk) 837 ztsw_jp1 = zwz(ji ,jj ,jk) + zwz(ji-1,jj ,jk) + zwz(ji-1,jj+1,jk) 838 ztsw_ip1 = zwz(ji+1,jj-1,jk) + zwz(ji ,jj-1,jk) + zwz(ji ,jj ,jk) 839 ! 840 zua = + r1_12 * r1_e1u(ji,jj) * ( ztne * zwy + ztnw_ip1 * zwy_ip1 & 841 & + ztse * zwy_jm1 + ztsw_ip1 * zwy_ip1_jm1 ) 842 zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw_jp1 * zwx_im1_jp1 + ztse_jp1 * zwx_jp1 & 843 & + ztnw * zwx_im1 + ztne * zwx ) 844 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua 845 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 846 END_3D 847 #else 848 DO jk = 1, jpkm1 849 ! 850 ! !== horizontal fluxes ==! 851 DO_2D( 1, 1, 1, 1 ) 852 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 853 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 854 END_2D 791 855 ! 792 856 ! !== compute and add the vorticity term trend =! … … 806 870 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva 807 871 END_2D 808 ! ! =============== 809 END DO ! End of slab 872 END DO 873 #endif 874 ! ! =============== 875 ! ! End of slab 810 876 ! ! =============== 811 877 END SUBROUTINE vor_een … … 839 905 REAL(wp) :: zua, zva ! local scalars 840 906 REAL(wp) :: zmsk, z1_e3t ! local scalars 841 REAL(wp), DIMENSION(jpi,jpj) :: zwx , zwy 842 REAL(wp), DIMENSION(jpi,jpj) :: ztnw, ztne, ztsw, ztse 843 REAL(wp), DIMENSION(jpi,jpj,jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 844 !!---------------------------------------------------------------------- 845 ! 846 IF( kt == nit000 ) THEN 847 IF(lwp) WRITE(numout,*) 848 IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme' 849 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 907 REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy 908 REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse 909 REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined 910 !!---------------------------------------------------------------------- 911 ! 912 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 913 IF( kt == nit000 ) THEN 914 IF(lwp) WRITE(numout,*) 915 IF(lwp) WRITE(numout,*) 'dyn:vor_eeT : vorticity term: energy and enstrophy conserving scheme' 916 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 917 ENDIF 850 918 ENDIF 851 919 ! … … 857 925 SELECT CASE( kvor ) !== vorticity considered ==! 858 926 CASE ( np_COR ) !* Coriolis (planetary vorticity) 859 DO_2D( 1, 0, 1, 0)927 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 860 928 zwz(ji,jj,jk) = ff_f(ji,jj) 861 929 END_2D 862 930 CASE ( np_RVO ) !* relative vorticity 863 DO_2D( 1, 0, 1, 0 ) 864 zwz(ji,jj,jk) = ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 865 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & 931 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 932 ! round brackets added to fix the order of floating point operations 933 ! needed to ensure halo 1 - halo 2 compatibility 934 zwz(ji,jj,jk) = ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) & 935 & - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) & 866 936 & * r1_e1e2f(ji,jj) 867 937 END_2D 868 938 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 869 DO_2D( 1, 0, 1, 0)939 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 870 940 zwz(ji,jj,jk) = zwz(ji,jj,jk) * fmask(ji,jj,jk) 871 941 END_2D 872 942 ENDIF 873 943 CASE ( np_MET ) !* metric term 874 DO_2D( 1, 0, 1, 0)944 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 875 945 zwz(ji,jj,jk) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 876 946 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) 877 947 END_2D 878 948 CASE ( np_CRV ) !* Coriolis + relative vorticity 879 DO_2D( 1, 0, 1, 0 ) 880 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk) & 881 & - e1u(ji ,jj+1) * pu(ji ,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) & 882 & * r1_e1e2f(ji,jj) ) 949 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 950 ! round brackets added to fix the order of floating point operations 951 ! needed to ensure halo 1 - halo 2 compatibility 952 zwz(ji,jj,jk) = ( ff_f(ji,jj) + ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) & 953 & - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) & 954 & * r1_e1e2f(ji,jj) ) 883 955 END_2D 884 956 IF( ln_dynvor_msk ) THEN ! mask the relative vorticity 885 DO_2D( 1, 0, 1, 0)957 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 886 958 zwz(ji,jj,jk) = ( zwz(ji,jj,jk) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj) 887 959 END_2D 888 960 ENDIF 889 961 CASE ( np_CME ) !* Coriolis + metric 890 DO_2D( 1, 0, 1, 0)962 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 891 963 zwz(ji,jj,jk) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) & 892 964 & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) … … 900 972 ! ! =============== 901 973 ! 902 CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp )974 IF (nn_hls==1) CALL lbc_lnk( 'dynvor', zwz, 'F', 1.0_wp ) 903 975 ! 904 976 ! ! =============== … … 907 979 ! 908 980 ! !== horizontal fluxes ==! 909 zwx(:,:) = e2u(:,:) * e3u(:,:,jk,Kmm) * pu(:,:,jk) 910 zwy(:,:) = e1v(:,:) * e3v(:,:,jk,Kmm) * pv(:,:,jk) 981 DO_2D( 1, 1, 1, 1 ) 982 zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk) 983 zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk) 984 END_2D 911 985 ! 912 986 ! !== compute and add the vorticity term trend =! … … 1021 1095 dj_e1v_2(ji,jj) = ( e1v(ji,jj) - e1v(ji ,jj-1) ) * 0.5_wp 1022 1096 END_2D 1023 CALL lbc_lnk _multi( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp ) ! Lateral boundary conditions1097 CALL lbc_lnk( 'dynvor', di_e2u_2, 'T', -1.0_wp , dj_e1v_2, 'T', -1.0_wp ) ! Lateral boundary conditions 1024 1098 ! 1025 1099 CASE DEFAULT !* F-point metric term : pre-compute di(e2u)/(2*e1e2f) and dj(e1v)/(2*e1e2f) … … 1029 1103 dj_e1u_2e1e2f(ji,jj) = ( e1u(ji ,jj+1) - e1u(ji,jj) ) * 0.5 * r1_e1e2f(ji,jj) 1030 1104 END_2D 1031 CALL lbc_lnk _multi( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp ) ! Lateral boundary conditions1105 CALL lbc_lnk( 'dynvor', di_e2v_2e1e2f, 'F', -1.0_wp , dj_e1u_2e1e2f, 'F', -1.0_wp ) ! Lateral boundary conditions 1032 1106 END SELECT 1033 1107 !
Note: See TracChangeset
for help on using the changeset viewer.
