Changeset 8135 for branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
- Timestamp:
- 2017-06-05T11:01:20+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r7646 r8135 576 576 !!---------------------------------------------------------------------- 577 577 ! 578 return 579 578 580 zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 579 581 IF( zalpha > 1. ) zalpha = 1. … … 603 605 REAL(wp) :: zrhox , zalpha1, zalpha2, zalpha3 604 606 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 605 LOGICAL :: western_side, eastern_side,northern_side,southern_side 606 !!---------------------------------------------------------------------- 607 ! 607 LOGICAL :: western_side, eastern_side,northern_side,southern_side 608 ! VERTICAL REFINEMENT BEGIN 609 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 610 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 611 REAL(wp) :: h_in(k1:k2) 612 REAL(wp) :: h_out(1:jpk) 613 INTEGER :: N_in, N_out 614 REAL(wp) :: h_diff 615 REAL(wp) :: zrhoxy 616 ! VERTICAL REFINEMENT END 617 618 zrhoxy = Agrif_rhox()*Agrif_rhoy() 608 619 IF (before) THEN 609 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 610 ELSE 620 DO jn = n1,n2-1 621 DO jk=k1,k2 622 DO jj=j1,j2 623 DO ji=i1,i2 624 ptab(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 625 END DO 626 END DO 627 END DO 628 END DO 629 DO jk=k1,k2 630 DO jj=j1,j2 631 DO ji=i1,i2 632 ptab(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 633 END DO 634 END DO 635 END DO 636 637 ELSE 638 ! VERTICAL REFINEMENT BEGIN 639 640 ptab_child(:,:,:,:) = 0. 641 do jj=j1,j2 642 do ji=i1,i2 643 N_in = 0 644 DO jk=k1,k2 !k2 = jpk of parent grid 645 IF (ptab(ji,jj,jk,n2) == 0) EXIT 646 N_in = N_in + 1 647 tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)/ptab(ji,jj,jk,n2) 648 h_in(N_in) = ptab(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 649 END DO 650 N_out = 0 651 DO jk=1,jpk ! jpk of child grid 652 IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 653 N_out = N_out + 1 654 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 655 ENDDO 656 IF (N_in > 0) THEN 657 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 658 ! if (h_diff > 0) then 659 ! h_in(N_in+1) = h_diff 660 ! N_in = N_in + 1 661 ! else 662 ! h_out(N_out+1) = -h_diff 663 ! N_out = N_out + 1 664 ! endif 665 ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) !what is this line for????? 666 do jn=1,jpts 667 call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 668 enddo 669 ! if (abs(h_diff) > 1000.) then 670 ! do jn=1,jpts 671 ! do jk=1,N_out 672 ! print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 673 ! enddo 674 ! enddo 675 ! endif 676 ENDIF 677 enddo 678 enddo 679 680 681 ! VERTICAL REFINEMENT END 682 611 683 ! 612 684 western_side = (nb == 1).AND.(ndir == 1) … … 640 712 IF( eastern_side ) THEN 641 713 DO jn = 1, jpts 642 tsa(nlci,j1:j2, k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)714 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) 643 715 DO jk = 1, jpkm1 644 716 DO jj = jmin,jmax … … 660 732 IF( northern_side ) THEN 661 733 DO jn = 1, jpts 662 tsa(i1:i2,nlcj, k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)734 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) 663 735 DO jk = 1, jpkm1 664 736 DO ji = imin,imax … … 680 752 IF( western_side ) THEN 681 753 DO jn = 1, jpts 682 tsa(1,j1:j2, k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)754 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) 683 755 DO jk = 1, jpkm1 684 756 DO jj = jmin,jmax … … 699 771 IF( southern_side ) THEN 700 772 DO jn = 1, jpts 701 tsa(i1:i2,1, k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)773 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) 702 774 DO jk = 1, jpk 703 775 DO ji=imin,imax … … 720 792 ! East south 721 793 IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 722 tsa(nlci-1,2,:,:) = ptab (nlci-1,2,:,:)794 tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts) 723 795 ENDIF 724 796 ! East north 725 797 IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 726 tsa(nlci-1,nlcj-1,:,:) = ptab (nlci-1,nlcj-1,:,:)798 tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts) 727 799 ENDIF 728 800 ! West south 729 801 IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 730 tsa(2,2,:,:) = ptab (2,2,:,:)802 tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts) 731 803 ENDIF 732 804 ! West north 733 805 IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 734 tsa(2,nlcj-1,:,:) = ptab (2,nlcj-1,:,:)806 tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts) 735 807 ENDIF 736 808 ! … … 771 843 !!---------------------------------------------------------------------- 772 844 !! *** ROUTINE interpun *** 773 !!---------------------------------------------------------------------- 774 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 775 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 776 LOGICAL , INTENT(in ) :: before 777 ! 778 INTEGER :: ji, jj, jk 779 REAL(wp) :: zrhoy 780 !!---------------------------------------------------------------------- 781 ! 782 IF( before ) THEN 783 DO jk = k1, jpk 784 ptab(i1:i2,j1:j2,jk) = e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk) 845 !!--------------------------------------------- 846 !! 847 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 848 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 849 LOGICAL, INTENT(in) :: before 850 INTEGER, INTENT(in) :: nb , ndir 851 !! 852 INTEGER :: ji,jj,jk 853 REAL(wp) :: zrhoy 854 ! VERTICAL REFINEMENT BEGIN 855 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 856 REAL(wp), DIMENSION(k1:k2) :: tabin 857 REAL(wp) :: h_in(k1:k2) 858 REAL(wp) :: h_out(1:jpk) 859 INTEGER :: N_in, N_out 860 REAL(wp) :: h_diff 861 LOGICAL :: western_side, eastern_side 862 INTEGER :: iref 863 864 ! VERTICAL REFINEMENT END 865 !!--------------------------------------------- 866 ! 867 zrhoy = Agrif_rhoy() 868 IF (before) THEN 869 !We can't use zero as the special value because we need to include zeros 870 !when interpolating the scale factors 871 IF(Agrif_UseSpecialValue) THEN 872 Agrif_SpecialValue = -999._wp 873 ELSE 874 Agrif_SpecialValue = 0._wp 875 ENDIF 876 DO jk=1,jpk 877 DO jj=j1,j2 878 DO ji=i1,i2 879 ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - & 880 & ((umask(ji,jj,jk)-1) * Agrif_SpecialValue) 881 ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 882 END DO 883 END DO 785 884 END DO 786 885 ELSE 787 zrhoy = Agrif_Rhoy() 886 ! VERTICAL REFINEMENT BEGIN 887 western_side = (nb == 1).AND.(ndir == 1) 888 eastern_side = (nb == 1).AND.(ndir == 2) 889 890 Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed 891 892 ptab_child(:,:,:) = 0. 893 DO jj=j1,j2 894 DO ji=i1,i2 895 iref = ji 896 IF (western_side) iref = 2 897 IF (eastern_side) iref = nlci-2 898 899 N_in = 0 900 DO jk=k1,k2 901 IF (ptab(ji,jj,jk,2) == 0) EXIT 902 N_in = N_in + 1 903 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 904 h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj) * zrhoy) 905 ENDDO 906 907 IF (N_in == 0) THEN 908 ptab_child(ji,jj,:) = 0. 909 CYCLE 910 ENDIF 911 912 N_out = 0 913 DO jk=1,jpk 914 if (umask(ji,jj,jk) == 0) EXIT 915 N_out = N_out + 1 916 h_out(N_out) = e3u_n(ji,jj,jk) 917 ENDDO 918 919 IF (N_out == 0) THEN 920 ptab_child(ji,jj,:) = 0. 921 CYCLE 922 ENDIF 923 924 ! IF (N_in * N_out > 0) THEN 925 ! h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 926 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 927 ! if (h_diff < 0.) then 928 ! print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 929 ! stop 930 ! endif 931 ! ENDIF 932 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 933 934 ENDDO 935 ENDDO 936 937 ! in the following 938 ! remove division of ua by fs e3u (already done) and also zrhoy and e2u 939 ! VERTICAL REFINEMENT END 940 788 941 DO jk = 1, jpkm1 789 942 DO jj=j1,j2 790 ua(i1:i2,jj,jk) = ptab(i1:i2,jj,jk) / ( zrhoy * e2u(i1:i2,jj) * e3u_n(i1:i2,jj,jk) ) 943 ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk) 944 !/(zrhoy*e2u(i1:i2,jj))) 791 945 END DO 792 946 END DO … … 800 954 !! *** ROUTINE interpvn *** 801 955 !!---------------------------------------------------------------------- 802 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2 803 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 804 LOGICAL , INTENT(in ) :: before 805 ! 806 INTEGER :: ji, jj, jk 807 REAL(wp) :: zrhox 808 !!---------------------------------------------------------------------- 956 ! 957 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 958 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 959 LOGICAL, INTENT(in) :: before 960 INTEGER, INTENT(in) :: nb , ndir 961 ! 962 INTEGER :: ji,jj,jk 963 REAL(wp) :: zrhox 964 ! VERTICAL REFINEMENT BEGIN 965 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 966 REAL(wp), DIMENSION(k1:k2) :: tabin 967 REAL(wp) :: h_in(k1:k2) 968 REAL(wp) :: h_out(1:jpk) 969 INTEGER :: N_in, N_out 970 REAL(wp) :: h_diff 971 LOGICAL :: northern_side,southern_side 972 INTEGER :: jref 973 974 ! VERTICAL REFINEMENT END 975 !!--------------------------------------------- 809 976 ! 810 IF( before ) THEN !interpv entre 1 et k2 et interpv2d en jpkp1 811 DO jk = k1, jpk 812 ptab(i1:i2,j1:j2,jk) = e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) * vn(i1:i2,j1:j2,jk) 813 END DO 814 ELSE 815 zrhox= Agrif_Rhox() 816 DO jk = 1, jpkm1 817 va(i1:i2,j1:j2,jk) = ptab(i1:i2,j1:j2,jk) / ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) ) 977 zrhox = Agrif_rhox() 978 IF (before) THEN 979 IF(Agrif_UseSpecialValue) THEN 980 Agrif_SpecialValue = -999._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) 988 ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk) 989 END DO 990 END DO 991 END DO 992 ELSE 993 ! VERTICAL REFINEMENT BEGIN 994 ptab_child(:,:,:) = 0. 995 southern_side = (nb == 2).AND.(ndir == 1) 996 northern_side = (nb == 2).AND.(ndir == 2) 997 998 Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done 999 1000 do jj=j1,j2 1001 jref = jj 1002 IF (southern_side) jref = 2 1003 IF (northern_side) jref = nlcj-2 1004 do ji=i1,i2 1005 1006 N_in = 0 1007 do jk=k1,k2 1008 if (ptab(ji,jj,jk,2) == 0) EXIT 1009 N_in = N_in + 1 1010 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1011 h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox) 1012 enddo 1013 IF (N_in == 0) THEN 1014 ptab_child(ji,jj,:) = 0. 1015 CYCLE 1016 ENDIF 1017 1018 N_out = 0 1019 do jk=1,jpk 1020 if (vmask(ji,jref,jk) == 0) EXIT 1021 N_out = N_out + 1 1022 h_out(N_out) = e3v_n(ji,jj,jk) 1023 enddo 1024 IF (N_out == 0) THEN 1025 ptab_child(ji,jj,:) = 0. 1026 CYCLE 1027 ENDIF 1028 1029 ! IF (N_in * N_out > 0) THEN 1030 ! h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 1031 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 1032 ! if (h_diff < 0.) then 1033 ! print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 1034 ! stop 1035 ! endif 1036 ! ENDIF 1037 call reconstructandremap(tabin(1:N_in),h_in(1:N_in),ptab_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 1038 1039 enddo 1040 enddo 1041 ! in the following 1042 ! remove division of va by fs e3v, zrhox and e1v (already done) 1043 ! VERTICAL REFINEMENT END 1044 DO jk=1,jpkm1 1045 DO jj=j1,j2 1046 va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk) 1047 END DO 818 1048 END DO 819 1049 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
