Changeset 6258 for branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
- Timestamp:
- 2016-01-15T13:11:56+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90
r5819 r6258 38 38 39 39 PUBLIC Agrif_tra, Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts 40 ! VERTICAL REFINEMENT BEGIN 41 PUBLIC Agrif_Init_InterpScales 42 ! VERTICAL REFINEMENT END 40 43 PUBLIC interpun, interpvn, interpun2d, interpvn2d 41 44 PUBLIC interptsn, interpsshn … … 50 53 !!---------------------------------------------------------------------- 51 54 !! NEMO/NST 3.6 , NEMO Consortium (2010) 52 !! $Id $55 !! $Id: agrif_opa_interp.F90 5081 2015-02-13 09:51:27Z smasson $ 53 56 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 54 57 !!---------------------------------------------------------------------- 55 58 59 ! VERTICAL REFINEMENT BEGIN 60 REAL, DIMENSION(:,:,:), ALLOCATABLE :: interp_scales_t, interp_scales_u, interp_scales_v 61 !$AGRIF_DO_NOT_TREAT 62 LOGICAL :: scaleT, scaleU, scaleV = .FALSE. 63 !$AGRIF_END_DO_NOT_TREAT 64 ! VERTICAL REFINEMENT END 65 56 66 CONTAINS 67 68 ! VERTICAL REFINEMENT BEGIN 69 70 SUBROUTINE Agrif_Init_InterpScales 71 72 scaleT = .TRUE. 73 Call Agrif_Bc_Variable(scales_t_id,calledweight=1.,procname=interpscales) 74 scaleT = .FALSE. 75 76 scaleU = .TRUE. 77 Call Agrif_Bc_Variable(scales_u_id,calledweight=1.,procname=interpscales) 78 scaleU = .FALSE. 79 80 scaleV = .TRUE. 81 Call Agrif_Bc_Variable(scales_v_id,calledweight=1.,procname=interpscales) 82 scaleV = .FALSE. 83 84 END SUBROUTINE Agrif_Init_InterpScales 85 86 SUBROUTINE interpscales(ptab,i1,i2,j1,j2,k1,k2,before) 87 !!--------------------------------------------- 88 !! *** ROUTINE interpscales *** 89 !!--------------------------------------------- 90 91 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 92 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 93 94 INTEGER :: ji, jj, jk 95 LOGICAL :: before 96 97 IF (before) THEN 98 IF (scaleT ) THEN 99 DO jk=k1,k2 100 DO jj=j1,j2 101 DO ji=i1,i2 102 ! ptab(ji,jj,jk) = fse3t_n(ji,jj,jk) * tmask(ji,jj,jk) 103 ptab(ji,jj,jk) = fse3t_n(ji,jj,jk) 104 END DO 105 END DO 106 END DO 107 ELSEIF (scaleU) THEN 108 DO jk=k1,k2 109 DO jj=j1,j2 110 DO ji=i1,i2 111 ! ptab(ji,jj,jk) = fse3u_n(ji,jj,jk) * umask(ji,jj,jk) 112 ! ptab(ji,jj,jk) = fse3u_n(ji,jj,jk) 113 ptab(ji,jj,jk) = umask(ji,jj,jk) 114 END DO 115 END DO 116 END DO 117 ELSEIF (scaleV) THEN 118 DO jk=k1,k2 119 DO jj=j1,j2 120 DO ji=i1,i2 121 ! ptab(ji,jj,jk) = fse3v_n(ji,jj,jk) * vmask(ji,jj,jk) 122 ! ptab(ji,jj,jk) = fse3v_n(ji,jj,jk) 123 ptab(ji,jj,jk) = vmask(ji,jj,jk) 124 END DO 125 END DO 126 END DO 127 ENDIF 128 ELSE 129 IF (scaleT ) THEN 130 IF (.not.allocated(interp_scales_t)) allocate(interp_scales_t(jpi,jpj,k1:k2)) 131 DO jk=k1,k2 132 DO jj=j1,j2 133 DO ji=i1,i2 134 interp_scales_t(ji,jj,jk) = ptab(ji,jj,jk) 135 END DO 136 END DO 137 END DO 138 ELSEIF (scaleU) THEN 139 IF (.not.allocated(interp_scales_u)) allocate(interp_scales_u(jpi,jpj,k1:k2)) 140 DO jk=k1,k2 141 DO jj=j1,j2 142 DO ji=i1,i2 143 interp_scales_u(ji,jj,jk) = ptab(ji,jj,jk) 144 END DO 145 END DO 146 END DO 147 ELSEIF (scaleV) THEN 148 IF (.not.allocated(interp_scales_v)) allocate(interp_scales_v(jpi,jpj,k1:k2)) 149 DO jk=k1,k2 150 DO jj=j1,j2 151 DO ji=i1,i2 152 interp_scales_v(ji,jj,jk) = ptab(ji,jj,jk) 153 END DO 154 END DO 155 END DO 156 ENDIF 157 ENDIF 158 159 END SUBROUTINE interpscales 160 161 ! VERTICAL REFINEMENT END 57 162 58 163 SUBROUTINE Agrif_tra … … 611 716 REAL(wp) :: zalpha 612 717 ! 718 return 719 613 720 zalpha = REAL( Agrif_NbStepint() + Agrif_IRhot() - 1, wp ) / REAL( Agrif_IRhot(), wp ) 614 721 IF( zalpha > 1. ) zalpha = 1. … … 638 745 REAL(wp) :: zalpha4, zalpha5, zalpha6, zalpha7 639 746 LOGICAL :: western_side, eastern_side,northern_side,southern_side 747 ! VERTICAL REFINEMENT BEGIN 748 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child 749 REAL(wp) :: h_in(k1:k2) 750 REAL(wp) :: h_out(1:jpk) 751 INTEGER :: N_in, N_out 752 REAL(wp) :: h_diff 753 ! VERTICAL REFINEMENT END 640 754 641 755 IF (before) THEN 642 756 ptab(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 643 ELSE 757 ELSE 758 ! VERTICAL REFINEMENT BEGIN 759 760 ptab_child(:,:,:,:) = 0. 761 do jj=j1,j2 762 do ji=i1,i2 763 h_in(k1:k2) = interp_scales_t(ji,jj,k1:k2) 764 h_out(1:jpk) = fse3t(ji,jj,1:jpk) 765 h_diff = sum(h_out(1:jpk-1))-sum(h_in(k1:k2-1)) 766 N_in = k2-1 767 N_out = jpk-1 768 if (h_diff > 0) then 769 h_in(N_in+1) = h_diff 770 N_in = N_in + 1 771 else 772 h_out(N_out+1) = -h_diff 773 N_out = N_out + 1 774 endif 775 ptab(ji,jj,k2,:) = ptab(ji,jj,k2-1,:) 776 do jn=1,jpts 777 call reconstructandremap(ptab(ji,jj,1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 778 enddo 779 ! if (abs(h_diff) > 1000.) then 780 ! do jn=1,jpts 781 ! do jk=1,N_out 782 ! print *,'AVANT APRES = ',ji,jj,jk,N_out,ptab(ji,jj,jk,jn),ptab_child(ji,jj,jk,jn) 783 ! enddo 784 ! enddo 785 ! endif 786 enddo 787 enddo 788 789 ! VERTICAL REFINEMENT END 790 644 791 ! 645 792 western_side = (nb == 1).AND.(ndir == 1) … … 671 818 IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2 672 819 ! 820 ! VERTICAL REFINEMENT BEGIN 821 822 ! WARNING : 823 ! ptab replaced by ptab_child in the following 824 ! k1 k2 replaced by 1 jpk 825 ! VERTICAL REFINEMENT END 673 826 IF( eastern_side) THEN 674 827 DO jn = 1, jpts 675 tsa(nlci,j1:j2, k1:k2,jn) = zalpha1 * ptab(nlci,j1:j2,k1:k2,jn) + zalpha2 * ptab(nlci-1,j1:j2,k1:k2,jn)828 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) 676 829 DO jk = 1, jpkm1 677 830 DO jj = jmin,jmax … … 692 845 IF( northern_side ) THEN 693 846 DO jn = 1, jpts 694 tsa(i1:i2,nlcj, k1:k2,jn) = zalpha1 * ptab(i1:i2,nlcj,k1:k2,jn) + zalpha2 * ptab(i1:i2,nlcj-1,k1:k2,jn)847 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) 695 848 DO jk = 1, jpkm1 696 849 DO ji = imin,imax … … 711 864 IF( western_side) THEN 712 865 DO jn = 1, jpts 713 tsa(1,j1:j2, k1:k2,jn) = zalpha1 * ptab(1,j1:j2,k1:k2,jn) + zalpha2 * ptab(2,j1:j2,k1:k2,jn)866 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) 714 867 DO jk = 1, jpkm1 715 868 DO jj = jmin,jmax … … 729 882 IF( southern_side ) THEN 730 883 DO jn = 1, jpts 731 tsa(i1:i2,1, k1:k2,jn) = zalpha1 * ptab(i1:i2,1,k1:k2,jn) + zalpha2 * ptab(i1:i2,2,k1:k2,jn)884 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) 732 885 DO jk=1,jpk 733 886 DO ji=imin,imax … … 749 902 ! East south 750 903 IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 751 tsa(nlci-1,2,:,:) = ptab (nlci-1,2,:,:)904 tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,:) 752 905 ENDIF 753 906 ! East north 754 907 IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 755 tsa(nlci-1,nlcj-1,:,:) = ptab (nlci-1,nlcj-1,:,:)908 tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,:) 756 909 ENDIF 757 910 ! West south 758 911 IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN 759 tsa(2,2,:,:) = ptab (2,2,:,:)912 tsa(2,2,:,:) = ptab_child(2,2,:,:) 760 913 ENDIF 761 914 ! West north 762 915 IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN 763 tsa(2,nlcj-1,:,:) = ptab (2,nlcj-1,:,:)916 tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,:) 764 917 ENDIF 765 918 ! … … 794 947 END SUBROUTINE interpsshn 795 948 796 SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2, before)949 SUBROUTINE interpun(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir) 797 950 !!--------------------------------------------- 798 951 !! *** ROUTINE interpun *** 799 952 !!--------------------------------------------- 800 953 !! 801 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 802 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: ptab954 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 955 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 803 956 LOGICAL, INTENT(in) :: before 957 INTEGER, INTENT(in) :: nb , ndir 804 958 !! 805 959 INTEGER :: ji,jj,jk 806 REAL(wp) :: zrhoy 960 REAL(wp) :: zrhoy 961 ! VERTICAL REFINEMENT BEGIN 962 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 963 REAL(wp), DIMENSION(k1:k2) :: tabin 964 REAL(wp) :: h_in(k1:k2) 965 REAL(wp) :: h_out(1:jpk) 966 INTEGER :: N_in, N_out 967 REAL(wp) :: h_diff 968 LOGICAL :: western_side, eastern_side 969 INTEGER :: iref 970 971 ! VERTICAL REFINEMENT END 807 972 !!--------------------------------------------- 808 973 ! … … 811 976 DO jj=j1,j2 812 977 DO ji=i1,i2 813 ptab(ji,jj,jk) = e2u(ji,jj) * un(ji,jj,jk) 814 ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3u(ji,jj,jk) 978 ptab(ji,jj,jk,1) = e2u(ji,jj) * un(ji,jj,jk) 979 ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3u(ji,jj,jk) 980 ptab(ji,jj,jk,2) = fse3u(ji,jj,jk) 815 981 END DO 816 982 END DO 817 983 END DO 818 984 ELSE 985 ! VERTICAL REFINEMENT BEGIN 986 western_side = (nb == 1).AND.(ndir == 1) 987 eastern_side = (nb == 1).AND.(ndir == 2) 988 989 ptab_child(:,:,:) = 0. 990 do jj=j1,j2 991 do ji=i1,i2 992 iref = ji 993 IF (western_side) iref = 2 994 IF (eastern_side) iref = nlci-2 995 996 N_in = 0 997 do jk=k1,k2 998 if (interp_scales_u(ji,jj,jk) == 0) EXIT 999 N_in = N_in + 1 1000 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1001 h_in(N_in) = ptab(ji,jj,jk,2) 1002 enddo 1003 1004 IF (N_in == 0) THEN 1005 ptab_child(ji,jj,:) = 0. 1006 CYCLE 1007 ENDIF 1008 1009 N_out = 0 1010 do jk=1,jpk 1011 if (umask(iref,jj,jk) == 0) EXIT 1012 N_out = N_out + 1 1013 h_out(N_out) = fse3u(ji,jj,jk) 1014 enddo 1015 1016 IF (N_out == 0) THEN 1017 ptab_child(ji,jj,:) = 0. 1018 CYCLE 1019 ENDIF 1020 1021 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 1022 IF (h_diff > 0.) THEN 1023 N_in = N_in + 1 1024 h_in(N_in) = h_diff 1025 tabin(N_in) = 0. 1026 ELSE 1027 h_out(N_out) = h_out(N_out) - h_diff 1028 ENDIF 1029 1030 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) 1031 1032 ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3u(ji,jj,N_out) 1033 1034 ENDDO 1035 ENDDO 1036 1037 ! in the following 1038 ! remove division of ua by fs e3u (already done) 1039 ! VERTICAL REFINEMENT END 1040 819 1041 zrhoy = Agrif_Rhoy() 820 1042 DO jk=1,jpkm1 821 1043 DO jj=j1,j2 822 ua(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj))) 823 ua(i1:i2,jj,jk) = ua(i1:i2,jj,jk) / fse3u(i1:i2,jj,jk) 1044 ua(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhoy*e2u(i1:i2,jj))) 824 1045 END DO 825 1046 END DO … … 861 1082 862 1083 863 SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2, before)1084 SUBROUTINE interpvn(ptab,i1,i2,j1,j2,k1,k2,m1,m2,before,nb,ndir) 864 1085 !!--------------------------------------------- 865 1086 !! *** ROUTINE interpvn *** 866 1087 !!--------------------------------------------- 867 1088 ! 868 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 869 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: ptab1089 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 1090 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab 870 1091 LOGICAL, INTENT(in) :: before 1092 INTEGER, INTENT(in) :: nb , ndir 871 1093 ! 872 1094 INTEGER :: ji,jj,jk 873 REAL(wp) :: zrhox 1095 REAL(wp) :: zrhox 1096 ! VERTICAL REFINEMENT BEGIN 1097 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child 1098 REAL(wp), DIMENSION(k1:k2) :: tabin 1099 REAL(wp) :: h_in(k1:k2) 1100 REAL(wp) :: h_out(1:jpk) 1101 INTEGER :: N_in, N_out 1102 REAL(wp) :: h_diff 1103 LOGICAL :: northern_side,southern_side 1104 INTEGER :: jref 1105 1106 ! VERTICAL REFINEMENT END 874 1107 !!--------------------------------------------- 875 1108 ! … … 879 1112 DO jj=j1,j2 880 1113 DO ji=i1,i2 881 ptab(ji,jj,jk) = e1v(ji,jj) * vn(ji,jj,jk) 882 ptab(ji,jj,jk) = ptab(ji,jj,jk) * fse3v(ji,jj,jk) 1114 ptab(ji,jj,jk,1) = e1v(ji,jj) * vn(ji,jj,jk) 1115 ptab(ji,jj,jk,1) = ptab(ji,jj,jk,1) * fse3v(ji,jj,jk) 1116 ptab(ji,jj,jk,2) = fse3v(ji,jj,jk) 883 1117 END DO 884 1118 END DO 885 1119 END DO 886 ELSE 1120 ELSE 1121 ! VERTICAL REFINEMENT BEGIN 1122 ptab_child(:,:,:) = 0. 1123 southern_side = (nb == 2).AND.(ndir == 1) 1124 northern_side = (nb == 2).AND.(ndir == 2) 1125 do jj=j1,j2 1126 jref = jj 1127 IF (southern_side) jref = 2 1128 IF (northern_side) jref = nlcj-2 1129 do ji=i1,i2 1130 1131 N_in = 0 1132 do jk=k1,k2 1133 if (interp_scales_v(ji,jj,jk) == 0) EXIT 1134 N_in = N_in + 1 1135 tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2) 1136 h_in(N_in) = ptab(ji,jj,jk,2) 1137 enddo 1138 IF (N_in == 0) THEN 1139 ptab_child(ji,jj,:) = 0. 1140 CYCLE 1141 ENDIF 1142 1143 N_out = 0 1144 do jk=1,jpk 1145 if (vmask(ji,jref,jk) == 0) EXIT 1146 N_out = N_out + 1 1147 h_out(N_out) = fse3v(ji,jj,jk) 1148 enddo 1149 IF (N_out == 0) THEN 1150 ptab_child(ji,jj,:) = 0. 1151 CYCLE 1152 ENDIF 1153 1154 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 1155 IF (h_diff > 0.) THEN 1156 N_in = N_in + 1 1157 h_in(N_in) = h_diff 1158 tabin(N_in) = 0. 1159 ELSE 1160 h_out(N_out) = h_out(N_out) - h_diff 1161 ENDIF 1162 1163 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) 1164 1165 ptab_child(ji,jj,N_out) = ptab_child(ji,jj,N_out) * h_out(N_out) / fse3v(ji,jj,N_out) 1166 1167 enddo 1168 enddo 1169 ! in the following 1170 ! remove division of va by fs e3v (already done) 1171 ! VERTICAL REFINEMENT END 887 1172 zrhox= Agrif_Rhox() 888 1173 DO jk=1,jpkm1 889 1174 DO jj=j1,j2 890 va(i1:i2,jj,jk) = (ptab(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj))) 891 va(i1:i2,jj,jk) = va(i1:i2,jj,jk) / fse3v(i1:i2,jj,jk) 1175 va(i1:i2,jj,jk) = (ptab_child(i1:i2,jj,jk)/(zrhox*e1v(i1:i2,jj))) 892 1176 END DO 893 1177 END DO
Note: See TracChangeset
for help on using the changeset viewer.