Changeset 8998 for branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
- Timestamp:
- 2017-12-13T09:34:57+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r8993 r8998 118 118 ! 119 119 IF (.NOT.lk_agrif_clp) THEN 120 DO jk=1,jpkm1 ! Smooth120 DO jk=1,jpkm1 ! Smooth 121 121 DO jj=j1,j2 122 122 ua(2,jj,jk) = 0.25_wp*(ua(1,jj,jk)+2._wp*ua(2,jj,jk)+ua(3,jj,jk)) … … 124 124 END DO 125 125 END DO 126 END 126 ENDIF 127 127 ! 128 128 zub(2,:) = 0._wp ! Correct transport … … 189 189 190 190 IF (.NOT.lk_agrif_clp) THEN 191 DO jk = 1, jpkm1 ! Smooth191 DO jk = 1, jpkm1 ! Smooth 192 192 DO jj = j1, j2 193 193 ua(nlci-2,jj,jk) = 0.25_wp * umask(nlci-2,jj,jk) & … … 196 196 END DO 197 197 ENDIF 198 199 198 zub(nlci-2,:) = 0._wp ! Correct transport 200 199 DO jk = 1, jpkm1 … … 331 330 ! 332 331 IF (.NOT.lk_agrif_clp) THEN 333 DO jk = 1, jpkm1 ! Smooth332 DO jk = 1, jpkm1 ! Smooth 334 333 DO ji = i1, i2 335 334 va(ji,nlcj-2,jk) = 0.25_wp * vmask(ji,nlcj-2,jk) & … … 337 336 END DO 338 337 END DO 339 END IF338 END IF 340 339 ! 341 340 zvb(:,nlcj-2) = 0._wp ! Correct transport … … 614 613 INTEGER , INTENT(in ) :: nb , ndir 615 614 ! 616 INTEGER :: ji, jj, jk, jn ! dummy loop indices617 INTEGER :: imin, imax, jmin, jmax 615 INTEGER :: ji, jj, jk, jn, iref, jref ! dummy loop indices 616 INTEGER :: imin, imax, jmin, jmax, N_in, N_out 618 617 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 619 618 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 620 LOGICAL :: western_side, eastern_side,northern_side,southern_side 621 !!---------------------------------------------------------------------- 622 ! 619 LOGICAL :: western_side, eastern_side,northern_side,southern_side 620 ! vertical interpolation: 621 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 622 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 623 REAL(wp), DIMENSION(k1:k2) :: h_in 624 REAL(wp), DIMENSION(1:jpk) :: h_out(1:jpk) 625 REAL(wp) :: h_diff, zrhoxy 626 627 zrhoxy = Agrif_rhox()*Agrif_rhoy() 623 628 IF (before) THEN 624 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 625 ELSE 629 DO jn = 1,jpts 630 DO jk=k1,k2 631 DO jj=j1,j2 632 DO ji=i1,i2 633 ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 634 END DO 635 END DO 636 END DO 637 END DO 638 639 # if defined key_vertical 640 DO jk=k1,k2 641 DO jj=j1,j2 642 DO ji=i1,i2 643 ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 644 END DO 645 END DO 646 END DO 647 # endif 648 ELSE 649 650 western_side = (nb == 1).AND.(ndir == 1) 651 eastern_side = (nb == 1).AND.(ndir == 2) 652 southern_side = (nb == 2).AND.(ndir == 1) 653 northern_side = (nb == 2).AND.(ndir == 2) 654 655 # if defined key_vertical 656 DO jj=j1,j2 657 DO ji=i1,i2 658 iref = ji 659 jref = jj 660 if(western_side) iref=MAX(2,ji) 661 if(eastern_side) iref=MIN(nlci-1,ji) 662 if(southern_side) jref=MAX(2,jj) 663 if(northern_side) jref=MIN(nlcj-1,jj) 664 N_in = 0 665 DO jk=k1,k2 !k2 = jpk of parent grid 666 IF (ptab(ji,jj,jk,n2) == 0) EXIT 667 N_in = N_in + 1 668 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1) 669 h_in(N_in) = ptab(ji,jj,jk,n2) 670 END DO 671 N_out = 0 672 DO jk=1,jpk ! jpk of child grid 673 IF (tmask(iref,jref,jk) == 0) EXIT 674 N_out = N_out + 1 675 h_out(jk) = e3t_n(iref,jref,jk) 676 ENDDO 677 IF (N_in > 0) THEN 678 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 679 DO jn=1,jpts 680 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 681 ENDDO 682 ENDIF 683 ENDDO 684 ENDDO 685 # else 686 ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts) 687 # endif 688 ! 626 689 IF (lk_agrif_clp) THEN 627 690 DO jn = 1, jpts … … 629 692 DO ji = i1,i2 630 693 DO jj = j1,j2 631 tsa(ji,jj,jk,jn) = ptab (ji,jj,jk,jn)694 tsa(ji,jj,jk,jn) = ptab_child(ji,jj,jk,jn) 632 695 END DO 633 696 END DO … … 636 699 return 637 700 ENDIF 638 !639 western_side = (nb == 1).AND.(ndir == 1)640 eastern_side = (nb == 1).AND.(ndir == 2)641 southern_side = (nb == 2).AND.(ndir == 1)642 northern_side = (nb == 2).AND.(ndir == 2)643 701 ! 644 702 zrhox = Agrif_Rhox() … … 667 725 IF( eastern_side ) THEN 668 726 DO jn = 1, jpts 669 tsa(nlci,j1:j2, k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)727 tsa(nlci,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(nlci,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(nlci-1,j1:j2,1:jpk,jn) 670 728 DO jk = 1, jpkm1 671 729 DO jj = jmin,jmax … … 687 745 IF( northern_side ) THEN 688 746 DO jn = 1, jpts 689 tsa(i1:i2,nlcj, k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)747 tsa(i1:i2,nlcj,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,nlcj,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,nlcj-1,1:jpk,jn) 690 748 DO jk = 1, jpkm1 691 749 DO ji = imin,imax … … 707 765 IF( western_side ) THEN 708 766 DO jn = 1, jpts 709 tsa(1,j1:j2, k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)767 tsa(1,j1:j2,1:jpk,jn) = zalpha1 * ptab_child(1,j1:j2,1:jpk,jn) + zalpha2 * ptab_child(2,j1:j2,1:jpk,jn) 710 768 DO jk = 1, jpkm1 711 769 DO jj = jmin,jmax … … 726 784 IF( southern_side ) THEN 727 785 DO jn = 1, jpts 728 tsa(i1:i2,1, k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)786 tsa(i1:i2,1,1:jpk,jn) = zalpha1 * ptab_child(i1:i2,1,1:jpk,jn) + zalpha2 * ptab_child(i1:i2,2,1:jpk,jn) 729 787 DO jk = 1, jpk 730 788 DO ji=imin,imax … … 747 805 ! East south 748 806 IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 749 tsa(nlci-1,2,:,:) = ptab (nlci-1,2,:,:)807 tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts) 750 808 ENDIF 751 809 ! East north 752 810 IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 753 tsa(nlci-1,nlcj-1,:,:) = ptab (nlci-1,nlcj-1,:,:)811 tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts) 754 812 ENDIF 755 813 ! West south 756 814 IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 757 tsa(2,2,:,:) = ptab (2,2,:,:)815 tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts) 758 816 ENDIF 759 817 ! West north 760 818 IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 761 tsa(2,nlcj-1,:,:) = ptab (2,nlcj-1,:,:)819 tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts) 762 820 ENDIF 763 821 ! … … 765 823 ! 766 824 END SUBROUTINE interptsn 767 768 825 769 826 SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir ) … … 794 851 END SUBROUTINE interpsshn 795 852 796 797 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, before ) 853 SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 798 854 !!---------------------------------------------------------------------- 799 855 !! *** ROUTINE interpun *** 800 !!---------------------------------------------------------------------- 801 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 802 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 803 LOGICAL , INTENT(in ) :: before 804 ! 805 INTEGER :: ji, jj, jk 806 REAL(wp) :: zrhoy 807 !!---------------------------------------------------------------------- 808 ! 809 IF( before ) THEN 810 DO jk = 1, jpkm1 811 ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 856 !!--------------------------------------------- 857 !! 858 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 859 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 860 LOGICAL, INTENT(in) :: before 861 INTEGER, INTENT(in) :: nb , ndir 862 !! 863 INTEGER :: ji,jj,jk 864 REAL(wp) :: zrhoy 865 ! vertical interpolation: 866 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 867 REAL(wp), DIMENSION(1:jpk) :: h_out 868 INTEGER :: N_in, N_out, iref 869 REAL(wp) :: h_diff 870 LOGICAL :: western_side, eastern_side 871 !!--------------------------------------------- 872 ! 873 zrhoy = Agrif_rhoy() 874 IF (before) THEN 875 !We can't use zero as the special value because we need to include zeros 876 !when interpolating the scale factors 877 IF(Agrif_UseSpecialValue) THEN 878 ! Agrif_SpecialValue = -999._wp 879 Agrif_SpecialValue = 0._wp 880 ELSE 881 Agrif_SpecialValue = 0._wp 882 ENDIF 883 DO jk=1,jpk 884 DO jj=j1,j2 885 DO ji=i1,i2 886 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - & 887 & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue) 888 # if defined key_vertical 889 ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 890 # endif 891 END DO 892 END DO 812 893 END DO 813 894 ELSE 814 zrhoy = Agrif_Rhoy() 895 zrhoy = Agrif_rhoy() 896 # if defined key_vertical 897 ! VERTICAL REFINEMENT BEGIN 898 western_side = (nb == 1).AND.(ndir == 1) 899 eastern_side = (nb == 1).AND.(ndir == 2) 900 901 Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed 902 903 DO ji=i1,i2 904 iref = ji 905 IF (western_side) iref = MAX(2,ji) 906 IF (eastern_side) iref = MIN(nlci-2,ji) 907 DO jj=j1,j2 908 N_in = 0 909 DO jk=k1,k2 910 IF (ptab(ji,jj,jk,2) == 0) EXIT 911 N_in = N_in + 1 912 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 913 h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy) 914 ENDDO 915 916 IF (N_in == 0) THEN 917 ua(ji,jj,:) = 0._wp 918 CYCLE 919 ENDIF 920 921 N_out = 0 922 DO jk=1,jpk 923 if (umask(iref,jj,jk) == 0) EXIT 924 N_out = N_out + 1 925 h_out(N_out) = e3u_a(iref,jj,jk) 926 ENDDO 927 928 IF (N_out == 0) THEN 929 ua(ji,jj,:) = 0._wp 930 CYCLE 931 ENDIF 932 933 IF (N_in * N_out > 0) THEN 934 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 935 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 936 if (h_diff < -1.e4) then 937 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 938 ! stop 939 endif 940 ENDIF 941 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ua(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 942 ENDDO 943 ENDDO 944 945 # else 815 946 DO jk = 1, jpkm1 816 947 DO jj=j1,j2 817 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 818 END DO 819 END DO 948 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u_a(i1:i2,jj,jk) ) 949 END DO 950 END DO 951 # endif 952 820 953 ENDIF 821 954 ! 822 955 END SUBROUTINE interpun 823 956 824 825 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, before ) 957 SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir ) 826 958 !!---------------------------------------------------------------------- 827 959 !! *** ROUTINE interpvn *** 828 960 !!---------------------------------------------------------------------- 829 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 830 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 831 LOGICAL , INTENT(in ) :: before 832 ! 833 INTEGER :: ji, jj, jk 834 REAL(wp) :: zrhox 835 !!---------------------------------------------------------------------- 961 ! 962 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 963 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 964 LOGICAL, INTENT(in) :: before 965 INTEGER, INTENT(in) :: nb , ndir 966 ! 967 INTEGER :: ji,jj,jk 968 REAL(wp) :: zrhox 969 ! vertical interpolation: 970 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 971 REAL(wp), DIMENSION(1:jpk) :: h_out 972 INTEGER :: N_in, N_out, jref 973 REAL(wp) :: h_diff 974 LOGICAL :: northern_side,southern_side 975 !!--------------------------------------------- 836 976 ! 837 IF( before ) THEN 977 IF (before) THEN 978 IF(Agrif_UseSpecialValue) THEN 979 ! Agrif_SpecialValue = -999._wp 980 Agrif_SpecialValue = 0._wp 981 ELSE 982 Agrif_SpecialValue = 0._wp 983 ENDIF 984 DO jk=k1,k2 985 DO jj=j1,j2 986 DO ji=i1,i2 987 ptab(ji,jj,jk,1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) - & 988 & ((vmask(ji,jj,jk)-1) * Agrif_SpecialValue) 989 # if defined key_vertical 990 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 991 # endif 992 END DO 993 END DO 994 END DO 995 ELSE 996 zrhox = Agrif_rhox() 997 # if defined key_vertical 998 Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done 999 1000 southern_side = (nb == 2).AND.(ndir == 1) 1001 northern_side = (nb == 2).AND.(ndir == 2) 1002 1003 DO jj=j1,j2 1004 jref = jj 1005 IF (southern_side) jref = MAX(2,jj) 1006 IF (northern_side) jref = MIN(nlcj-2,jj) 1007 DO ji=i1,i2 1008 N_in = 0 1009 DO jk=k1,k2 1010 if (ptab(ji,jj,jk,2) == 0) EXIT 1011 N_in = N_in + 1 1012 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1013 h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1014 END DO 1015 IF (N_in == 0) THEN 1016 va(ji,jj,:) = 0._wp 1017 CYCLE 1018 ENDIF 1019 1020 N_out = 0 1021 DO jk=1,jpk 1022 if (vmask(ji,jref,jk) == 0) EXIT 1023 N_out = N_out + 1 1024 h_out(N_out) = e3v_a(ji,jref,jk) 1025 END DO 1026 IF (N_out == 0) THEN 1027 va(ji,jj,:) = 0._wp 1028 CYCLE 1029 ENDIF 1030 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),va(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 1031 END DO 1032 END DO 1033 # else 838 1034 DO jk = 1, jpkm1 839 ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 840 END DO 841 ELSE 842 zrhox= Agrif_Rhox() 843 DO jk = 1, jpkm1 844 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 845 END DO 1035 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_a(i1:i2,j1:j2,jk) ) 1036 END DO 1037 # endif 846 1038 ENDIF 847 1039 ! 848 1040 END SUBROUTINE interpvn 849 850 1041 851 1042 SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
Note: See TracChangeset
for help on using the changeset viewer.
