Ignore:
Timestamp:
2018-09-05T15:33:44+02:00 (2 years ago)
Author:
rblod
Message:

update AGRIF library

File:
1 edited

Legend:

Unmodified
Added
Removed
  • vendors/AGRIF/CMEMS_2020/AGRIF_FILES/modinterpbasic.F90

    r5656 r10087  
    4141    integer, dimension(:), allocatable      :: indparentppm_1d, indchildppm_1d 
    4242! 
     43 
    4344    private :: Agrif_limiter_vanleer 
    4445! 
     
    5657    integer,             intent(in)     :: np           !< Length of input array 
    5758    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) 
    6263! 
    6364    integer :: i, coeffraf, locind_parent_left 
    64     real    :: globind_parent_left, globind_parent_right 
    65     real    :: invds, invds2, ypos, ypos2, diff 
     65    real(kind=8)    :: globind_parent_left, globind_parent_right 
     66    real(kind=8)    :: invds, invds2, ypos, ypos2, diff 
    6667! 
    6768    coeffraf = nint(ds_parent/ds_child) 
     
    9293! 
    9394        diff = globind_parent_right - ypos2 
     95! quick fix for roundoff error 
     96        diff=nint(diff*coeffraf)/real(coeffraf) 
     97 
    9498        y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1)) 
     99 
    95100        ypos2 = ypos2 + invds2 
    96101! 
     
    104109    else 
    105110        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 
    109120!--------------------------------------------------------------------------------------------------- 
    110121end subroutine Agrif_basicinterp_linear1D 
     
    120131!--------------------------------------------------------------------------------------------------- 
    121132    integer, intent(in) :: np,nc,np2 
    122     real,    intent(in) :: s_parent, s_child 
    123     real,    intent(in) :: ds_parent, ds_child 
     133    real(kind=8),    intent(in) :: s_parent, s_child 
     134    real(kind=8),    intent(in) :: ds_parent, ds_child 
    124135    integer, intent(in) :: dir 
    125136! 
     
    127138    integer, dimension(:,:), allocatable :: indparent_tmp 
    128139    real, dimension(:,:), allocatable :: coeffparent_tmp 
    129     real    :: ypos,globind_parent_left,globind_parent_right 
    130     real    :: invds, invds2, invds3 
    131     real :: ypos2,diff 
     140    real(kind=8)    :: ypos,globind_parent_left,globind_parent_right 
     141    real(kind=8)    :: invds, invds2, invds3 
     142    real(kind=8) :: ypos2,diff 
    132143! 
    133144    coeffraf = nint(ds_parent/ds_child) 
     
    164175        if (ypos2 > globind_parent_right) then 
    165176            locind_parent_left = locind_parent_left + 1 
    166             globind_parent_right = globind_parent_right + 1. 
     177            globind_parent_right = globind_parent_right + 1.d0 
    167178            ypos2 = ypos*invds+(i-1)*invds2 
    168179        endif 
     
    239250    real, dimension(np), intent(in)     :: x 
    240251    real, dimension(nc), intent(out)    :: y 
    241     real,                intent(in)     :: s_parent, s_child 
    242     real,                intent(in)     :: ds_parent, ds_child 
     252    real(kind=8),                intent(in)     :: s_parent, s_child 
     253    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    243254! 
    244255    integer :: i, coeffraf, locind_parent_left 
    245     real    :: ypos,globind_parent_left 
    246     real    :: deltax, invdsparent 
     256    real(kind=8)    :: ypos,globind_parent_left 
     257    real(kind=8)    :: deltax, invdsparent 
    247258    real    :: t2,t3,t4,t5,t6,t7,t8 
    248259! 
     
    304315    real, dimension(np), intent(in)     :: x 
    305316    real, dimension(nc), intent(out)    :: y 
    306     real,                intent(in)     :: s_parent, s_child 
    307     real,                intent(in)     :: ds_parent, ds_child 
     317    real(kind=8),                intent(in)     :: s_parent, s_child 
     318    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    308319! 
    309320    integer :: i, coeffraf, locind_parent 
    310     real    :: ypos 
     321    real(kind=8)    :: ypos 
    311322! 
    312323    coeffraf = nint(ds_parent/ds_child) 
     
    342353    real, dimension(np), intent(in)     :: x 
    343354    real, dimension(nc), intent(out)    :: y 
    344     real,                intent(in)     :: s_parent, s_child 
    345     real,                intent(in)     :: ds_parent, ds_child 
     355    real(kind=8),                intent(in)     :: s_parent, s_child 
     356    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    346357! 
    347358    real, dimension(:), allocatable :: ytemp 
    348359    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    349     real    :: ypos,xdiffmod,xpmin,xpmax,slope 
     360    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope 
    350361    integer :: i1,i2,ii 
    351362    integer :: diffmod 
     
    429440    real, dimension(np), intent(in)     :: x 
    430441    real, dimension(nc), intent(out)    :: y 
    431     real,                intent(in)     :: s_parent, s_child 
    432     real,                intent(in)     :: ds_parent, ds_child 
     442    real(kind=8),                intent(in)     :: s_parent, s_child 
     443    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    433444! 
    434445    real, dimension(:), allocatable :: ytemp 
    435446    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    436     real    :: ypos,xdiffmod,xpmin,xpmax,slope 
     447    real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope 
    437448    integer :: i1,i2,ii 
    438449    integer :: diffmod 
     
    524535    real, dimension(np), intent(in)     :: x 
    525536    real, dimension(nc), intent(out)    :: y 
    526     real,                intent(in)     :: s_parent, s_child 
    527     real,                intent(in)     :: ds_parent, ds_child 
     537    real(kind=8),                intent(in)     :: s_parent, s_child 
     538    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    528539! 
    529540    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    530541    integer :: iparent,ipos,pos,nmin,nmax 
    531     real    :: ypos 
     542    real(kind=8)    :: ypos 
    532543    integer :: i1,jj 
    533     real :: xpmin,a 
     544    real(kind=8) :: xpmin 
     545    real :: a 
    534546! 
    535547    real, dimension(np) :: xl,delta,a6,slope 
     
    646658!--------------------------------------------------------------------------------------------------- 
    647659    integer,             intent(in)     :: np2, np, nc 
    648     real,                intent(in)     :: s_parent, s_child 
    649     real,                intent(in)     :: ds_parent, ds_child 
     660    real(kind=8),                intent(in)     :: s_parent, s_child 
     661    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    650662    integer,             intent(in)     :: dir 
    651663! 
     
    655667    integer :: iparent,ipos,pos 
    656668    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 
    659672! 
    660673    integer :: diffmod 
     
    738751    enddo 
    739752! 
    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 
    743763    enddo 
    744764!--------------------------------------------------------------------------------------------------- 
     
    746766!=================================================================================================== 
    747767! 
    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 
    856769!=================================================================================================== 
    857770! 
     
    863776! Use precomputed coefficient and index. 
    864777!--------------------------------------------------------------------------------------------------- 
    865 subroutine PPM1dAfterCompute ( x, y, np, nc, dir ) 
    866 !--------------------------------------------------------------------------------------------------- 
     778subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc ) 
     779!--------------------------------------------------------------------------------------------------- 
     780    integer,             intent(in)     :: np, nc 
    867781    real, dimension(np), intent(in)     :: x 
    868782    real, dimension(nc), intent(out)    :: y 
    869     integer,             intent(in)     :: np, nc 
    870783    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     
    874790    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 
    881803!--------------------------------------------------------------------------------------------------- 
    882804end subroutine PPM1dAfterCompute 
     805 
    883806!=================================================================================================== 
    884807! 
     
    1069992    real, dimension(np), intent(in)     :: x 
    1070993    real, dimension(nc), intent(out)    :: y 
    1071     real,                intent(in)     :: s_parent, s_child 
    1072     real,                intent(in)     :: ds_parent, ds_child 
     994    real(kind=8),                intent(in)     :: s_parent, s_child 
     995    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    1073996! 
    1074997    real, dimension(:), allocatable :: ytemp 
    1075998    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    1076999    integer :: iparent,ipos,pos,nmin,nmax 
    1077     real    :: ypos 
     1000    real(kind=8)    :: ypos 
    10781001    integer :: i1,jj 
    1079     real :: xpmin 
     1002    real(kind=8) :: xpmin 
    10801003! 
    10811004    real, dimension(np) :: slope 
     
    11661089    real, dimension(np), intent(in)     :: x 
    11671090    real, dimension(nc), intent(out)    :: y 
    1168     real,                intent(in)     :: s_parent, s_child 
    1169     real,                intent(in)     :: ds_parent, ds_child 
     1091    real(kind=8),                intent(in)     :: s_parent, s_child 
     1092    real(kind=8),                intent(in)     :: ds_parent, ds_child 
    11701093! 
    11711094    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    11721095    integer :: ipos, pos 
    1173     real    :: ypos,xi 
     1096    real(kind=8)    :: ypos,xi 
    11741097    integer :: i1,jj 
    1175     real :: xpmin 
     1098    real(kind=8) :: xpmin 
    11761099! 
    11771100    real, dimension(:),   allocatable  :: ytemp 
     
    12761199      Real, Dimension(nc) :: y 
    12771200      Real, Dimension(:),Allocatable :: ytemp 
    1278       Real                :: s_parent,s_child,ds_parent,ds_child 
     1201      Real(kind=8)        :: s_parent,s_child,ds_parent,ds_child 
    12791202! 
    12801203!     Local scalars 
    12811204      Integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    12821205      Integer :: iparent,ipos,pos,nmin,nmax 
    1283       Real    :: ypos 
     1206      Real(kind=8)    :: ypos 
    12841207      integer :: i1,jj 
    1285       Real :: xpmin,cavg,a,b 
     1208      Real(kind=8) :: xpmin 
     1209      real :: cavg,a,b 
    12861210!       
    12871211      Real :: xrmin,xrmax,am3,s2,s1   
Note: See TracChangeset for help on using the changeset viewer.