Changeset 10725 for vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90
- Timestamp:
- 2019-02-27T14:55:54+01:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90
r10087 r10725 41 41 integer, dimension(:), allocatable :: indparentppm_1d, indchildppm_1d 42 42 ! 43 44 43 private :: Agrif_limiter_vanleer 45 44 ! … … 57 56 integer, intent(in) :: np !< Length of input array 58 57 integer, intent(in) :: nc !< Length of output array 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)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) 63 62 ! 64 63 integer :: i, coeffraf, locind_parent_left 65 real (kind=8):: globind_parent_left, globind_parent_right66 real (kind=8):: invds, invds2, ypos, ypos2, diff64 real :: globind_parent_left, globind_parent_right 65 real :: invds, invds2, ypos, ypos2, diff 67 66 ! 68 67 coeffraf = nint(ds_parent/ds_child) … … 93 92 ! 94 93 diff = globind_parent_right - ypos2 95 ! quick fix for roundoff error96 diff=nint(diff*coeffraf)/real(coeffraf)97 98 94 y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1)) 99 100 95 ypos2 = ypos2 + invds2 101 96 ! … … 109 104 else 110 105 globind_parent_left = s_parent + (locind_parent_left - 1)*ds_parent 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 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 120 109 !--------------------------------------------------------------------------------------------------- 121 110 end subroutine Agrif_basicinterp_linear1D … … 131 120 !--------------------------------------------------------------------------------------------------- 132 121 integer, intent(in) :: np,nc,np2 133 real (kind=8), intent(in) :: s_parent, s_child134 real (kind=8), intent(in) :: ds_parent, ds_child122 real, intent(in) :: s_parent, s_child 123 real, intent(in) :: ds_parent, ds_child 135 124 integer, intent(in) :: dir 136 125 ! … … 138 127 integer, dimension(:,:), allocatable :: indparent_tmp 139 128 real, dimension(:,:), allocatable :: coeffparent_tmp 140 real (kind=8):: ypos,globind_parent_left,globind_parent_right141 real (kind=8):: invds, invds2, invds3142 real (kind=8):: ypos2,diff129 real :: ypos,globind_parent_left,globind_parent_right 130 real :: invds, invds2, invds3 131 real :: ypos2,diff 143 132 ! 144 133 coeffraf = nint(ds_parent/ds_child) … … 175 164 if (ypos2 > globind_parent_right) then 176 165 locind_parent_left = locind_parent_left + 1 177 globind_parent_right = globind_parent_right + 1. d0166 globind_parent_right = globind_parent_right + 1. 178 167 ypos2 = ypos*invds+(i-1)*invds2 179 168 endif … … 250 239 real, dimension(np), intent(in) :: x 251 240 real, dimension(nc), intent(out) :: y 252 real (kind=8), intent(in) :: s_parent, s_child253 real (kind=8), intent(in) :: ds_parent, ds_child241 real, intent(in) :: s_parent, s_child 242 real, intent(in) :: ds_parent, ds_child 254 243 ! 255 244 integer :: i, coeffraf, locind_parent_left 256 real (kind=8):: ypos,globind_parent_left257 real (kind=8):: deltax, invdsparent245 real :: ypos,globind_parent_left 246 real :: deltax, invdsparent 258 247 real :: t2,t3,t4,t5,t6,t7,t8 259 248 ! … … 315 304 real, dimension(np), intent(in) :: x 316 305 real, dimension(nc), intent(out) :: y 317 real (kind=8), intent(in) :: s_parent, s_child318 real (kind=8), intent(in) :: ds_parent, ds_child306 real, intent(in) :: s_parent, s_child 307 real, intent(in) :: ds_parent, ds_child 319 308 ! 320 309 integer :: i, coeffraf, locind_parent 321 real (kind=8):: ypos310 real :: ypos 322 311 ! 323 312 coeffraf = nint(ds_parent/ds_child) … … 353 342 real, dimension(np), intent(in) :: x 354 343 real, dimension(nc), intent(out) :: y 355 real (kind=8), intent(in) :: s_parent, s_child356 real (kind=8), intent(in) :: ds_parent, ds_child344 real, intent(in) :: s_parent, s_child 345 real, intent(in) :: ds_parent, ds_child 357 346 ! 358 347 real, dimension(:), allocatable :: ytemp 359 348 integer :: i,coeffraf,locind_parent_left,locind_parent_last 360 real (kind=8):: ypos,xdiffmod,xpmin,xpmax,slope349 real :: ypos,xdiffmod,xpmin,xpmax,slope 361 350 integer :: i1,i2,ii 362 351 integer :: diffmod … … 440 429 real, dimension(np), intent(in) :: x 441 430 real, dimension(nc), intent(out) :: y 442 real (kind=8), intent(in) :: s_parent, s_child443 real (kind=8), intent(in) :: ds_parent, ds_child431 real, intent(in) :: s_parent, s_child 432 real, intent(in) :: ds_parent, ds_child 444 433 ! 445 434 real, dimension(:), allocatable :: ytemp 446 435 integer :: i,coeffraf,locind_parent_left,locind_parent_last 447 real (kind=8):: ypos,xdiffmod,xpmin,xpmax,slope436 real :: ypos,xdiffmod,xpmin,xpmax,slope 448 437 integer :: i1,i2,ii 449 438 integer :: diffmod … … 535 524 real, dimension(np), intent(in) :: x 536 525 real, dimension(nc), intent(out) :: y 537 real (kind=8), intent(in) :: s_parent, s_child538 real (kind=8), intent(in) :: ds_parent, ds_child526 real, intent(in) :: s_parent, s_child 527 real, intent(in) :: ds_parent, ds_child 539 528 ! 540 529 integer :: i,coeffraf,locind_parent_left,locind_parent_last 541 530 integer :: iparent,ipos,pos,nmin,nmax 542 real (kind=8):: ypos531 real :: ypos 543 532 integer :: i1,jj 544 real(kind=8) :: xpmin 545 real :: a 533 real :: xpmin,a 546 534 ! 547 535 real, dimension(np) :: xl,delta,a6,slope … … 658 646 !--------------------------------------------------------------------------------------------------- 659 647 integer, intent(in) :: np2, np, nc 660 real (kind=8), intent(in) :: s_parent, s_child661 real (kind=8), intent(in) :: ds_parent, ds_child648 real, intent(in) :: s_parent, s_child 649 real, intent(in) :: ds_parent, ds_child 662 650 integer, intent(in) :: dir 663 651 ! … … 667 655 integer :: iparent,ipos,pos 668 656 real :: ypos 669 integer :: i1,jj,k,l,j 670 real(kind=8) :: xpmin 671 real :: a 657 integer :: i1,jj 658 real :: xpmin,a 672 659 ! 673 660 integer :: diffmod … … 751 738 enddo 752 739 ! 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 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) 763 743 enddo 764 744 !--------------------------------------------------------------------------------------------------- … … 766 746 !=================================================================================================== 767 747 ! 768 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 769 856 !=================================================================================================== 770 857 ! … … 776 863 ! Use precomputed coefficient and index. 777 864 !--------------------------------------------------------------------------------------------------- 778 subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc ) 779 !--------------------------------------------------------------------------------------------------- 780 integer, intent(in) :: np, nc 865 subroutine PPM1dAfterCompute ( x, y, np, nc, dir ) 866 !--------------------------------------------------------------------------------------------------- 781 867 real, dimension(np), intent(in) :: x 782 868 real, dimension(nc), intent(out) :: y 869 integer, intent(in) :: np, nc 783 870 integer, intent(in) :: dir 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 871 ! 872 integer :: i 873 ! 790 874 do i = 1,nc 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 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 803 881 !--------------------------------------------------------------------------------------------------- 804 882 end subroutine PPM1dAfterCompute 805 806 883 !=================================================================================================== 807 884 ! … … 992 1069 real, dimension(np), intent(in) :: x 993 1070 real, dimension(nc), intent(out) :: y 994 real (kind=8), intent(in) :: s_parent, s_child995 real (kind=8), intent(in) :: ds_parent, ds_child1071 real, intent(in) :: s_parent, s_child 1072 real, intent(in) :: ds_parent, ds_child 996 1073 ! 997 1074 real, dimension(:), allocatable :: ytemp 998 1075 integer :: i,coeffraf,locind_parent_left,locind_parent_last 999 1076 integer :: iparent,ipos,pos,nmin,nmax 1000 real (kind=8):: ypos1077 real :: ypos 1001 1078 integer :: i1,jj 1002 real (kind=8):: xpmin1079 real :: xpmin 1003 1080 ! 1004 1081 real, dimension(np) :: slope … … 1089 1166 real, dimension(np), intent(in) :: x 1090 1167 real, dimension(nc), intent(out) :: y 1091 real (kind=8), intent(in) :: s_parent, s_child1092 real (kind=8), intent(in) :: ds_parent, ds_child1168 real, intent(in) :: s_parent, s_child 1169 real, intent(in) :: ds_parent, ds_child 1093 1170 ! 1094 1171 integer :: i,coeffraf,locind_parent_left,locind_parent_last 1095 1172 integer :: ipos, pos 1096 real (kind=8):: ypos,xi1173 real :: ypos,xi 1097 1174 integer :: i1,jj 1098 real (kind=8):: xpmin1175 real :: xpmin 1099 1176 ! 1100 1177 real, dimension(:), allocatable :: ytemp … … 1199 1276 Real, Dimension(nc) :: y 1200 1277 Real, Dimension(:),Allocatable :: ytemp 1201 Real (kind=8):: s_parent,s_child,ds_parent,ds_child1278 Real :: s_parent,s_child,ds_parent,ds_child 1202 1279 ! 1203 1280 ! Local scalars 1204 1281 Integer :: i,coeffraf,locind_parent_left,locind_parent_last 1205 1282 Integer :: iparent,ipos,pos,nmin,nmax 1206 Real (kind=8):: ypos1283 Real :: ypos 1207 1284 integer :: i1,jj 1208 Real(kind=8) :: xpmin 1209 real :: cavg,a,b 1285 Real :: xpmin,cavg,a,b 1210 1286 ! 1211 1287 Real :: xrmin,xrmax,am3,s2,s1
Note: See TracChangeset
for help on using the changeset viewer.