Changeset 8135
- Timestamp:
- 2017-06-05T11:01:20+02:00 (8 years ago)
- Location:
- branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 4 edited
-
agrif_oce.F90 (modified) (2 diffs)
-
agrif_opa_interp.F90 (modified) (9 diffs)
-
agrif_opa_update.F90 (modified) (6 diffs)
-
agrif_user.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_oce.F90
r5656 r8135 17 17 18 18 PUBLIC agrif_oce_alloc ! routine called by nemo_init in nemogcm.F90 19 PUBLIC reconstructandremap 19 20 20 21 ! !!* Namelist namagrif: AGRIF parameters … … 104 105 END FUNCTION agrif_oce_alloc 105 106 107 subroutine reconstructandremap(tabin,hin,tabout,hout,N,Nout) 108 implicit none 109 integer N, Nout 110 real tabin(N), tabout(Nout) 111 real hin(N), hout(Nout) 112 real coeffremap(N,3),zwork(N,3) 113 real zwork2(N+1,3) 114 integer k 115 double precision, parameter :: dsmll=1.0d-8 116 real q,q01,q02,q001,q002,q0 117 real z_win(1:N+1), z_wout(1:Nout+1) 118 real,parameter :: dpthin = 1.D-3 119 integer :: k1, kbox, ktop, ka, kbot 120 real :: tsum, qbot, rpsum, zbox, ztop, zthk, zbot, offset, qtop 121 122 z_win(1)=0.; z_wout(1)= 0. 123 do k=1,N 124 z_win(k+1)=z_win(k)+hin(k) 125 enddo 126 127 do k=1,Nout 128 z_wout(k+1)=z_wout(k)+hout(k) 129 enddo 130 131 do k=2,N 132 zwork(k,1)=1./(hin(k-1)+hin(k)) 133 enddo 134 135 do k=2,N-1 136 q0 = 1./(hin(k-1)+hin(k)+hin(k+1)) 137 zwork(k,2)=hin(k-1)+2.*hin(k)+hin(k+1) 138 zwork(k,3)=q0 139 enddo 140 141 do k= 2,N 142 zwork2(k,1)=zwork(k,1)*(tabin(k)-tabin(k-1)) 143 enddo 144 145 coeffremap(:,1) = tabin(:) 146 147 do k=2,N-1 148 q001 = hin(k)*zwork2(k+1,1) 149 q002 = hin(k)*zwork2(k,1) 150 if (q001*q002 < 0) then 151 q001 = 0. 152 q002 = 0. 153 endif 154 q=zwork(k,2) 155 q01=q*zwork2(k+1,1) 156 q02=q*zwork2(k,1) 157 if (abs(q001) > abs(q02)) q001 = q02 158 if (abs(q002) > abs(q01)) q002 = q01 159 160 q=(q001-q002)*zwork(k,3) 161 q001=q001-q*hin(k+1) 162 q002=q002+q*hin(k-1) 163 164 coeffremap(k,3)=coeffremap(k,1)+q001 165 coeffremap(k,2)=coeffremap(k,1)-q002 166 167 zwork2(k,1)=(2.*q001-q002)**2 168 zwork2(k,2)=(2.*q002-q001)**2 169 enddo 170 171 do k=1,N 172 if (k.eq.1 .or. k.eq.N .or. hin(k).le.dpthin) then 173 coeffremap(k,3) = coeffremap(k,1) 174 coeffremap(k,2) = coeffremap(k,1) 175 zwork2(k,1) = 0. 176 zwork2(k,2) = 0. 177 endif 178 enddo 179 180 do k=2,N 181 q002=max(zwork2(k-1,2),dsmll) 182 q001=max(zwork2(k,1),dsmll) 183 zwork2(k,3)=(q001*coeffremap(k-1,3)+q002*coeffremap(k,2))/(q001+q002) 184 enddo 185 186 zwork2(1,3) = 2*coeffremap(1,1)-zwork2(2,3) 187 zwork2(N+1,3)=2*coeffremap(N,1)-zwork2(N,3) 188 189 do k=1,N 190 q01=zwork2(k+1,3)-coeffremap(k,1) 191 q02=coeffremap(k,1)-zwork2(k,3) 192 q001=2.*q01 193 q002=2.*q02 194 if (q01*q02<0) then 195 q01=0. 196 q02=0. 197 elseif (abs(q01)>abs(q002)) then 198 q01=q002 199 elseif (abs(q02)>abs(q001)) then 200 q02=q001 201 endif 202 coeffremap(k,2)=coeffremap(k,1)-q02 203 coeffremap(k,3)=coeffremap(k,1)+q01 204 enddo 205 206 zbot=0.0 207 kbot=1 208 do k=1,Nout 209 ztop=zbot !top is bottom of previous layer 210 ktop=kbot 211 if (ztop.ge.z_win(ktop+1)) then 212 ktop=ktop+1 213 endif 214 215 zbot=z_wout(k+1) 216 zthk=zbot-ztop 217 218 if (zthk.gt.dpthin .and. ztop.lt.z_wout(Nout+1)) then 219 220 kbot=ktop 221 do while (z_win(kbot+1).lt.zbot.and.kbot.lt.N) 222 kbot=kbot+1 223 enddo 224 zbox=zbot 225 do k1= k+1,Nout 226 if (z_wout(k1+1)-z_wout(k1).gt.dpthin) then 227 exit !thick layer 228 else 229 zbox=z_wout(k1+1) !include thin adjacent layers 230 if (zbox.eq.z_wout(Nout+1)) then 231 exit !at bottom 232 endif 233 endif 234 enddo 235 zthk=zbox-ztop 236 237 kbox=ktop 238 do while (z_win(kbox+1).lt.zbox.and.kbox.lt.N) 239 kbox=kbox+1 240 enddo 241 242 if (ktop.eq.kbox) then 243 244 245 if (z_wout(k) .ne.z_win(kbox) .or.z_wout(k+1).ne.z_win(kbox+1) ) then 246 247 if (hin(kbox).gt.dpthin) then 248 q001 = (zbox-z_win(kbox))/hin(kbox) 249 q002 = (ztop-z_win(kbox))/hin(kbox) 250 q01=q001**2+q002**2+q001*q002+1.-2.*(q001+q002) 251 q02=q01-1.+(q001+q002) 252 q0=1.-q01-q02 253 else 254 q0 = 1.0 255 q01 = 0. 256 q02 = 0. 257 endif 258 tabout(k)=q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3) 259 260 else 261 tabout(k) = tabin(kbox) 262 263 endif 264 265 else 266 267 if (ktop.le.k .and. kbox.ge.k) then 268 ka = k 269 elseif (kbox-ktop.ge.3) then 270 ka = (kbox+ktop)/2 271 elseif (hin(ktop).ge.hin(kbox)) then 272 ka = ktop 273 else 274 ka = kbox 275 endif !choose ka 276 277 offset=coeffremap(ka,1) 278 279 qtop = z_win(ktop+1)-ztop !partial layer thickness 280 if (hin(ktop).gt.dpthin) then 281 q=(ztop-z_win(ktop))/hin(ktop) 282 q01=q*(q-1.) 283 q02=q01+q 284 q0=1-q01-q02 285 else 286 q0 = 1. 287 q01 = 0. 288 q02 = 0. 289 endif 290 291 tsum =((q0*coeffremap(ktop,1)+q01*coeffremap(ktop,2)+q02*coeffremap(ktop,3))-offset)*qtop 292 293 do k1= ktop+1,kbox-1 294 tsum =tsum +(coeffremap(k1,1)-offset)*hin(k1) 295 enddo !k1 296 297 298 qbot = zbox-z_win(kbox) !partial layer thickness 299 if (hin(kbox).gt.dpthin) then 300 q=qbot/hin(kbox) 301 q01=(q-1.)**2 302 q02=q01-1.+q 303 q0=1-q01-q02 304 else 305 q0 = 1.0 306 q01 = 0. 307 q02 = 0. 308 endif 309 310 tsum = tsum +((q0*coeffremap(kbox,1)+q01*coeffremap(kbox,2)+q02*coeffremap(kbox,3))-offset)*qbot 311 312 rpsum=1.0d0/zthk 313 tabout(k)=offset+tsum*rpsum 314 315 endif !single or multiple layers 316 else 317 if (k==1) then 318 write(*,'(a7,i4,i4,3f12.5)')'problem = ',N,Nout,zthk,z_wout(k+1),hout(1) 319 endif 320 tabout(k) = tabout(k-1) 321 322 endif !normal:thin layer 323 enddo !k 324 325 return 326 end subroutine reconstructandremap 327 106 328 #endif 107 329 !!====================================================================== -
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 -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r7646 r8135 158 158 INTEGER, INTENT(in) :: kt 159 159 ! 160 return 161 160 162 IF( ( Agrif_NbStepint() .NE. 0 ) .AND. (kt /= 0) ) RETURN 161 163 # if defined TWO_WAY … … 185 187 !! 186 188 INTEGER :: ji,jj,jk,jn 187 !!--------------------------------------------- 188 ! 189 IF (before) THEN 190 DO jn = n1,n2 189 ! VERTICAL REFINEMENT BEGIN 190 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: tabres_child 191 REAL(wp) :: h_in(k1:k2) 192 REAL(wp) :: h_out(1:jpk) 193 INTEGER :: N_in, N_out 194 REAL(wp) :: h_diff 195 REAL(wp) :: zrho_xy 196 REAL(wp) :: tabin(k1:k2,n1:n2) 197 ! VERTICAL REFINEMENT END 198 !!--------------------------------------------- 199 ! 200 IF (before) THEN 201 zrho_xy = Agrif_rhox() * Agrif_rhoy() 202 DO jn = n1,n2-1 191 203 DO jk=k1,k2 192 204 DO jj=j1,j2 193 205 DO ji=i1,i2 194 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn)206 tabres(ji,jj,jk,jn) = zrho_xy * tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 195 207 END DO 196 208 END DO 197 209 END DO 198 210 END DO 199 ELSE 211 DO jk=k1,k2 212 DO jj=j1,j2 213 DO ji=i1,i2 214 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * zrho_xy * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 215 END DO 216 END DO 217 END DO 218 219 ELSE 220 ! VERTICAL REFINEMENT BEGIN 221 tabres_child(:,:,:,:) = 0. 222 223 DO jj=j1,j2 224 DO ji=i1,i2 225 N_in = 0 226 DO jk=k1,k2 !k2 = jpk of child grid 227 IF (tabres(ji,jj,jk,n2) == 0) EXIT 228 N_in = N_in + 1 229 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 230 h_in(N_in) = tabres(ji,jj,jk,n2)/e1e2t(ji,jj) 231 ENDDO 232 N_out = 0 233 DO jk=1,jpk ! jpk of parent grid 234 IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF 235 N_out = N_out + 1 236 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 237 ENDDO 238 IF (N_in > 0) THEN 239 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 240 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 241 ! IF abs(h_diff > 1.e-8) THEN 242 ! N_in = N_in + 1 243 ! h_in(N_in) = h_diff 244 ! tabin(N_in,:) = tabin(N_in-1,:) 245 IF (h_diff < 0) THEN 246 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,sum(h_in(1:N_in)),sum(h_out(1:N_out)) 247 print *,'Nval = ',N_out,mbathy(ji,jj) 248 print *,'BATHY = ',gdepw_0(ji,jj,mbathy(ji,jj)+1),sum(e3t_0(ji,jj,1:mbathy(ji,jj))) 249 STOP 250 ! N_out = N_out + 1 251 ! h_out(N_out) = - h_diff 252 ENDIF 253 DO jn=n1,n2-1 254 CALL reconstructandremap(tabin(1:N_in,jn),h_in(1:N_in),tabres_child(ji,jj,1:N_out,jn),h_out(1:N_out),N_in,N_out) 255 ENDDO 256 ENDIF 257 ENDDO 258 ENDDO 259 260 ! WARNING : 261 ! tabres replaced by tabres_child in the following 262 ! k1 k2 replaced by 1 jpk 263 ! VERTICAL REFINEMENT END 264 200 265 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 201 266 ! Add asselin part 202 DO jn = n1,n2 203 DO jk= k1,k2267 DO jn = n1,n2-1 268 DO jk=1,jpk 204 269 DO jj=j1,j2 205 270 DO ji=i1,i2 206 IF( tabres (ji,jj,jk,jn) .NE. 0. ) THEN271 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 207 272 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 208 & + atfp * ( tabres (ji,jj,jk,jn) &273 & + atfp * ( tabres_child(ji,jj,jk,jn) & 209 274 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 210 275 ENDIF … … 214 279 ENDDO 215 280 ENDIF 216 DO jn = n1,n2 217 DO jk= k1,k2281 DO jn = n1,n2-1 282 DO jk=1,jpk 218 283 DO jj=j1,j2 219 284 DO ji=i1,i2 220 IF( tabres (ji,jj,jk,jn) .NE. 0. ) THEN221 tsn(ji,jj,jk,jn) = tabres (ji,jj,jk,jn) * tmask(ji,jj,jk)285 IF( tabres_child(ji,jj,jk,jn) .NE. 0. ) THEN 286 tsn(ji,jj,jk,jn) = tabres_child(ji,jj,jk,jn) * tmask(ji,jj,jk) 222 287 END IF 223 288 END DO … … 230 295 231 296 232 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, before )297 SUBROUTINE updateu( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 233 298 !!--------------------------------------------- 234 299 !! *** ROUTINE updateu *** 235 300 !!--------------------------------------------- 236 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2237 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres238 LOGICAL , INTENT(in ) :: before301 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 302 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2), INTENT(inout) :: tabres 303 LOGICAL , INTENT(in ) :: before 239 304 ! 240 305 INTEGER :: ji, jj, jk 241 306 REAL(wp) :: zrhoy 307 ! VERTICAL REFINEMENT BEGIN 308 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 309 REAL(wp) :: h_in(k1:k2) 310 REAL(wp) :: h_out(1:jpk) 311 INTEGER :: N_in, N_out 312 REAL(wp) :: h_diff 313 REAL(wp) :: tabin(k1:k2) 314 ! VERTICAL REFINEMENT END 242 315 !!--------------------------------------------- 243 316 ! 244 317 IF( before ) THEN 245 318 zrhoy = Agrif_Rhoy() 246 DO jk = k1, k2247 tabres(i1:i2,j1:j2,jk) = zrhoy * e2u(i1:i2,j1:j2) * e3u_n(i1:i2,j1:j2,jk) * un(i1:i2,j1:j2,jk)248 END DO249 ELSE250 319 DO jk=k1,k2 251 320 DO jj=j1,j2 252 321 DO ji=i1,i2 253 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 322 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) 323 tabres(ji,jj,jk,2) = umask(ji,jj,jk) * zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) 324 END DO 325 END DO 326 END DO 327 ELSE 328 ! VERTICAL REFINEMENT BEGIN 329 tabres_child(:,:,:) = 0. 330 331 DO jj=j1,j2 332 DO ji=i1,i2 333 N_in = 0 334 DO jk=k1,k2 !k2=jpk of child grid 335 IF (tabres(ji,jj,jk,2) == 0) EXIT 336 N_in = N_in + 1 337 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 338 h_in(N_in) = tabres(ji,jj,jk,2)/e2u(ji,jj) 339 ENDDO 340 N_out = 0 341 DO jk=1,jpk 342 IF (umask(ji,jj,jk) == 0) EXIT 343 N_out = N_out + 1 344 h_out(N_out) = e3u_n(ji,jj,jk) 345 ENDDO 346 IF (N_in * N_out > 0) THEN 347 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 348 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 349 if (h_diff < 0.) then 350 print *,'CHECK YOUR BATHY ...' 351 stop 352 ! else ! Extends with 0 353 ! N_in = N_in + 1 354 ! tabin(N_in) = 0. 355 ! h_in(N_in) = h_diff 356 endif 357 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 358 ENDIF 359 ENDDO 360 ENDDO 361 362 ! WARNING : 363 ! tabres replaced by tabres_child in the following 364 ! k1 k2 replaced by 1 jpk 365 ! remove division by fs e3u (already done) 366 ! VERTICAL REFINEMENT END 367 368 DO jk=1,jpk 369 DO jj=j1,j2 370 DO ji=i1,i2 371 !Following line now replaced by division higher up in vertical refinement case 372 ! tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk) 254 373 ! 255 374 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 256 375 ub(ji,jj,jk) = ub(ji,jj,jk) & 257 & + atfp * ( tabres (ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk)376 & + atfp * ( tabres_child(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 258 377 ENDIF 259 378 ! 260 un(ji,jj,jk) = tabres (ji,jj,jk) * umask(ji,jj,jk)379 un(ji,jj,jk) = tabres_child(ji,jj,jk) * umask(ji,jj,jk) 261 380 END DO 262 381 END DO … … 271 390 !! *** ROUTINE updatev *** 272 391 !!--------------------------------------------- 273 INTEGER :: i1,i2,j1,j2,k1,k2 392 INTEGER :: i1,i2,j1,j2,k1,k2,n1,n2 274 393 INTEGER :: ji,jj,jk 275 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ) :: tabres394 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,2) :: tabres 276 395 LOGICAL :: before 277 396 !! 278 397 REAL(wp) :: zrhox 398 ! VERTICAL REFINEMENT BEGIN 399 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 400 REAL(wp) :: h_in(k1:k2) 401 REAL(wp) :: h_out(1:jpk) 402 INTEGER :: N_in, N_out 403 REAL(wp) :: h_diff 404 REAL(wp) :: tabin(k1:k2) 405 ! VERTICAL REFINEMENT END 279 406 !!--------------------------------------------- 280 407 ! … … 284 411 DO jj=j1,j2 285 412 DO ji=i1,i2 286 tabres(ji,jj,jk) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk) 287 END DO 288 END DO 289 END DO 290 ELSE 291 DO jk=k1,k2 413 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 414 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) 415 END DO 416 END DO 417 END DO 418 ELSE 419 ! VERTICAL REFINEMENT BEGIN 420 tabres_child(:,:,:) = 0. 421 422 DO jj=j1,j2 423 DO ji=i1,i2 424 N_in = 0 425 DO jk=k1,k2 426 IF (tabres(ji,jj,jk,2) == 0) EXIT 427 N_in = N_in + 1 428 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) 429 h_in(N_in) = tabres(ji,jj,jk,2)/e1v(ji,jj) 430 ENDDO 431 N_out = 0 432 DO jk=1,jpk 433 IF (vmask(ji,jj,jk) == 0) EXIT 434 N_out = N_out + 1 435 h_out(N_out) = e3v_n(ji,jj,jk) 436 ENDDO 437 IF (N_in * N_out > 0) THEN 438 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 439 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 440 if (h_diff < 0.) then 441 print *,'CHECK YOUR BATHY ...' 442 stop 443 ! else ! Extends with 0 444 ! N_in = N_in + 1 445 ! tabin(N_in) = 0. 446 ! h_in(N_in) = h_diff 447 endif 448 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out) 449 ENDIF 450 ENDDO 451 ENDDO 452 453 ! WARNING : 454 ! tabres replaced by tabres_child in the following 455 ! k1 k2 replaced by 1 jpk 456 ! remove division by fs e3v (already done) 457 ! VERTICAL REFINEMENT END 458 459 DO jk=k1,jpk 292 460 DO jj=j1,j2 293 461 DO ji=i1,i2 294 tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 462 !Following line now replaced by division higher up in vertical refinement case 463 ! tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk) 295 464 ! 296 465 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 297 466 vb(ji,jj,jk) = vb(ji,jj,jk) & 298 & + atfp * ( tabres (ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk)467 & + atfp * ( tabres_child(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 299 468 ENDIF 300 469 ! 301 vn(ji,jj,jk) = tabres (ji,jj,jk) * vmask(ji,jj,jk)470 vn(ji,jj,jk) = tabres_child(ji,jj,jk) * vmask(ji,jj,jk) 302 471 END DO 303 472 END DO -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r7761 r8135 346 346 ! 1. Declaration of the type of variable which have to be interpolated 347 347 !--------------------------------------------------------------------- 348 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts /),tsn_id)348 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 349 349 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts/),tsn_sponge_id) 350 350 351 CALL agrif_declare_variable((/1,2,0 /),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_interp_id)352 CALL agrif_declare_variable((/2,1,0 /),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_interp_id)353 CALL agrif_declare_variable((/1,2,0 /),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_update_id)354 CALL agrif_declare_variable((/2,1,0 /),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_update_id)351 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) 352 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_interp_id) 353 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 354 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 355 355 CALL agrif_declare_variable((/1,2,0/),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id) 356 356 CALL agrif_declare_variable((/2,1,0/),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id)
Note: See TracChangeset
for help on using the changeset viewer.
