Ignore:
Timestamp:
2019-02-27T14:55:54+01:00 (19 months ago)
Author:
rblod
Message:

Update agrif library and conv see ticket #2129

File:
1 edited

Legend:

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

    r10087 r10725  
    4141    integer, dimension(:), allocatable      :: indparentppm_1d, indchildppm_1d 
    4242! 
    43  
    4443    private :: Agrif_limiter_vanleer 
    4544! 
     
    5756    integer,             intent(in)     :: np           !< Length of input array 
    5857    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) 
    6362! 
    6463    integer :: i, coeffraf, locind_parent_left 
    65     real(kind=8)    :: globind_parent_left, globind_parent_right 
    66     real(kind=8)    :: invds, invds2, ypos, ypos2, diff 
     64    real    :: globind_parent_left, globind_parent_right 
     65    real    :: invds, invds2, ypos, ypos2, diff 
    6766! 
    6867    coeffraf = nint(ds_parent/ds_child) 
     
    9392! 
    9493        diff = globind_parent_right - ypos2 
    95 ! quick fix for roundoff error 
    96         diff=nint(diff*coeffraf)/real(coeffraf) 
    97  
    9894        y(i) = (diff*x(locind_parent_left) + (1.-diff)*x(locind_parent_left+1)) 
    99  
    10095        ypos2 = ypos2 + invds2 
    10196! 
     
    109104    else 
    110105        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 
    120109!--------------------------------------------------------------------------------------------------- 
    121110end subroutine Agrif_basicinterp_linear1D 
     
    131120!--------------------------------------------------------------------------------------------------- 
    132121    integer, intent(in) :: np,nc,np2 
    133     real(kind=8),    intent(in) :: s_parent, s_child 
    134     real(kind=8),    intent(in) :: ds_parent, ds_child 
     122    real,    intent(in) :: s_parent, s_child 
     123    real,    intent(in) :: ds_parent, ds_child 
    135124    integer, intent(in) :: dir 
    136125! 
     
    138127    integer, dimension(:,:), allocatable :: indparent_tmp 
    139128    real, dimension(:,:), allocatable :: coeffparent_tmp 
    140     real(kind=8)    :: ypos,globind_parent_left,globind_parent_right 
    141     real(kind=8)    :: invds, invds2, invds3 
    142     real(kind=8) :: ypos2,diff 
     129    real    :: ypos,globind_parent_left,globind_parent_right 
     130    real    :: invds, invds2, invds3 
     131    real :: ypos2,diff 
    143132! 
    144133    coeffraf = nint(ds_parent/ds_child) 
     
    175164        if (ypos2 > globind_parent_right) then 
    176165            locind_parent_left = locind_parent_left + 1 
    177             globind_parent_right = globind_parent_right + 1.d0 
     166            globind_parent_right = globind_parent_right + 1. 
    178167            ypos2 = ypos*invds+(i-1)*invds2 
    179168        endif 
     
    250239    real, dimension(np), intent(in)     :: x 
    251240    real, dimension(nc), intent(out)    :: y 
    252     real(kind=8),                intent(in)     :: s_parent, s_child 
    253     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     241    real,                intent(in)     :: s_parent, s_child 
     242    real,                intent(in)     :: ds_parent, ds_child 
    254243! 
    255244    integer :: i, coeffraf, locind_parent_left 
    256     real(kind=8)    :: ypos,globind_parent_left 
    257     real(kind=8)    :: deltax, invdsparent 
     245    real    :: ypos,globind_parent_left 
     246    real    :: deltax, invdsparent 
    258247    real    :: t2,t3,t4,t5,t6,t7,t8 
    259248! 
     
    315304    real, dimension(np), intent(in)     :: x 
    316305    real, dimension(nc), intent(out)    :: y 
    317     real(kind=8),                intent(in)     :: s_parent, s_child 
    318     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     306    real,                intent(in)     :: s_parent, s_child 
     307    real,                intent(in)     :: ds_parent, ds_child 
    319308! 
    320309    integer :: i, coeffraf, locind_parent 
    321     real(kind=8)    :: ypos 
     310    real    :: ypos 
    322311! 
    323312    coeffraf = nint(ds_parent/ds_child) 
     
    353342    real, dimension(np), intent(in)     :: x 
    354343    real, dimension(nc), intent(out)    :: y 
    355     real(kind=8),                intent(in)     :: s_parent, s_child 
    356     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     344    real,                intent(in)     :: s_parent, s_child 
     345    real,                intent(in)     :: ds_parent, ds_child 
    357346! 
    358347    real, dimension(:), allocatable :: ytemp 
    359348    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    360     real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope 
     349    real    :: ypos,xdiffmod,xpmin,xpmax,slope 
    361350    integer :: i1,i2,ii 
    362351    integer :: diffmod 
     
    440429    real, dimension(np), intent(in)     :: x 
    441430    real, dimension(nc), intent(out)    :: y 
    442     real(kind=8),                intent(in)     :: s_parent, s_child 
    443     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     431    real,                intent(in)     :: s_parent, s_child 
     432    real,                intent(in)     :: ds_parent, ds_child 
    444433! 
    445434    real, dimension(:), allocatable :: ytemp 
    446435    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    447     real(kind=8)    :: ypos,xdiffmod,xpmin,xpmax,slope 
     436    real    :: ypos,xdiffmod,xpmin,xpmax,slope 
    448437    integer :: i1,i2,ii 
    449438    integer :: diffmod 
     
    535524    real, dimension(np), intent(in)     :: x 
    536525    real, dimension(nc), intent(out)    :: y 
    537     real(kind=8),                intent(in)     :: s_parent, s_child 
    538     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     526    real,                intent(in)     :: s_parent, s_child 
     527    real,                intent(in)     :: ds_parent, ds_child 
    539528! 
    540529    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    541530    integer :: iparent,ipos,pos,nmin,nmax 
    542     real(kind=8)    :: ypos 
     531    real    :: ypos 
    543532    integer :: i1,jj 
    544     real(kind=8) :: xpmin 
    545     real :: a 
     533    real :: xpmin,a 
    546534! 
    547535    real, dimension(np) :: xl,delta,a6,slope 
     
    658646!--------------------------------------------------------------------------------------------------- 
    659647    integer,             intent(in)     :: np2, np, nc 
    660     real(kind=8),                intent(in)     :: s_parent, s_child 
    661     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     648    real,                intent(in)     :: s_parent, s_child 
     649    real,                intent(in)     :: ds_parent, ds_child 
    662650    integer,             intent(in)     :: dir 
    663651! 
     
    667655    integer :: iparent,ipos,pos 
    668656    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 
    672659! 
    673660    integer :: diffmod 
     
    751738    enddo 
    752739! 
    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) 
    763743    enddo 
    764744!--------------------------------------------------------------------------------------------------- 
     
    766746!=================================================================================================== 
    767747! 
    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 
    769856!=================================================================================================== 
    770857! 
     
    776863! Use precomputed coefficient and index. 
    777864!--------------------------------------------------------------------------------------------------- 
    778 subroutine PPM1dAfterCompute ( x, y, np, nc, dir, indchildppmloc, tabppmloc ) 
    779 !--------------------------------------------------------------------------------------------------- 
    780     integer,             intent(in)     :: np, nc 
     865subroutine PPM1dAfterCompute ( x, y, np, nc, dir ) 
     866!--------------------------------------------------------------------------------------------------- 
    781867    real, dimension(np), intent(in)     :: x 
    782868    real, dimension(nc), intent(out)    :: y 
     869    integer,             intent(in)     :: np, nc 
    783870    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! 
    790874    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 
    803881!--------------------------------------------------------------------------------------------------- 
    804882end subroutine PPM1dAfterCompute 
    805  
    806883!=================================================================================================== 
    807884! 
     
    9921069    real, dimension(np), intent(in)     :: x 
    9931070    real, dimension(nc), intent(out)    :: y 
    994     real(kind=8),                intent(in)     :: s_parent, s_child 
    995     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     1071    real,                intent(in)     :: s_parent, s_child 
     1072    real,                intent(in)     :: ds_parent, ds_child 
    9961073! 
    9971074    real, dimension(:), allocatable :: ytemp 
    9981075    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    9991076    integer :: iparent,ipos,pos,nmin,nmax 
    1000     real(kind=8)    :: ypos 
     1077    real    :: ypos 
    10011078    integer :: i1,jj 
    1002     real(kind=8) :: xpmin 
     1079    real :: xpmin 
    10031080! 
    10041081    real, dimension(np) :: slope 
     
    10891166    real, dimension(np), intent(in)     :: x 
    10901167    real, dimension(nc), intent(out)    :: y 
    1091     real(kind=8),                intent(in)     :: s_parent, s_child 
    1092     real(kind=8),                intent(in)     :: ds_parent, ds_child 
     1168    real,                intent(in)     :: s_parent, s_child 
     1169    real,                intent(in)     :: ds_parent, ds_child 
    10931170! 
    10941171    integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    10951172    integer :: ipos, pos 
    1096     real(kind=8)    :: ypos,xi 
     1173    real    :: ypos,xi 
    10971174    integer :: i1,jj 
    1098     real(kind=8) :: xpmin 
     1175    real :: xpmin 
    10991176! 
    11001177    real, dimension(:),   allocatable  :: ytemp 
     
    11991276      Real, Dimension(nc) :: y 
    12001277      Real, Dimension(:),Allocatable :: ytemp 
    1201       Real(kind=8)        :: s_parent,s_child,ds_parent,ds_child 
     1278      Real                :: s_parent,s_child,ds_parent,ds_child 
    12021279! 
    12031280!     Local scalars 
    12041281      Integer :: i,coeffraf,locind_parent_left,locind_parent_last 
    12051282      Integer :: iparent,ipos,pos,nmin,nmax 
    1206       Real(kind=8)    :: ypos 
     1283      Real    :: ypos 
    12071284      integer :: i1,jj 
    1208       Real(kind=8) :: xpmin 
    1209       real :: cavg,a,b 
     1285      Real :: xpmin,cavg,a,b 
    12101286!       
    12111287      Real :: xrmin,xrmax,am3,s2,s1   
Note: See TracChangeset for help on using the changeset viewer.