Changeset 10087 for vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90
- Timestamp:
- 2018-09-05T15:33:44+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90
r5656 r10087 41 41 integer, dimension(:), allocatable :: indparentppm_1d, indchildppm_1d 42 42 ! 43 43 44 private :: Agrif_limiter_vanleer 44 45 ! … … 56 57 integer, intent(in) :: np !< Length of input array 57 58 integer, intent(in) :: nc !< Length of output array 58 real , intent(in) :: s_parent !< Parent grid position (s_root = 0)59 real , intent(in) :: s_child !< Child grid position (s_root = 0)60 real , intent(in) :: ds_parent !< Parent grid dx (ds_root = 1)61 real , intent(in) :: ds_child !< Child grid dx (ds_root = 1)59 real(kind=8), intent(in) :: s_parent !< Parent grid position (s_root = 0) 60 real(kind=8), intent(in) :: s_child !< Child grid position (s_root = 0) 61 real(kind=8), intent(in) :: ds_parent !< Parent grid dx (ds_root = 1) 62 real(kind=8), intent(in) :: ds_child !< Child grid dx (ds_root = 1) 62 63 ! 63 64 integer :: i, coeffraf, locind_parent_left 64 real :: globind_parent_left, globind_parent_right65 real :: invds, invds2, ypos, ypos2, diff65 real(kind=8) :: globind_parent_left, globind_parent_right 66 real(kind=8) :: invds, invds2, ypos, ypos2, diff 66 67 ! 67 68 coeffraf = nint(ds_parent/ds_child) … … 92 93 ! 93 94 diff = globind_parent_right - ypos2 95 ! quick fix for roundoff error 96 diff=nint(diff*coeffraf)/real(coeffraf) 97 94 98 y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1)) 99 95 100 ypos2 = ypos2 + invds2 96 101 ! … … 104 109 else 105 110 globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent 106 y(nc) = ((globind_parent_left + ds_parent - ypos)*x(locind_parent_left) & 107 + (ypos - globind_parent_left)*x(locind_parent_left+1))*invds 108 endif 111 diff=(globind_parent_left + ds_parent - ypos)*invds 112 113 ! quick fix for roundoff error 114 diff=nint(diff*coeffraf)/real(coeffraf) 115 ! y(nc) = ((globind_parent_left + ds_parent - ypos)*x(locind_parent_left) & 116 ! + (ypos - globind_parent_left)*x(locind_parent_left+1))*invds 117 y(nc) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1)) 118 endif 119 109 120 !--------------------------------------------------------------------------------------------------- 110 121 end subroutine Agrif_basicinterp_linear1D … … 120 131 !--------------------------------------------------------------------------------------------------- 121 132 integer, intent(in) :: np,nc,np2 122 real , intent(in) :: s_parent, s_child123 real , intent(in) :: ds_parent, ds_child133 real(kind=8), intent(in) :: s_parent, s_child 134 real(kind=8), intent(in) :: ds_parent, ds_child 124 135 integer, intent(in) :: dir 125 136 ! … … 127 138 integer, dimension(:,:), allocatable :: indparent_tmp 128 139 real, dimension(:,:), allocatable :: coeffparent_tmp 129 real :: ypos,globind_parent_left,globind_parent_right130 real :: invds, invds2, invds3131 real :: ypos2,diff140 real(kind=8) :: ypos,globind_parent_left,globind_parent_right 141 real(kind=8) :: invds, invds2, invds3 142 real(kind=8) :: ypos2,diff 132 143 ! 133 144 coeffraf = nint(ds_parent/ds_child) … … 164 175 if (ypos2 > globind_parent_right) then 165 176 locind_parent_left = locind_parent_left + 1 166 globind_parent_right = globind_parent_right + 1. 177 globind_parent_right = globind_parent_right + 1.d0 167 178 ypos2 = ypos*invds+(i-1)*invds2 168 179 endif … … 239 250 real, dimension(np), intent(in) :: x 240 251 real, dimension(nc), intent(out) :: y 241 real , intent(in) :: s_parent, s_child242 real , intent(in) :: ds_parent, ds_child252 real(kind=8), intent(in) :: s_parent, s_child 253 real(kind=8), intent(in) :: ds_parent, ds_child 243 254 ! 244 255 integer :: i, coeffraf, locind_parent_left 245 real :: ypos,globind_parent_left246 real :: deltax, invdsparent256 real(kind=8) :: ypos,globind_parent_left 257 real(kind=8) :: deltax, invdsparent 247 258 real :: t2,t3,t4,t5,t6,t7,t8 248 259 ! … … 304 315 real, dimension(np), intent(in) :: x 305 316 real, dimension(nc), intent(out) :: y 306 real , intent(in) :: s_parent, s_child307 real , intent(in) :: ds_parent, ds_child317 real(kind=8), intent(in) :: s_parent, s_child 318 real(kind=8), intent(in) :: ds_parent, ds_child 308 319 ! 309 320 integer :: i, coeffraf, locind_parent 310 real :: ypos321 real(kind=8) :: ypos 311 322 ! 312 323 coeffraf = nint(ds_parent/ds_child) … … 342 353 real, dimension(np), intent(in) :: x 343 354 real, dimension(nc), intent(out) :: y 344 real , intent(in) :: s_parent, s_child345 real , intent(in) :: ds_parent, ds_child355 real(kind=8), intent(in) :: s_parent, s_child 356 real(kind=8), intent(in) :: ds_parent, ds_child 346 357 ! 347 358 real, dimension(:), allocatable :: ytemp 348 359 integer :: i,coeffraf,locind_parent_left,locind_parent_last 349 real :: ypos,xdiffmod,xpmin,xpmax,slope360 real(kind=8) :: ypos,xdiffmod,xpmin,xpmax,slope 350 361 integer :: i1,i2,ii 351 362 integer :: diffmod … … 429 440 real, dimension(np), intent(in) :: x 430 441 real, dimension(nc), intent(out) :: y 431 real , intent(in) :: s_parent, s_child432 real , intent(in) :: ds_parent, ds_child442 real(kind=8), intent(in) :: s_parent, s_child 443 real(kind=8), intent(in) :: ds_parent, ds_child 433 444 ! 434 445 real, dimension(:), allocatable :: ytemp 435 446 integer :: i,coeffraf,locind_parent_left,locind_parent_last 436 real :: ypos,xdiffmod,xpmin,xpmax,slope447 real(kind=8) :: ypos,xdiffmod,xpmin,xpmax,slope 437 448 integer :: i1,i2,ii 438 449 integer :: diffmod … … 524 535 real, dimension(np), intent(in) :: x 525 536 real, dimension(nc), intent(out) :: y 526 real , intent(in) :: s_parent, s_child527 real , intent(in) :: ds_parent, ds_child537 real(kind=8), intent(in) :: s_parent, s_child 538 real(kind=8), intent(in) :: ds_parent, ds_child 528 539 ! 529 540 integer :: i,coeffraf,locind_parent_left,locind_parent_last 530 541 integer :: iparent,ipos,pos,nmin,nmax 531 real :: ypos542 real(kind=8) :: ypos 532 543 integer :: i1,jj 533 real :: xpmin,a 544 real(kind=8) :: xpmin 545 real :: a 534 546 ! 535 547 real, dimension(np) :: xl,delta,a6,slope … … 646 658 !--------------------------------------------------------------------------------------------------- 647 659 integer, intent(in) :: np2, np, nc 648 real , intent(in) :: s_parent, s_child649 real , intent(in) :: ds_parent, ds_child660 real(kind=8), intent(in) :: s_parent, s_child 661 real(kind=8), intent(in) :: ds_parent, ds_child 650 662 integer, intent(in) :: dir 651 663 ! … … 655 667 integer :: iparent,ipos,pos 656 668 real :: ypos 657 integer :: i1,jj 658 real :: xpmin,a 669 integer :: i1,jj,k,l,j 670 real(kind=8) :: xpmin 671 real :: a 659 672 ! 660 673 integer :: diffmod … … 738 751 enddo 739 752 ! 740 do i = 1,np2 741 indparentppm(1+(i-1)*nc:i*nc,dir) = indparentppm_1d(1:nc) + (i-1)*np 742 indchildppm (1+(i-1)*nc:i*nc,dir) = indchildppm_1d (1:nc) 753 754 k=1 755 l=0 756 do i=1,np2 757 do j=1,nc 758 indchildppm(k,dir) = indchildppm_1d(j) 759 indparentppm(k,dir) = indparentppm_1d(j) + l 760 k=k+1 761 enddo 762 l=l+np 743 763 enddo 744 764 !--------------------------------------------------------------------------------------------------- … … 746 766 !=================================================================================================== 747 767 ! 748 !=================================================================================================== 749 !subroutine PPM1dPrecompute(np,nc,& 750 ! s_parent,s_child,ds_parent,ds_child) 751 !! 752 !!CC Description: 753 !!CC subroutine to compute coefficient and index for a 1D interpolation 754 !!CC using piecewise parabolic method 755 !!C Method: 756 !! 757 !! Declarations: 758 !! 759 ! Implicit none 760 !! 761 !! Arguments 762 ! Integer :: np,nc 763 !! Real, Dimension(:),Allocatable :: ytemp 764 ! Real :: s_parent,s_child,ds_parent,ds_child 765 !! 766 !! Local scalars 767 ! Integer :: i,coeffraf,locind_parent_left,locind_parent_last 768 ! Integer :: iparent,ipos,pos,nmin,nmax 769 ! Real :: ypos 770 ! integer :: i1,jj 771 ! Real :: xpmin,a 772 !! 773 ! Real :: xrmin,xrmax,am3,s2,s1 774 ! Real, Dimension(np) :: xl,delta,a6,slope 775 !! Real, Dimension(:),Allocatable :: diff,diff2,diff3 776 ! INTEGER :: diffmod 777 ! REAL :: invcoeffraf 778 !! 779 ! coeffraf = nint(ds_parent/ds_child) 780 !! 781 ! If (coeffraf == 1) Then 782 ! return 783 ! End If 784 ! invcoeffraf = ds_child/ds_parent 785 !! 786 ! 787 ! if (.not.allocated(indparentppm)) then 788 ! allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1),& 789 ! indchildppm(-2*coeffraf:nc+2*coeffraf,1)) 790 ! else 791 ! if (size(indparentppm,1)<nc+4*coeffraf+1) then 792 ! deallocate(indparentppm,indchildppm) 793 ! allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1),& 794 ! indchildppm(-2*coeffraf:nc+2*coeffraf,1)) 795 ! endif 796 ! endif 797 ! 798 ! ypos = s_child 799 !! 800 ! locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 801 ! locind_parent_last = 1 +& 802 ! agrif_ceiling((ypos +(nc - 1)& 803 ! *ds_child - s_parent)/ds_parent) 804 !! 805 ! xpmin = s_parent + (locind_parent_left-1)*ds_parent 806 ! i1 = 1+agrif_int((xpmin-s_child)/ds_child) 807 !! 808 !! 809 ! 810 ! Do i=1,coeffraf 811 ! tabdiff2(i)=(real(i)-0.5)*invcoeffraf 812 ! EndDo 813 ! 814 ! a = invcoeffraf**2 815 ! tabdiff3(1) = (1./3.)*a 816 ! a=2.*a 817 !!CDIR ALTCODE 818 !!!!CDIR SHORTLOOP 819 ! Do i=2,coeffraf 820 ! tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a 821 ! EndDo 822 ! 823 !!CDIR ALTCODE 824 !!!!CDIR SHORTLOOP 825 ! Do i=1,coeffraf 826 ! tabppm(1,i,1) = 0.08333333333333*(-1.+4*tabdiff2(i)-3*tabdiff3(i)) 827 ! tabppm(2,i,1) = 0.08333333333333*& 828 ! (7.-26.*tabdiff2(i)+18.*tabdiff3(i)) 829 ! tabppm(3,i,1) =0.08333333333333*(7.+30*tabdiff2(i)-30*tabdiff3(i)) 830 ! tabppm(4,i,1) = 0.08333333333333*& 831 ! (-1.-10.*tabdiff2(i)+18.*tabdiff3(i)) 832 ! tabppm(5,i,1) = 0.08333333333333*(2*tabdiff2(i)-3*tabdiff3(i)) 833 ! End Do 834 !! 835 !! 836 ! diffmod = 0 837 ! IF (mod(coeffraf,2) == 0) diffmod = 1 838 !! 839 ! ipos = i1 840 !! 841 ! Do iparent = locind_parent_left,locind_parent_last 842 ! pos=1 843 !!CDIR ALTCODE 844 !!CDIR SHORTLOOP 845 ! Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2 846 ! indparentppm(jj,1) = iparent-2 847 ! indchildppm(jj,1) = pos 848 ! pos = pos+1 849 ! End do 850 ! ipos = ipos + coeffraf 851 !! 852 ! End do 853 ! 854 ! Return 855 ! End subroutine ppm1dprecompute 768 856 769 !=================================================================================================== 857 770 ! … … 863 776 ! Use precomputed coefficient and index. 864 777 !--------------------------------------------------------------------------------------------------- 865 subroutine PPM1dAfterCompute ( x, y, np, nc, dir ) 866 !--------------------------------------------------------------------------------------------------- 778 subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc ) 779 !--------------------------------------------------------------------------------------------------- 780 integer, intent(in) :: np, nc 867 781 real, dimension(np), intent(in) :: x 868 782 real, dimension(nc), intent(out) :: y 869 integer, intent(in) :: np, nc870 783 integer, intent(in) :: dir 871 ! 872 integer :: i 873 ! 784 integer, dimension(1:),intent(in) :: indchildppmloc 785 real, dimension(1:,1:),intent(in) :: tabppmloc 786 ! 787 integer :: i,j,jp1,jp2,jp3,jp4,k 788 789 874 790 do i = 1,nc 875 y(i) = tabppm(1,indchildppm(i,dir),dir) * x(indparentppm(i,dir) ) + & 876 tabppm(2,indchildppm(i,dir),dir) * x(indparentppm(i,dir)+1) + & 877 tabppm(3,indchildppm(i,dir),dir) * x(indparentppm(i,dir)+2) + & 878 tabppm(4,indchildppm(i,dir),dir) * x(indparentppm(i,dir)+3) + & 879 tabppm(5,indchildppm(i,dir),dir) * x(indparentppm(i,dir)+4) 880 enddo 791 j = indparentppm(i,dir) 792 jp1=j+1 793 jp2=jp1+1 794 jp3=jp2+1 795 jp4=jp3+1 796 y(i) = tabppmloc(1,indchildppmloc(i)) * x(j) + & 797 tabppmloc(2,indchildppmloc(i)) * x(jp1) + & 798 tabppmloc(3,indchildppmloc(i)) * x(jp2) + & 799 tabppmloc(4,indchildppmloc(i)) * x(jp3) + & 800 tabppmloc(5,indchildppmloc(i)) * x(jp4) 801 enddo 802 881 803 !--------------------------------------------------------------------------------------------------- 882 804 end subroutine PPM1dAfterCompute 805 883 806 !=================================================================================================== 884 807 ! … … 1069 992 real, dimension(np), intent(in) :: x 1070 993 real, dimension(nc), intent(out) :: y 1071 real , intent(in) :: s_parent, s_child1072 real , intent(in) :: ds_parent, ds_child994 real(kind=8), intent(in) :: s_parent, s_child 995 real(kind=8), intent(in) :: ds_parent, ds_child 1073 996 ! 1074 997 real, dimension(:), allocatable :: ytemp 1075 998 integer :: i,coeffraf,locind_parent_left,locind_parent_last 1076 999 integer :: iparent,ipos,pos,nmin,nmax 1077 real :: ypos1000 real(kind=8) :: ypos 1078 1001 integer :: i1,jj 1079 real :: xpmin1002 real(kind=8) :: xpmin 1080 1003 ! 1081 1004 real, dimension(np) :: slope … … 1166 1089 real, dimension(np), intent(in) :: x 1167 1090 real, dimension(nc), intent(out) :: y 1168 real , intent(in) :: s_parent, s_child1169 real , intent(in) :: ds_parent, ds_child1091 real(kind=8), intent(in) :: s_parent, s_child 1092 real(kind=8), intent(in) :: ds_parent, ds_child 1170 1093 ! 1171 1094 integer :: i,coeffraf,locind_parent_left,locind_parent_last 1172 1095 integer :: ipos, pos 1173 real :: ypos,xi1096 real(kind=8) :: ypos,xi 1174 1097 integer :: i1,jj 1175 real :: xpmin1098 real(kind=8) :: xpmin 1176 1099 ! 1177 1100 real, dimension(:), allocatable :: ytemp … … 1276 1199 Real, Dimension(nc) :: y 1277 1200 Real, Dimension(:),Allocatable :: ytemp 1278 Real 1201 Real(kind=8) :: s_parent,s_child,ds_parent,ds_child 1279 1202 ! 1280 1203 ! Local scalars 1281 1204 Integer :: i,coeffraf,locind_parent_left,locind_parent_last 1282 1205 Integer :: iparent,ipos,pos,nmin,nmax 1283 Real :: ypos1206 Real(kind=8) :: ypos 1284 1207 integer :: i1,jj 1285 Real :: xpmin,cavg,a,b 1208 Real(kind=8) :: xpmin 1209 real :: cavg,a,b 1286 1210 ! 1287 1211 Real :: xrmin,xrmax,am3,s2,s1
Note: See TracChangeset
for help on using the changeset viewer.