- 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/ICE/icedyn_adv_umx.F90 (modified) (55 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/ICE/icedyn_adv_umx.F90
r14215 r15574 119 119 CALL icemax3D( ph_ip, zhip_max) 120 120 CALL icemax3D( zs_i , zsi_max ) 121 CALL lbc_lnk _multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp )121 CALL lbc_lnk( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 122 122 ! 123 123 ! enthalpies … … 164 164 ! 165 165 ! --- define velocity for advection: u*grad(H) --- ! 166 DO_2D( 0, 0, 0, 0)166 DO_2D( nn_hls-1, nn_hls, nn_hls, nn_hls ) 167 167 IF ( pu_ice(ji,jj) * pu_ice(ji-1,jj) <= 0._wp ) THEN ; zcu_box(ji,jj) = 0._wp 168 168 ELSEIF( pu_ice(ji,jj) > 0._wp ) THEN ; zcu_box(ji,jj) = pu_ice(ji-1,jj) 169 169 ELSE ; zcu_box(ji,jj) = pu_ice(ji ,jj) 170 170 ENDIF 171 171 END_2D 172 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls ) 172 173 IF ( pv_ice(ji,jj) * pv_ice(ji,jj-1) <= 0._wp ) THEN ; zcv_box(ji,jj) = 0._wp 173 174 ELSEIF( pv_ice(ji,jj) > 0._wp ) THEN ; zcv_box(ji,jj) = pv_ice(ji,jj-1) … … 204 205 IF( .NOT. ALLOCATED(jmsk_small) ) ALLOCATE( jmsk_small(jpi,jpj,jpl) ) 205 206 DO jl = 1, jpl 206 DO_2D( 1, 0, 1, 0)207 DO_2D( 1, 0, nn_hls, nn_hls ) 207 208 zvi_cen = 0.5_wp * ( pv_i(ji+1,jj,jl) + pv_i(ji,jj,jl) ) 208 209 IF( zvi_cen < epsi06) THEN ; imsk_small(ji,jj,jl) = 0 209 210 ELSE ; imsk_small(ji,jj,jl) = 1 ; ENDIF 211 END_2D 212 DO_2D( nn_hls, nn_hls, 1, 0 ) 210 213 zvi_cen = 0.5_wp * ( pv_i(ji,jj+1,jl) + pv_i(ji,jj,jl) ) 211 214 IF( zvi_cen < epsi06) THEN ; jmsk_small(ji,jj,jl) = 0 … … 360 363 ! --- Lateral boundary conditions --- ! 361 364 IF ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 362 CALL lbc_lnk _multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &363 & 365 CALL lbc_lnk( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 366 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 364 367 ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 365 CALL lbc_lnk _multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp &366 & 368 CALL lbc_lnk( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 369 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 367 370 ELSE 368 CALL lbc_lnk _multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp )371 CALL lbc_lnk( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 369 372 ENDIF 370 373 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) … … 583 586 ! 584 587 DO jl = 1, jpl 585 DO_2D( 1, 0, 1, 0)588 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) 586 589 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) 587 590 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) … … 594 597 ! 595 598 DO jl = 1, jpl !-- flux in x-direction 596 DO_2D( 1, 0, 1, 1)599 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls ) 597 600 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) 598 601 END_2D … … 600 603 ! 601 604 DO jl = 1, jpl !-- first guess of tracer from u-flux 602 DO_2D( 0, 0, 1, 1)605 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls ) 603 606 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 604 607 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 609 612 ! 610 613 DO jl = 1, jpl !-- flux in y-direction 611 DO_2D( 0, 0, 1, 0)614 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 612 615 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) 613 616 END_2D … … 617 620 ! 618 621 DO jl = 1, jpl !-- flux in y-direction 619 DO_2D( 1, 1, 1, 0)622 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 ) 620 623 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) 621 624 END_2D … … 623 626 ! 624 627 DO jl = 1, jpl !-- first guess of tracer from v-flux 625 DO_2D( 1, 1, 0, 0)628 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) 626 629 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 627 630 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 632 635 ! 633 636 DO jl = 1, jpl !-- flux in x-direction 634 DO_2D( 1, 0, 0, 0)637 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 635 638 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) 636 639 END_2D … … 642 645 ! 643 646 DO jl = 1, jpl !-- after tracer with upstream scheme 644 DO_2D( 0, 0, 0, 0)647 DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) 645 648 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj ,jl) & 646 649 & + pfv_ups(ji,jj,jl) - pfv_ups(ji ,jj-1,jl) ) & … … 651 654 END_2D 652 655 END DO 653 CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp )656 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', pt_ups, 'T', 1.0_wp ) 654 657 655 658 END SUBROUTINE upstream … … 681 684 ! 682 685 DO jl = 1, jpl 683 DO_2D( 1, 0, 1, 1)686 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls ) 684 687 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 685 688 END_2D 686 DO_2D( 1, 1, 1, 0)689 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 ) 687 690 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 688 691 END_2D … … 701 704 ! 702 705 DO jl = 1, jpl !-- flux in x-direction 703 DO_2D( 1, 0, 1, 1)706 DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls ) 704 707 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 705 708 END_2D … … 708 711 709 712 DO jl = 1, jpl !-- first guess of tracer from u-flux 710 DO_2D( 0, 0, 1, 1)713 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls ) 711 714 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 712 715 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 717 720 718 721 DO jl = 1, jpl !-- flux in y-direction 719 DO_2D( 0, 0, 1, 0)722 DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 ) 720 723 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 721 724 END_2D … … 726 729 ! 727 730 DO jl = 1, jpl !-- flux in y-direction 728 DO_2D( 1, 1, 1, 0)731 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls-1 ) 729 732 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 730 733 END_2D … … 733 736 ! 734 737 DO jl = 1, jpl !-- first guess of tracer from v-flux 735 DO_2D( 1, 1, 0, 0)738 DO_2D( nn_hls, nn_hls, nn_hls-1, nn_hls-1 ) 736 739 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 737 740 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 742 745 ! 743 746 DO jl = 1, jpl !-- flux in x-direction 744 DO_2D( 1, 0, 0, 0)747 DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 ) 745 748 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 746 749 END_2D … … 785 788 ! 786 789 ! !-- ultimate interpolation of pt at u-point --! 787 CALL ultimate_x( pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho )790 CALL ultimate_x( nn_hls, pamsk, kn_umx, pdt, pt, pu, zt_u, pfu_ho ) 788 791 ! !-- limiter in x --! 789 792 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_x( pdt, pu, pt, pfu_ups, pfu_ho ) 790 793 ! !-- advective form update in zpt --! 791 794 DO jl = 1, jpl 792 DO_2D( 0, 0, 0, 0)795 DO_2D( 0, 0, nn_hls, nn_hls ) 793 796 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) & 794 797 & + pt (ji,jj,jl) * ( pu (ji,jj ) - pu (ji-1,jj ) ) * r1_e1e2t(ji,jj) & … … 797 800 END_2D 798 801 END DO 799 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )800 802 ! 801 803 ! !-- ultimate interpolation of pt at v-point --! 802 804 IF( ll_hoxy ) THEN 803 CALL ultimate_y( pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho )805 CALL ultimate_y( 0, pamsk, kn_umx, pdt, zpt, pv, zt_v, pfv_ho ) 804 806 ELSE 805 CALL ultimate_y( pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho )807 CALL ultimate_y( 0, pamsk, kn_umx, pdt, pt , pv, zt_v, pfv_ho ) 806 808 ENDIF 807 809 ! !-- limiter in y --! … … 812 814 ! 813 815 ! !-- ultimate interpolation of pt at v-point --! 814 CALL ultimate_y( pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho )816 CALL ultimate_y( nn_hls, pamsk, kn_umx, pdt, pt, pv, zt_v, pfv_ho ) 815 817 ! !-- limiter in y --! 816 818 IF( np_limiter == 2 .OR. np_limiter == 3 ) CALL limiter_y( pdt, pv, pt, pfv_ups, pfv_ho ) 817 819 ! !-- advective form update in zpt --! 818 820 DO jl = 1, jpl 819 DO_2D( 0, 0, 0, 0 )821 DO_2D( nn_hls, nn_hls, 0, 0 ) 820 822 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) & 821 823 & + pt (ji,jj,jl) * ( pv (ji,jj ) - pv (ji,jj-1 ) ) * r1_e1e2t(ji,jj) & … … 824 826 END_2D 825 827 END DO 826 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )827 828 ! 828 829 ! !-- ultimate interpolation of pt at u-point --! 829 830 IF( ll_hoxy ) THEN 830 CALL ultimate_x( pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho )831 CALL ultimate_x( 0, pamsk, kn_umx, pdt, zpt, pu, zt_u, pfu_ho ) 831 832 ELSE 832 CALL ultimate_x( pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho )833 CALL ultimate_x( 0, pamsk, kn_umx, pdt, pt , pu, zt_u, pfu_ho ) 833 834 ENDIF 834 835 ! !-- limiter in x --! … … 842 843 843 844 844 SUBROUTINE ultimate_x( pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho )845 SUBROUTINE ultimate_x( kloop, pamsk, kn_umx, pdt, pt, pu, pt_u, pfu_ho ) 845 846 !!--------------------------------------------------------------------- 846 847 !! *** ROUTINE ultimate_x *** … … 852 853 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 853 854 !!---------------------------------------------------------------------- 855 INTEGER , INTENT(in ) :: kloop ! either 0 or nn_hls depending on the order of the call 854 856 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 855 857 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) … … 867 869 ! !-- Laplacian in i-direction --! 868 870 DO jl = 1, jpl 869 DO jj = 2, jpjm1 ! First derivative (gradient) 870 DO ji = 1, jpim1 871 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 872 END DO 873 ! ! Second derivative (Laplacian) 874 DO ji = 2, jpim1 875 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 876 END DO 877 END DO 878 END DO 879 CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp ) 871 DO_2D( nn_hls, nn_hls-1, kloop, kloop ) ! First derivative (gradient) 872 ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 873 END_2D 874 DO_2D( nn_hls-1, nn_hls-1, kloop, kloop ) ! Second derivative (Laplacian) 875 ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 876 END_2D 877 !!$ DO jj = 2, jpjm1 ! First derivative (gradient) 878 !!$ DO ji = 1, jpim1 879 !!$ ztu1(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 880 !!$ END DO 881 !!$ ! ! Second derivative (Laplacian) 882 !!$ DO ji = 2, jpim1 883 !!$ ztu2(ji,jj,jl) = ( ztu1(ji,jj,jl) - ztu1(ji-1,jj,jl) ) * r1_e1t(ji,jj) 884 !!$ END DO 885 !!$ END DO 886 END DO 887 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', ztu2, 'T', 1.0_wp ) 880 888 ! 881 889 ! !-- BiLaplacian in i-direction --! 882 890 DO jl = 1, jpl 883 DO jj = 2, jpjm1 ! Third derivative 884 DO ji = 1, jpim1 885 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 886 END DO 887 ! ! Fourth derivative 888 DO ji = 2, jpim1 889 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 890 END DO 891 END DO 892 END DO 893 CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp ) 891 DO_2D( 1, 0, kloop, kloop ) ! Third derivative 892 ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 893 END_2D 894 DO_2D( 0, 0, kloop, kloop ) ! Fourth derivative 895 ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 896 END_2D 897 !!$ DO jj = 2, jpjm1 ! Third derivative 898 !!$ DO ji = 1, jpim1 899 !!$ ztu3(ji,jj,jl) = ( ztu2(ji+1,jj,jl) - ztu2(ji,jj,jl) ) * r1_e1u(ji,jj) * umask(ji,jj,1) 900 !!$ END DO 901 !!$ ! ! Fourth derivative 902 !!$ DO ji = 2, jpim1 903 !!$ ztu4(ji,jj,jl) = ( ztu3(ji,jj,jl) - ztu3(ji-1,jj,jl) ) * r1_e1t(ji,jj) 904 !!$ END DO 905 !!$ END DO 906 END DO 894 907 ! 895 908 ! … … 899 912 ! 900 913 DO jl = 1, jpl 901 DO_2D( 1, 0, 0, 0)914 DO_2D( 1, 0, kloop, kloop ) 902 915 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 903 916 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 908 921 ! 909 922 DO jl = 1, jpl 910 DO_2D( 1, 0, 0, 0)923 DO_2D( 1, 0, kloop, kloop ) 911 924 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 912 925 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 918 931 ! 919 932 DO jl = 1, jpl 920 DO_2D( 1, 0, 0, 0)933 DO_2D( 1, 0, kloop, kloop ) 921 934 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 922 935 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 932 945 ! 933 946 DO jl = 1, jpl 934 DO_2D( 1, 0, 0, 0)947 DO_2D( 1, 0, kloop, kloop ) 935 948 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 936 949 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 945 958 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 946 959 ! 947 DO jl = 1, jpl 948 DO_2D( 1, 0, 0, 0 ) 960 CALL lbc_lnk( 'icedyn_adv_umx', ztu4, 'T', 1.0_wp ) 961 ! 962 DO jl = 1, jpl 963 DO_2D( 1, 0, kloop, kloop ) 949 964 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 950 965 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 967 982 IF( ll_neg ) THEN 968 983 DO jl = 1, jpl 969 DO_2D( 1, 0, 0, 0)984 DO_2D( 1, 0, kloop, kloop ) 970 985 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 971 986 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 985 1000 986 1001 987 SUBROUTINE ultimate_y( pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho )1002 SUBROUTINE ultimate_y( kloop, pamsk, kn_umx, pdt, pt, pv, pt_v, pfv_ho ) 988 1003 !!--------------------------------------------------------------------- 989 1004 !! *** ROUTINE ultimate_y *** … … 995 1010 !! Reference : Leonard, B.P., 1991, Comput. Methods Appl. Mech. Eng., 88, 17-74. 996 1011 !!---------------------------------------------------------------------- 1012 INTEGER , INTENT(in ) :: kloop ! either 0 or nn_hls depending on the order of the call 997 1013 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) 998 1014 INTEGER , INTENT(in ) :: kn_umx ! order of the scheme (1-5=UM or 20=CEN2) … … 1010 1026 ! !-- Laplacian in j-direction --! 1011 1027 DO jl = 1, jpl 1012 DO_2D( 0, 0, 1, 0 )! First derivative (gradient)1028 DO_2D( kloop, kloop, nn_hls, nn_hls-1 ) ! First derivative (gradient) 1013 1029 ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1014 1030 END_2D 1015 DO_2D( 0, 0, 0, 0 )! Second derivative (Laplacian)1031 DO_2D( kloop, kloop, nn_hls-1, nn_hls-1 ) ! Second derivative (Laplacian) 1016 1032 ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1017 1033 END_2D 1018 1034 END DO 1019 CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp )1035 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', ztv2, 'T', 1.0_wp ) 1020 1036 ! 1021 1037 ! !-- BiLaplacian in j-direction --! 1022 1038 DO jl = 1, jpl 1023 DO_2D( 0, 0, 1, 0 ) ! Firstderivative1039 DO_2D( kloop, kloop, 1, 0 ) ! Third derivative 1024 1040 ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 1025 1041 END_2D 1026 DO_2D( 0, 0, 0, 0 ) ! Secondderivative1042 DO_2D( kloop, kloop, 0, 0 ) ! Fourth derivative 1027 1043 ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 1028 1044 END_2D 1029 1045 END DO 1030 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp )1031 1046 ! 1032 1047 ! … … 1035 1050 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1036 1051 DO jl = 1, jpl 1037 DO_2D( 0, 0, 1, 0 )1052 DO_2D( kloop, kloop, 1, 0 ) 1038 1053 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1039 1054 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1043 1058 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1044 1059 DO jl = 1, jpl 1045 DO_2D( 0, 0, 1, 0 )1060 DO_2D( kloop, kloop, 1, 0 ) 1046 1061 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1047 1062 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1052 1067 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1053 1068 DO jl = 1, jpl 1054 DO_2D( 0, 0, 1, 0 )1069 DO_2D( kloop, kloop, 1, 0 ) 1055 1070 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1056 1071 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1065 1080 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1066 1081 DO jl = 1, jpl 1067 DO_2D( 0, 0, 1, 0 )1082 DO_2D( kloop, kloop, 1, 0 ) 1068 1083 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1069 1084 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1077 1092 ! 1078 1093 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1079 DO jl = 1, jpl 1080 DO_2D( 0, 0, 1, 0 ) 1094 ! 1095 CALL lbc_lnk( 'icedyn_adv_umx', ztv4, 'T', 1.0_wp ) 1096 ! 1097 DO jl = 1, jpl 1098 DO_2D( kloop, kloop, 1, 0 ) 1081 1099 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1082 1100 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1099 1117 IF( ll_neg ) THEN 1100 1118 DO jl = 1, jpl 1101 DO_2D( 0, 0, 1, 0 )1119 DO_2D( kloop, kloop, 1, 0 ) 1102 1120 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1103 1121 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1169 1187 END_2D 1170 1188 END DO 1171 CALL lbc_lnk _multi( 'icedyn_adv_umx', zti_ups, 'T', 1.0_wp, ztj_ups, 'T', 1.0_wp )1189 CALL lbc_lnk( 'icedyn_adv_umx', zti_ups, 'T', 1.0_wp, ztj_ups, 'T', 1.0_wp ) 1172 1190 1173 1191 DO jl = 1, jpl … … 1191 1209 END_2D 1192 1210 END DO 1193 CALL lbc_lnk _multi( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp, pfv_ho, 'V', -1.0_wp ) ! lateral boundary cond.1211 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp, pfv_ho, 'V', -1.0_wp ) ! lateral boundary cond. 1194 1212 1195 1213 ENDIF … … 1248 1266 END_2D 1249 1267 END DO 1250 CALL lbc_lnk _multi( 'icedyn_adv_umx', zbetup, 'T', 1.0_wp, zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign)1268 CALL lbc_lnk( 'icedyn_adv_umx', zbetup, 'T', 1.0_wp, zbetdo, 'T', 1.0_wp ) ! lateral boundary cond. (unchanged sign) 1251 1269 1252 1270 … … 1299 1317 ! 1300 1318 DO jl = 1, jpl 1301 DO_2D( 0, 0, 0, 0 )1319 DO_2D( nn_hls, nn_hls-1, 0, 0 ) 1302 1320 zslpx(ji,jj,jl) = ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) * umask(ji,jj,1) 1303 1321 END_2D 1304 1322 END DO 1305 CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp) ! lateral boundary cond.1306 1307 DO jl = 1, jpl 1308 DO_2D( 0, 0, 0, 0 )1323 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', zslpx, 'U', -1.0_wp) ! lateral boundary cond. 1324 1325 DO jl = 1, jpl 1326 DO_2D( nn_hls-1, 0, 0, 0 ) 1309 1327 uCFL = pdt * ABS( pu(ji,jj) ) * r1_e1e2t(ji,jj) 1310 1328 … … 1367 1385 END_2D 1368 1386 END DO 1369 CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp) ! lateral boundary cond.1387 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', pfu_ho, 'U', -1.0_wp) ! lateral boundary cond. 1370 1388 ! 1371 1389 END SUBROUTINE limiter_x … … 1390 1408 ! 1391 1409 DO jl = 1, jpl 1392 DO_2D( 0, 0, 0, 0)1410 DO_2D( 0, 0, nn_hls, nn_hls-1 ) 1393 1411 zslpy(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * vmask(ji,jj,1) 1394 1412 END_2D 1395 1413 END DO 1396 CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp) ! lateral boundary cond.1397 1398 DO jl = 1, jpl 1399 DO_2D( 0, 0, 0, 0 )1414 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', zslpy, 'V', -1.0_wp) ! lateral boundary cond. 1415 1416 DO jl = 1, jpl 1417 DO_2D( 0, 0, nn_hls-1, 0 ) 1400 1418 vCFL = pdt * ABS( pv(ji,jj) ) * r1_e1e2t(ji,jj) 1401 1419 … … 1459 1477 END_2D 1460 1478 END DO 1461 CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp) ! lateral boundary cond.1479 IF( nn_hls == 1 ) CALL lbc_lnk( 'icedyn_adv_umx', pfv_ho, 'V', -1.0_wp) ! lateral boundary cond. 1462 1480 ! 1463 1481 END SUBROUTINE limiter_y … … 1494 1512 ! 1495 1513 DO jl = 1, jpl 1496 DO_2D( 1, 1, 1, 1)1514 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1497 1515 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1498 1516 ! … … 1541 1559 ! ! -- check e_i/v_i -- ! 1542 1560 DO jl = 1, jpl 1543 DO_3D( 1, 1, 1, 1, 1, nlay_i )1561 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_i ) 1544 1562 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1545 1563 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean … … 1555 1573 ! ! -- check e_s/v_s -- ! 1556 1574 DO jl = 1, jpl 1557 DO_3D( 1, 1, 1, 1, 1, nlay_s )1575 DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, nlay_s ) 1558 1576 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 1559 1577 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean … … 1598 1616 ! -- check snow load -- ! 1599 1617 DO jl = 1, jpl 1600 DO_2D( 1, 1, 1, 1)1618 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1601 1619 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 1602 1620 ! … … 1627 1645 !! ** Purpose : compute the max of the 9 points around 1628 1646 !!---------------------------------------------------------------------- 1629 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1630 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1631 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1647 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pice ! input 1648 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pmax ! output 1649 ! 1650 REAL(wp), DIMENSION(Nis0:Nie0) :: zmax1, zmax2 1651 REAL(wp) :: zmax3 1632 1652 INTEGER :: ji, jj, jl ! dummy loop indices 1633 1653 !!---------------------------------------------------------------------- 1634 DO jl = 1, jpl 1635 DO jj = Njs0-1, Nje0+1 1636 DO ji = Nis0, Nie0 1637 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1638 END DO 1639 END DO 1640 DO jj = Njs0, Nje0 1641 DO ji = Nis0, Nie0 1642 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1643 END DO 1644 END DO 1654 ! basic version: get the max of epsi20 + 9 neighbours 1655 !!$ DO jl = 1, jpl 1656 !!$ DO_2D( 0, 0, 0, 0 ) 1657 !!$ pmax(ji,jj,jl) = MAX( epsi20, pice(ji-1,jj-1,jl), pice(ji,jj-1,jl), pice(ji+1,jj-1,jl), & 1658 !!$ & pice(ji-1,jj ,jl), pice(ji,jj ,jl), pice(ji+1,jj ,jl), & 1659 !!$ & pice(ji-1,jj+1,jl), pice(ji,jj+1,jl), pice(ji+1,jj+1,jl) ) 1660 !!$ END_2D 1661 !!$ END DO 1662 ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours 1663 DO jl = 1, jpl 1664 DO ji = Nis0, Nie0 1665 zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jl), pice(ji-1,Njs0-1,jl), pice(ji+1,Njs0-1,jl) ) 1666 zmax2(ji) = MAX( epsi20, pice(ji,Njs0 ,jl), pice(ji-1,Njs0 ,jl), pice(ji+1,Njs0 ,jl) ) 1667 END DO 1668 DO_2D( 0, 0, 0, 0 ) 1669 zmax3 = MAX( epsi20, pice(ji,jj+1,jl), pice(ji-1,jj+1,jl), pice(ji+1,jj+1,jl) ) 1670 pmax(ji,jj,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 ) 1671 zmax1(ji) = zmax2(ji) 1672 zmax2(ji) = zmax3 1673 END_2D 1645 1674 END DO 1646 1675 END SUBROUTINE icemax3D … … 1651 1680 !! ** Purpose : compute the max of the 9 points around 1652 1681 !!---------------------------------------------------------------------- 1653 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1654 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1655 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1682 REAL(wp), DIMENSION(:,:,:,:), INTENT(in ) :: pice ! input 1683 REAL(wp), DIMENSION(:,:,:,:), INTENT(out) :: pmax ! output 1684 ! 1685 REAL(wp), DIMENSION(Nis0:Nie0) :: zmax1, zmax2 1686 REAL(wp) :: zmax3 1656 1687 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1657 1688 !!---------------------------------------------------------------------- 1658 1689 jlay = SIZE( pice , 3 ) ! size of input arrays 1690 ! basic version: get the max of epsi20 + 9 neighbours 1691 !!$ DO jl = 1, jpl 1692 !!$ DO jk = 1, jlay 1693 !!$ DO_2D( 0, 0, 0, 0 ) 1694 !!$ pmax(ji,jj,jk,jl) = MAX( epsi20, pice(ji-1,jj-1,jk,jl), pice(ji,jj-1,jk,jl), pice(ji+1,jj-1,jk,jl), & 1695 !!$ & pice(ji-1,jj ,jk,jl), pice(ji,jj ,jk,jl), pice(ji+1,jj ,jk,jl), & 1696 !!$ & pice(ji-1,jj+1,jk,jl), pice(ji,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) ) 1697 !!$ END_2D 1698 !!$ END DO 1699 !!$ END DO 1700 ! optimized version : does a little bit more than 2 max of epsi20 + 3 neighbours 1659 1701 DO jl = 1, jpl 1660 1702 DO jk = 1, jlay 1661 DO j j = Njs0-1, Nje0+11662 DO ji = Nis0, Nie01663 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) )1664 1665 END DO1666 DO jj = Njs0, Nje01667 DO ji = Nis0, Nie01668 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1))1669 END DO1670 END DO1703 DO ji = Nis0, Nie0 1704 zmax1(ji) = MAX( epsi20, pice(ji,Njs0-1,jk,jl), pice(ji-1,Njs0-1,jk,jl), pice(ji+1,Njs0-1,jk,jl) ) 1705 zmax2(ji) = MAX( epsi20, pice(ji,Njs0 ,jk,jl), pice(ji-1,Njs0 ,jk,jl), pice(ji+1,Njs0 ,jk,jl) ) 1706 END DO 1707 DO_2D( 0, 0, 0, 0 ) 1708 zmax3 = MAX( epsi20, pice(ji,jj+1,jk,jl), pice(ji-1,jj+1,jk,jl), pice(ji+1,jj+1,jk,jl) ) 1709 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax1(ji), zmax2(ji), zmax3 ) 1710 zmax1(ji) = zmax2(ji) 1711 zmax2(ji) = zmax3 1712 END_2D 1671 1713 END DO 1672 1714 END DO
Note: See TracChangeset
for help on using the changeset viewer.
